Commits

Olivier Debeir committed 68a342a

initial commit

Comments (0)

Files changed (45)

+.idea/
+data/
+*.pyc
+*.zip

dependencies/pymic/__init__.py

+__all__ = ['chain','meanshift','rankfilter','snakes','tools','watershed','zvi','color','utils']

dependencies/pymic/chain/__init__.py

+from chain import chain

dependencies/pymic/chain/chain.c

+/*CHAIN.CPP*/
+
+void zero_borders(char * str,char len,int sizex,int sizey,int ext);
+int build_chain(int * str,int sizex,int sizey, int label,int *str1,int *str2);
+int build_chain4ex(int * str,int sizex,int sizey, int label,int *str1,int *str2);
+
+/*Conventions de directions...*/
+
+int ChainBufM[CHAINMAXLENGTH];
+int ChainBufN[CHAINMAXLENGTH];
+int Dy[16] = {-1, -1 , 0 ,+1 ,+1, +1, 0 ,-1,-1, -1 , 0 ,+1 ,+1, +1, 0 ,-1};
+int Dx[16] = {0, +1, +1, +1, 0, -1, -1, -1, 0, +1, +1, +1, 0, -1, -1, -1};
+int	Dx4[8]={ 0,-1,-1,-1, 0, 1, 1, 1};
+int	Dy4[8]={ 1, 1, 0,-1,-1,-1, 0, 1};
+
+/*
+fonction qui calcule une chaine8 sur l'image pi et sauve le resultat dans 'chain'
+retourne une structure Matlab contenant les [x,y] des points de la chaine
+de l'objet Label de l'image pi. La recherche est faite e partir du premier pixel == Label
+LabelImage est un plan de int
+label est un double
+ */
+
+
+void zero_borders(char * str,char len,int sizex,int sizey,int ext) {
+	//set borders pixels to 0 to avoid craches
+	int *p;
+	int i;
+
+	p = ( int *)str;
+
+	for(i=0;i<sizey;i++) {
+		*(p+i) = 0;
+		*(p+i+sizey) = 0;
+		*(p+(sizex-1)*sizey) = 0;
+		*(p+(sizex-2)*sizey) = 0;
+		if(ext){
+			*(p+i+2*sizey) = 0;
+			*(p+(sizex-3)*sizey) = 0;
+		}
+	}
+	for(i=0;i<sizex;i++) {
+		*(p+i*sizey) = 0;
+		*(p+i*sizey+1) = 0;
+		*(p+(i+1)*sizey-1) = 0;
+		*(p+(i+1)*sizey-2) = 0;
+		if(ext){
+			*(p+i*sizey+2) = 0;
+			*(p+(i+1)*sizey-3) = 0;
+		}
+	}
+}
+
+int build_chain(int * str,int sizex,int sizey, int label,int *str1,int *str2)
+{
+	int m,n,i,j,done;
+	int x,y,x0,y0,x1,y1,prev_dir,new_dir,d;
+	long int l;
+	int *pint;
+	int *X,*Y;
+
+	m = sizex;
+	n = sizey;
+	pint = ( int*)str;
+
+	x0=-1;
+	for(l=0;l<m*n;l++)
+	{
+		if(*pint==label)
+		{
+			x0=l/m;
+			y0=l%m;
+			break;
+		}
+		pint++;
+	}
+	if(x0<0)// pas d'objet de ce label
+	{
+		return 0;
+	}
+
+	pint = ( int*)str;
+
+	prev_dir = 3;
+	x=x0;
+	y=y0;
+	i = 0;
+	while(1)
+	{
+		done=0;
+		ChainBufM[i] = x+1;
+		ChainBufN[i] = y+1;
+		i++;
+		new_dir = (prev_dir+5)%8;
+		for(d = new_dir;d<new_dir+8;d++)
+		{
+			x1=x+Dx[d];
+			y1=y+Dy[d];
+			if(*(pint+x1*m+y1)==label)
+			{
+				x=x1;y=y1;
+				prev_dir = d%8;
+				done=1;
+				break;
+			}
+		}
+		if((i>=CHAINMAXLENGTH)||((x==x0)&&(y==y0))||(done==0)) break;
+	}
+
+	X = (int*)str1;
+	Y = (int*)str2;
+	for(j=0;j<i;j++)
+	{
+		*(Y+j) = (int)ChainBufM[j];
+		*(X+j) = (int)ChainBufN[j];
+	}
+	return i;
+
+};
+
+
+//_____________________________________________________________OBJECT
+int build_chain4ex(int * str,int sizex,int sizey, int label,int *str1,int *str2)
+//fonction qui calcule une chaine8 sur l'image pi et sauve le résultat dans 'chain'
+//retourne une structure Matlab contenant les [x,y] des points de la chaine
+//de l'objet Label de l'image pi. La recherche est faite à partir du premier pixel == Label
+//LabelImage est un plan de double
+//label est un double
+//le contour construit est extérieur et de connexité 4
+{
+	int m,n,i,j,dir;
+	int x,y,x0,y0;
+	long int l;
+	int *pint;
+	int *X,*Y;
+
+	m = sizex;
+	n = sizey;
+	pint = (int*)str;
+
+	//recherche du premier pixel
+	x0=-1;
+	for(l=0;l<m*n;l++)
+	{
+		if(*pint==label)
+		{
+			x0=l/m;
+			y0=l%m;
+			break;
+		}
+		pint++;
+	}
+	x0 = x0-1;
+	y0 = y0;
+	if(x0<0) //pas d'objet de ce label
+	{
+		return 0;
+	}
+
+	pint = ( int*)str;
+
+	dir = 2;
+	x = x0;
+	y = y0;
+	i = 0;
+	while(1)
+	{
+		switch(dir){
+		case 0:
+			if((*(pint+(x+Dx[2])*m+(y+Dy[2]))==label)){x=x+Dx[4];y=y+Dy[4];dir=4;break;}
+			if((*(pint+(x+Dx[1])*m+(y+Dy[1]))==label)){x=x+Dx[2];y=y+Dy[2];dir=2;break;}
+			if((*(pint+(x+Dx[0])*m+(y+Dy[0]))==label)){x=x+Dx[2];y=y+Dy[2];dir=2;break;}
+			if((*(pint+(x+Dx[6])*m+(y+Dy[6]))==label)){x=x+Dx[0];y=y+Dy[0];dir=0;break;}
+			x=x+Dx[6];y=y+Dy[6];dir=6;
+			break;
+		case 2:
+			if((*(pint+(x+Dx[4])*m+(y+Dy[4]))==label)){x=x+Dx[6];y=y+Dy[6];dir=6;break;}
+			if((*(pint+(x+Dx[3])*m+(y+Dy[3]))==label)){x=x+Dx[4];y=y+Dy[4];dir=4;break;}
+			if((*(pint+(x+Dx[2])*m+(y+Dy[2]))==label)){x=x+Dx[4];y=y+Dy[4];dir=4;break;}
+			if((*(pint+(x+Dx[0])*m+(y+Dy[0]))==label)){x=x+Dx[2];y=y+Dy[2];dir=2;break;}
+			x=x+Dx[0];y=y+Dy[0];dir=0;
+			break;
+		case 4:
+			if((*(pint+(x+Dx[6])*m+(y+Dy[6]))==label)){x=x+Dx[0];y=y+Dy[0];dir=0;break;}
+			if((*(pint+(x+Dx[5])*m+(y+Dy[5]))==label)){x=x+Dx[6];y=y+Dy[6];dir=6;break;}
+			if((*(pint+(x+Dx[4])*m+(y+Dy[4]))==label)){x=x+Dx[6];y=y+Dy[6];dir=6;break;}
+			if((*(pint+(x+Dx[2])*m+(y+Dy[2]))==label)){x=x+Dx[4];y=y+Dy[4];dir=4;break;}
+			x=x+Dx[2];y=y+Dy[2];dir=2;
+			break;
+		case 6:
+			if((*(pint+(x+Dx[0])*m+(y+Dy[0]))==label)){x=x+Dx[2];y=y+Dy[2];dir=2;break;}
+			if((*(pint+(x+Dx[7])*m+(y+Dy[7]))==label)){x=x+Dx[0];y=y+Dy[0];dir=0;break;}
+			if((*(pint+(x+Dx[6])*m+(y+Dy[6]))==label)){x=x+Dx[0];y=y+Dy[0];dir=0;break;}
+			if((*(pint+(x+Dx[4])*m+(y+Dy[4]))==label)){x=x+Dx[6];y=y+Dy[6];dir=6;break;}
+			x=x+Dx[4];y=y+Dy[4];dir=4;
+			break;
+		}
+		ChainBufM[i] = x+1;
+		ChainBufN[i] = y+1;
+		i++;
+		if((i>=CHAINMAXLENGTH)||((x==x0)&&(y==y0))) break;
+	}
+
+	X = (int*)str1;
+	Y = (int*)str2;
+
+	for(j=0;j<i;j++)
+	{
+		*(Y+j) = (int)ChainBufM[j];
+		*(X+j) = (int)ChainBufN[j];
+	}
+	//return chain;
+	return i;
+
+};
+

dependencies/pymic/chain/chain.py

+# -*- coding: utf-8 -*-
+'''
+Test the scipy weave module to the chain function
+use Weave for C compilation
+'''
+
+#Stuff related to the inline compilation required to return a value with weave.inline
+import sys, os
+
+import pylab
+import numpy as npy
+from scipy.weave import inline
+import scipy
+
+from scipy.ndimage import measurements
+
+def dtype2ctype(array):
+    """convert numpy type in C equivalent type
+    """
+    types = {
+        npy.dtype(npy.float64): 'double',
+        npy.dtype(npy.float32): 'float',
+        npy.dtype(npy.int32):   'int',
+        npy.dtype(npy.uint32):  'unsigned int',
+        npy.dtype(npy.int16):   'short',
+        npy.dtype(npy.uint8):   'unsigned char',
+        npy.dtype(npy.uint16):  'unsigned short',
+    }
+    return types.get(array.dtype)
+
+watershedMode = ['unmarked','marked','modifiedGradient']
+                   
+def chain(ima,labels = None,ext=0,verbose = False):
+    """compute the chain of the label image im
+
+    :param ima: image array
+    :type ima: uint32 numpy.ndarray
+    :param labels: list of the labels to process, if None, all the present labels are processed
+    :type labels: list
+    :param ext: if 1 external border is computed
+    :type ext: int[0-1]
+    :param verbose: True : displays details during C code executes (default is False)
+    :type verbose: Bool
+
+    :returns:  chain for each label
+    :raises: TypeError
+
+    example:
+
+    >>> im = scipy.misc.imread('../../test/data/blocs.png')
+    >>> m = measurements.label(im)
+    >>> label = m[0]
+    >>> contours = chain(label,labels = None,ext = 0)
+    """
+
+    if not isinstance(ima,npy.ndarray):
+        raise TypeError('2D numpy.array expected')
+    if not (ima.dtype == npy.int32) :
+        raise TypeError('int32 numpy.array expected')
+    if not(len(ima.shape) == 2):
+        raise TypeError('2D numpy.array expected')
+    if labels is None:
+        labels = range(1,npy.max(ima)+1)
+    
+    #buffer limit
+    CHAINMAXLENGTH = 100000
+    chainX   = npy.zeros(CHAINMAXLENGTH,dtype = 'int32', order='C')
+    chainY   = npy.zeros(CHAINMAXLENGTH,dtype = 'int32', order='C')
+    N   = npy.zeros(1,dtype = 'int32', order='C')
+        
+    #activate or de-activate verbose in C code
+    if verbose:
+        VERBOSE = '#define VERBOSE True'
+    else:
+        VERBOSE = '#undef VERBOSE'
+        
+    #code contains the main watershed C code
+    code = \
+"""
+%s
+#define CHAINMAXLENGTH %d 
+int *IN      = (int *) PyArray_GETPTR1(ima_array,0);
+int *CHAINX      = (int *) PyArray_GETPTR1(chainX_array,0);
+int *CHAINY      = (int *) PyArray_GETPTR1(chainY_array,0);
+int *n      = (int *) PyArray_GETPTR1(N_array,0);
+
+#ifdef  VERBOSE    
+    printf("(chain.c)\\n");
+#endif
+
+int N_IN = PyArray_SIZE(ima_array);
+int m_IN = PyArray_DIM(ima_array,0); 
+int n_IN = PyArray_DIM(ima_array,1);
+int sizex = n_IN;
+int sizey = m_IN;  
+
+// call the function defined in the chain.c file
+
+//return an int to Python
+if(ext)
+    *n = build_chain4ex(IN,sizex,sizey,label,CHAINX,CHAINY);
+else
+    *n = build_chain(IN,sizex,sizey,label,CHAINX,CHAINY);
+
+"""% (VERBOSE,CHAINMAXLENGTH)
+
+#    print code
+
+    contours = []
+    
+    extra_code = '#define CHAINMAXLENGTH %d '%(CHAINMAXLENGTH) + open(os.path.join(os.path.split(__file__)[0],'chain.c')).read()
+    for label in labels:
+        inline(code, ['ima','chainX','chainY','label','ext','N'],support_code=extra_code)
+        n = N[0]
+        contours.append((label,n,npy.copy(chainX[0:n]),npy.copy(chainY[0:n])))
+    return contours
+    
+def testChain():
+    """open a binarised test image, compute internal and external chain and display it in overlay
+    """
+    im = scipy.misc.imread('../../test/data/blocs.png')
+    m = measurements.label(im)
+    label = m[0]
+    print label
+    scipy.misc.imsave('../../test/temp/chainLabel.tif',label)
+    
+    contours = chain(label,labels = None,ext = 0)
+    contoursExt = chain(label,labels = None,ext = 1)
+    
+    
+    pylab.figure(1)
+    pylab.imshow(label, interpolation='nearest')
+    
+    for c in contours:
+        print c
+        
+        pylab.plot(c[2]-1,c[3]-1)
+        x0 = c[2][0]
+        y0 = c[3][0]
+        xg = c[2].mean()
+        yg = c[3].mean()
+        pylab.text(x0,y0,'label%d'%c[0],color=[.5,.5,.5])
+        print len(c[2])
+        
+#    for c in contoursExt:
+#        pylab.plot(c[2]-1,c[3]-1)
+        
+    pylab.show()
+
+
+if __name__ == "__main__":
+    import doctest
+    doctest.testmod()
+
+    testChain()
+#    help(chain)
+    

dependencies/pymic/chain/segmentation.py

+__author__ = 'olivier'
+import matplotlib.pyplot as plt
+import matplotlib.cm as cm
+import scipy as scp
+import numpy as np
+import scipy.ndimage.measurements  as measure
+from pprint import pprint
+from numpy.linalg import norm
+import itertools as its
+from scipy import ndimage
+
+import networkx as nx
+
+from chain import chain
+
+def lobe_fct(theta0,theta):
+    """return lobe function
+    rho=1.0 in theta=0.0
+    rho=0.0 outside [-theta0,+theta0]
+    e.g.
+
+    >>> import matplotlib.pyplot as plt
+    >>> import matplotlib.cm as cm
+    >>> import numpy as np
+
+    >>> theta = np.linspace(-np.pi,+np.pi,1000)
+
+    >>>for s in np.arange(1.0,5.0):
+    >>>    theta0 = np.pi/s
+    >>>    rho = lobe_fct(theta0*1.5,theta)
+    >>>    plt.polar(theta,rho,label='s=%d'%s)
+    >>>    plt.polar([theta0,theta0],[0,1],'k:')
+    >>>    plt.polar([-theta0,-theta0],[0,1],'k:')
+    >>>plt.legend()
+    >>>plt.show()
+
+    """
+    s = .5*np.pi/(theta0)
+    rho = np.cos(s*theta)**4
+    rho[theta<-theta0]=0.0
+    rho[theta>theta0]=0.0
+    return rho
+
+def lobe_f(theta0,theta):
+    s = .5*np.pi/(theta0)
+    rho = np.cos(s*theta)**4
+    if theta<-theta0:
+        rho = 0.0
+    if theta>theta0:
+        rho=0.0
+    return rho
+
+def test_lobe():
+    theta = np.linspace(-np.pi, +np.pi, 1000)
+    for s in np.arange(1.0, 5.0):
+        theta0 = np.pi / s
+        print theta0
+        rho = lobe_fct(theta0, theta)
+        #        plt.plot(theta,rho,label='s=%d'%s)
+        plt.polar(theta, rho, label='theta_0 = pi/%d' % s)
+        plt.polar([theta0, theta0], [0, 1], 'k:')
+        plt.polar([-theta0, -theta0], [0, 1], 'k:')
+    plt.legend()
+    plt.show()
+
+def zero_cross(x):
+    idx0 = np.arange(len(x))
+    idx1 = np.roll(idx0,-1)
+    return (x[idx1]*x[idx0])<=0
+
+
+def cyclic_derivative(x):
+    idx0 = np.arange(len(x))
+    idx1 = np.roll(idx0,-1)
+    idx2 = np.roll(idx0,+1)
+    return x[idx2]-x[idx1]
+
+class Node(object):
+    def __init__(self,label,x,y,parent_label,is_object):
+        self.x = x
+        self.y = y
+        self.parent_label = parent_label
+        self.parent = None
+        self.children = []
+        self.label = label
+        self.is_object = is_object
+
+def find_nodes(node_tree):
+    """returns all nodes from a tree
+    """
+    if node_tree.children == []:
+        return [node_tree]
+    r = [node_tree]
+    for c in node_tree.children:
+        r.extend(find_nodes(c))
+    return r
+
+def hierarchic_labelling(im):
+    """compute object(>0 pixels) labelling, hole and inclusions are detected
+    return contour list"""
+    label_objects = measure.label(im)
+    label_holes = measure.label(255-im)
+    n,m = im.shape
+    objects = measure.find_objects(label_objects[0])
+    holes = measure.find_objects(label_holes[0])
+    idx_valid_objects = range(1,label_objects[1]+1)
+    #reject bg obj touching the image border
+    idx_valid_holes = []
+    for i,obj in enumerate(holes):
+        slice_y,slice_x = obj
+        if (slice_x.start==0) | (slice_x.stop == m)\
+            |(slice_y.start==0) | (slice_y.stop == n):
+            pass
+        else:
+            idx_valid_holes.append(i+1)
+
+    #build contours
+    extern_contours = chain(label_objects[0],idx_valid_objects)
+    intern_contours = chain(label_holes[0],idx_valid_holes)
+
+    #build tree
+    plt.figure(4)
+    plt.imshow(label_holes[0],interpolation='nearest')
+    plt.colorbar()
+
+    nodes = {}
+    #contour graph get the object contour structure (object containing holes)
+    c_graph = nx.Graph()
+    for cont in extern_contours:
+        label,n,x,y = cont
+        sx = x[0]-2
+        sy = y[0]-1
+        plt.plot(x-1,y-1)
+        plt.plot(sx,sy,'o')
+        plt.text(x[0],y[0],'obj%d'%label)
+        neib_label = label_holes[0][sy,sx]
+        nodes['obj%d'%label] = Node(label,x,y,'bg%d'%neib_label,True)
+        contour_attr = {'x':x,'y':y,'parent':'bg%d'%neib_label,'is_object':True}
+        c_graph.add_node(label,contour_attr)
+
+    for cont in intern_contours:
+        label,n,x,y = cont
+        sx = x[0]-2
+        sy = y[0]-1
+        plt.plot(x-1,y-1)
+        plt.plot(sx,sy,'o')
+        plt.text(x[0],y[0],'bg%d'%label)
+        neib_label = label_objects[0][sy,sx]
+        nodes['bg%d'%label] = Node(label,x,y,'obj%d'%neib_label,False)
+        contour_attr = {'x':x,'y':y,'parent':'bg%d'%neib_label,'is_object':False}
+        c_graph.add_node(-label,contour_attr)
+        c_graph.add_edge(-label,neib_label)
+
+    plt.figure(14)
+    nx.draw(c_graph)
+    #extract connected component from the graph (i.e. object with its internal holes)
+    connected_components = nx.connected_component_subgraphs(c_graph)
+    
+
+    #build tree
+    roots = []
+    for k in nodes:
+        #if has parent,
+        n = nodes[k]
+        if n.parent_label in nodes:
+            nodes[k].parent = nodes[n.parent_label]
+            nodes[n.parent_label].children.append(nodes[k])
+        else:
+            roots.append(nodes[k])
+
+    plt.figure(5)
+    plt.title('roots')
+    for r in roots:
+        for n in find_nodes(r):
+            plt.plot(n.x,n.y)
+
+    #build connected parts
+    connect_parts = []
+    for k in nodes:
+        n = nodes[k]
+        if n.is_object:
+            part = [n]
+            for c in n.children:
+                part.append(c)
+            connect_parts.append(part)
+
+    plt.figure(6)
+    plt.title('connected parts')
+    for k in connect_parts:
+        x = k[0].x
+        y = k[0].y
+        for c in k[1:]:
+            x = np.hstack((x,c.x))
+            y = np.hstack((y,c.y))
+        plt.plot(x,y)
+
+    for k in connect_parts:
+        g = find_cp_graph(k)
+        links = compute_connectivity_graph(g,im)
+
+        plt.figure(6)
+
+        for link in links:
+            plt.plot([link[0][0],link[1][0]],[link[0][1],link[1][1]],'k',linewidth=3)
+
+        for i,sp in enumerate(g):
+            plt.plot(sp.x,sp.y,'or')
+            plt.plot(sp.x+sp.v1[0]*10,sp.y+sp.v1[1]*10,'+k')
+            plt.plot(sp.x+sp.v2[0]*10,sp.y+sp.v2[1]*10,'+k')
+            (x,y)  = sp.lobe()
+            plt.plot(x,y,'k')
+            dx = 20*sp.dx
+            dy = 20*sp.dy
+            plt.plot([sp.x,sp.x+dx],[sp.y,sp.y+dy],'k')
+            plt.text(sp.x,sp.y,'%d'%i)
+    plt.show()
+
+    
+def find_corners(w, x, y):
+    idx0 = np.array(range(len(x)))
+    idx1 = np.roll(idx0, w)
+    idx2 = np.roll(idx0, -w)
+    v1x = x[idx1] - x[idx0]
+    v2x = x[idx2] - x[idx0]
+    v1y = y[idx1] - y[idx0]
+    v2y = y[idx2] - y[idx0]
+    v1 = np.hstack((v1x[:, np.newaxis], v1y[:, np.newaxis]))
+    v2 = np.hstack((v2x[:, np.newaxis], v2y[:, np.newaxis]))
+    #norm v1,v2
+    v1 = norm2d(v1)
+    v2 = norm2d(v2)
+    v3 = np.cross(v1, v2)
+
+    plt.figure(100)
+    plt.plot(v3)
+
+    derv3 = cyclic_derivative(v3)
+    zeros = zero_cross(derv3)
+    p_in = zeros & (v3 < -0.4)
+    p_out = zeros & (v3 > 0.6)
+    temp,p_in = find_boolean_clusters(p_in,5)
+    temp,p_out = find_boolean_clusters(p_out,5)
+    return p_in, p_out, v1,v2,v3
+
+def norm2d(x):
+    """returns normalized vectors
+    """
+    n = np.sqrt(x[:,0]**2+x[:,1]**2)
+    return x/n[:,np.newaxis]
+
+class SalientPoint(object):
+    def __init__(self,x,y,v1,v2,v3):
+        self.x = x
+        self.y = y
+        self.v1 = v1
+        self.v2 = v2
+        self.v3 = v3
+        self.angle = np.arccos(np.dot(self.v1,self.v2))
+        if np.sign(self.v3) < 0:
+            self.direction = np.arctan2(-self.v2[1],-self.v2[0])+self.angle/2.0
+
+        else:
+            self.direction = np.arctan2(-self.v1[1],-self.v1[0])+self.angle/2.0
+        self.dx = np.cos(self.direction)
+        self.dy = np.sin(self.direction)
+    def lobe(self):
+        theta = np.linspace(-np.pi, +np.pi, 30)
+        rho = lobe_fct(self.angle, theta)
+        x = self.x + 10*rho * np.cos(self.direction+theta)
+        y = self.y + 10*rho * np.sin(self.direction+theta)
+        return (x,y)
+
+    def distance(self,x,y):
+        """returns a connectivity measure to a point depending on direction and distance
+        """
+        vx = x-self.x
+        vy = y-self.y
+        dist = norm([vx,vy])
+        direction = np.arctan2(vy,vx)
+        return dist
+
+    def mutual(self,sp):
+        vx = sp.x-self.x
+        vy = sp.y-self.y
+        dist = norm([vx,vy])
+        direction = np.array([vx,vy])/dist
+        alpha = np.dot(-direction,np.array([sp.dx,sp.dy]))
+        beta = np.dot(direction,np.array([self.dx,self.dy]))
+
+        print (2-min(alpha,beta))*dist,dist,alpha,beta
+        linked = (dist<150) & (min(alpha,beta)>.8)
+        return linked
+
+def profile(ima,p0,p1,num=None):
+    """extract interpolated profile along a segment in an image
+    """
+    if num is None:
+        num = norm([p0[0]-p1[0],p0[1]-p1[1]])
+    n = np.linspace(p0[0],p1[0],num)
+    m = np.linspace(p0[1],p1[1],num)
+    return [n,m,ndimage.map_coordinates(ima, [m,n], order=0)]
+
+def compute_connectivity_graph(sp_list,im):
+    """returns a list of connected sp
+    connection that touch the bg are discarded"""
+    con_list = []
+    c = its.combinations(sp_list,2)
+    for s1,s2 in c:
+        m,n,v =profile(im,(s1.x,s1.y),(s2.x,s2.y))
+        if 0 not in v[3:-2]:
+            con_list.append(((s1.x,s1.y),(s2.x,s2.y)))
+    return con_list
+
+def find_best_cuts(graph):
+    """returns a simplified graph
+    """
+
+#        linked = s1.mutual(s2)
+#        if linked:
+#            #check with image (no bg pixels allowed
+#            m,n,v =profile(im,(s1.x,s1.y),(s2.x,s2.y))
+#            if 0 not in v[3:-2]:
+#                con_list.append(((s1.x,s1.y),(s2.x,s2.y)))
+
+    return graph
+
+def find_cp_graph(cp):
+    """find corners for a connected part
+    """
+    w = 8
+    #object
+    pin,pout,v1,v2,v3 = find_corners(w, cp[0].x, cp[0].y)
+    x_corners = cp[0].x[pin]
+    y_corners = cp[0].y[pin]
+    v1_corners = v1[pin,:]
+    v2_corners = v2[pin,:]
+    v3_corners = v3[pin]
+    g = []
+
+    #holes
+    for c in cp[1:]:
+        pin,pout,v1,v2,v3 = find_corners(w, c.x, c.y)
+        x_corners = np.hstack((x_corners,c.x[pout]))
+        y_corners = np.hstack((y_corners,c.y[pout]))
+        v1_corners = np.vstack((v1_corners,v1[pout,:]))
+        v2_corners = np.vstack((v2_corners,v2[pout,:]))
+        v3_corners = np.hstack((v3_corners,v3[pout]))
+        
+    for x,y,v1,v2,v3 in zip(x_corners,y_corners,v1_corners,v2_corners,v3_corners):
+        g.append(SalientPoint(x,y,v1,v2,v3))
+
+    return g
+
+def test_contours():
+    im = plt.imread('../../test/data/h_bubble.tif')
+    label_obj = measure.label(im)
+    label_holes = measure.label(255 - im)
+    idx = range(1, label_obj[1] + 1)
+    plt.figure(0)
+    plt.imshow(label_obj[0], origin='lower', interpolation='nearest')
+    plt.colorbar()
+    contours = chain(label_obj[0], idx)
+    print contours
+    w = 10
+    for cont in contours:
+        i, n, x, y = cont
+        plt.figure(1)
+        plt.plot(x - 1, y - 1)
+
+        p_in, p_out = find_corners(w, x, y)
+        plt.figure(1)
+        plt.plot(x[p_in], y[p_in], 'o:')
+        plt.plot(x[p_out], y[p_out], 'o:')
+    plt.show()
+
+def find_boolean_clusters(x,th=2):
+    """groups series of consecutive 1, 0 gap <= th are replaced by 1
+    cyclic, find first 0 to 1 transition
+    """
+    nz = np.nonzero(x)[0]
+    try:
+        offset = nz[nz>0][0]
+        x = np.roll(x,-offset)
+
+        fx = np.ones(x.shape).astype('bool')
+        fxc = np.zeros(x.shape).astype('bool')
+        pos = 0
+        status = True
+        icount = 0
+        while pos<len(x):
+            if status:
+                if x[pos]:
+                    icount += 1
+                else:
+                    ocount = 1
+                    status = False
+            else:
+                if x[pos]:
+                    if ocount>=th:
+                        fx[pos-ocount:pos] = False
+                        fxc[pos-icount/2-ocount-1] = True
+                        icount = 1
+                    else:
+                        icount += ocount+1
+                    status = True
+                else:
+                    ocount += 1
+            pos += 1
+        if ocount>=th:
+            fx[pos-ocount:pos] = False
+            fxc[pos-icount/2-ocount-1] = True
+        fx = np.roll(fx,offset)
+        fxc = np.roll(fxc,offset)
+        return (fx,fxc)
+    except :
+        return (x,x)
+
+if __name__ == "__main__":
+
+#    test_lobe()
+#    test_contours()
+
+#    im = plt.imread('../../test/data/bubble_chain.tif')
+    im = plt.imread('../../test/data/h_bubble.tif')
+    hierarchic_labelling(im)
+    x = np.array([1,0,1,0,1,1,0,0,1,0,0,0,1,0])
+    x = np.array([1,0,0,1,0,1,1,0,0,1,1,1,0,0,0,1,0,0,1,1])
+#    x = np.array([0,0])
+    print x
+#    print group(x,1).astype('int')
+    print find_boolean_clusters(x,2)[0].astype('int')
+    print find_boolean_clusters(x,2)[1].astype('int')
+#    print group(x,3).astype('int')
+#    print group(x,4).astype('int')
+
+

dependencies/pymic/color/__init__.py

+from rgbcombine import combineRGB,blue_cmap,green_cmap,red_cmap

dependencies/pymic/color/rgbcombine.py

+import matplotlib
+
+__author__ = 'olivier'
+
+import numpy as np
+import os
+
+import matplotlib.cm as cm
+import matplotlib.pyplot as plt
+
+#definition of the red/green/blue palet from black to single channel maximum
+_red = {'red': [(0.0,  0.0, 0.0),(1.0 ,  1.0, 1.0)],
+     'green': [(0.0,  0.0, 0.0),(1.0,  0.0, 0.0)],
+     'blue':  [(0.0,  0.0, 0.0),(1.0,  0.0, 0.0)]}
+_green = {'red':[(0.0,  0.0, 0.0),(1.0,  0.0, 0.0)],
+     'green': [(0.0,  0.0, 0.0),(1.0,  1.0, 1.0)],
+     'blue':  [(0.0,  0.0, 0.0),(1.0,  0.0, 0.0)]}
+_blue = {'red': [(0.0,  0.0, 0.0),(1.0,  0.0, 0.0)],
+     'green': [(0.0,  0.0, 0.0),(1.0,  0.0, 0.0)],
+     'blue':  [(0.0,  0.0, 0.0),(1.0,  1.0, 1.0)]}
+
+red_cmap = matplotlib.colors.LinearSegmentedColormap('redmap', _red)
+green_cmap = matplotlib.colors.LinearSegmentedColormap('greenmap', _green)
+blue_cmap = matplotlib.colors.LinearSegmentedColormap('bluemap', _blue)
+
+def combineRGB(red=None,green=None,blue=None):
+    '''returns one RGB image using palet and 3 optional channels'''
+    if red is not None:
+        sh = red.shape
+    if green is not None:
+        sh = green.shape
+    if blue is not None:
+        sh = blue.shape
+
+    rgb = np.zeros((sh[0], sh[1], 4))
+
+    if red is not None:
+        rgb = rgb + matplotlib.cm.ScalarMappable(cmap=red_cmap).to_rgba(red)
+    if green is not None:
+        rgb = rgb + matplotlib.cm.ScalarMappable(cmap=green_cmap).to_rgba(green)
+    if blue is not None:
+        rgb = rgb + matplotlib.cm.ScalarMappable(cmap=blue_cmap).to_rgba(blue)
+
+    return rgb
+
+if __name__ == "__main__":
+    red,green = np.meshgrid(range(100),range(100))
+    blue = red+green
+    print red
+    print green
+    print blue
+
+    rgb = combineRGB(red=red,green=green,blue=blue)
+    print rgb
+    plt.imshow(rgb)
+    plt.show()

dependencies/pymic/examples/interactive_nucleus_detection.py

+__author__ = 'olivier'
+
+import struct
+import numpy as np
+import os
+
+from scipy.ndimage import measurements
+from scipy.fftpack import fft, ifft, fftshift, ifftshift
+from scipy.ndimage.filters import gaussian_filter
+from scipy.ndimage import measurements
+
+import matplotlib.cm as cm
+import matplotlib.pyplot as plt
+
+from enthought.traits.api import Range
+from enthought.traits.ui.api import Item, View, VSplit, HSplit, VGroup
+
+import sys
+sys.path.append('../')
+from rankfilter import rankfilter
+from watershed import watershed
+from zvi import zviread
+from uitools.interactive import Process, Tester
+
+
+def egal(ima,dtype='uint16',n = 1024):
+    fima = ima.flatten()
+    print max(fima)+1
+    h,bins = np.histogram(fima,bins = range(0,max(fima)+2))
+    ch = np.cumsum(h)
+    lut = n*ch/ch[-1]
+    return (lut[ima]).astype(dtype)
+
+def autolevel(ima,dtype='uint16',n = 1024):
+    fima = ima.flatten()
+    m = min(fima)
+    M = max(fima)
+    delta = float(M-m)/n
+    res = ((ima-m)/delta).astype(dtype)
+    return res
+
+def remove_small_obj(bw, area=10000):
+    labels, numfeat = measurements.label(bw)
+    area_obj = measurements.sum(bw, labels, index=range(numfeat + 1))
+    lut = area_obj > area
+    result = lut[labels]
+    return result
+
+def remove_small_holes(bw, area = 50000):
+    labels, numfeat = measurements.label(bw == 0)
+    area_obj = measurements.sum(bw==0, labels, index=range(numfeat + 1))
+    lut = area_obj > 50000
+    result = lut[labels]==0
+    return result
+
+class BgDetectionProcess(Process):
+    radius = Range(1.0,20,10)
+    threshold = Range(0.0,300,50)
+    traits_view = View(VGroup(Item('radius'),
+                    Item('threshold')))
+    def __init__(self):
+        #add here some initialisation if needed
+        pass
+
+    def apply(self,input):
+        #here is the code that return an image from the input image
+        al = autolevel(input)
+        gr = rankfilter(al,'gradient',self.radius,bitDepth=10)
+        #keep region with few borders
+        thresh = gr<self.threshold
+        #remove small foreground objects
+        bg = remove_small_obj(thresh,10000)
+        #remove small holes in background objects
+        bg = remove_small_holes(bg,50000)
+        return bg
+
+class NucleusDetectionProcess(Process):
+    radius = Range(1.0,20,10)
+    threshold = Range(0.0,2048,1024)
+    traits_view = View(VGroup(Item('radius'),
+                    Item('threshold')))
+    def __init__(self,input):
+        #add here some initialisation if needed
+        self.data = autolevel(input).astype('uint16')
+
+    def apply(self):
+        #here is the code that return an image from the input image
+
+#        gr = rankfilter(al,'gradient',self.radius,bitDepth=10)
+        #keep region with few borders
+#        thresh = gr<self.threshold
+        #remove small foreground objects
+#        bg = remove_small_obj(thresh,10000)
+        #remove small holes in background objects
+#        bg = remove_small_holes(bg,50000)
+        bg = self.data>self.threshold
+        return bg
+
+
+if __name__ == "__main__":
+    dir = ['../test/data/'+f for f in os.listdir('../test/data/') if 'zvi' in f]
+    print dir
+    dir = [dir[0]]
+    color = [cm.gray,cm.Blues,cm.Greens,cm.Oranges]
+
+
+    for i,fname in enumerate(dir):
+        dic = zviread(fname,0).Image.Array
+        dapi = zviread(fname,1).Image.Array
+
+#        tester = Tester(process,dic.astype('float'))
+        process = NucleusDetectionProcess(dapi.astype('float'))
+        tester = Tester(process)
+        tester.configure_traits()
+        print tester.output
+#        bg = background_detection(dic)
+#
+#
+#        plt.figure(i)
+#        plt.subplot(1,2,1)
+#        plt.imshow(dic,origin='lower',cmap=color[0],interpolation='nearest')
+#        plt.colorbar()
+#        plt.subplot(1,2,2)
+#        plt.imshow(bg,origin='lower',cmap=color[0],interpolation='nearest')
+#        plt.colorbar()
+    plt.show()
+

dependencies/pymic/examples/zvi_read_and_display.py

+'''
+Created on Nov 11, 2010
+
+test zviread python module
+
+@author: olivier
+'''
+import struct
+import numpy as np
+import matplotlib.cm as cm
+import matplotlib.pyplot as plt
+from scipy.ndimage import measurements
+
+import os
+import sys
+sys.path.append('../')
+from zvi import zviread
+
+import os
+
+
+
+
+def display_all():
+    
+#    dir = ['../../test/data/'+f for f in os.listdir('../../test/data/') if 'zvi' in f]
+    path = '/media/data/zeiss_zvi/Laure_tpima'
+#    path = '../../test/data'
+    dir = [os.path.join(path,f) for f in os.listdir(path) if 'zvi' in f]
+    print dir
+    color = [cm.gray,cm.Blues,cm.Greens,cm.Oranges]
+    
+    for i,fname in enumerate(dir):
+        print fname
+        plt.figure(i)
+        for p in range(4):
+            try:
+                b = zviread(fname,p).Image.Array
+                plt.subplot(2,2,p+1)
+                plt.imshow(b,origin='lower',cmap=color[p],interpolation='nearest')
+                plt.colorbar()
+            except :
+                print 'Plane %d unavailable'%p
+        plt.title(fname)
+
+    plt.show()
+
+
+
+if __name__ == "__main__":
+    display_all()

dependencies/pymic/fourier/__init__.py

+
+

dependencies/pymic/fourier/filter_build.py

+import wx
+
+from traits.api import HasTraits,Any, Instance, Button
+from mpl_figure_editor import MPLFigureEditor
+from traitsui.api import View, Item, HSplit, VGroup
+
+import numpy as np
+from numpy import fft
+import matplotlib.pyplot as plt
+from matplotlib.figure import Figure
+
+class FourierInteractive(HasTraits):
+
+    source_fig = Instance(Figure, ())
+    dest_fig = Instance(Figure, ())
+
+    transform_but = Button('FFT')
+
+    view = View(HSplit(
+                    HSplit(
+                        Item('source_fig', editor=MPLFigureEditor(), show_label=False),
+                        Item('dest_fig', editor=MPLFigureEditor(), show_label=False),
+                    ),
+                    VGroup(
+                        Item('transform_but', show_label=False),
+                        Item('transform_but', show_label=False),
+                        Item('transform_but', show_label=False),
+                    ),
+                ),
+        width=400,
+        height=300,
+        resizable=True)
+
+    def __init__(self,ima):
+        super(FourierInteractive, self).__init__()
+        self.ima = ima
+        self.source_axes = self.source_fig.add_subplot(111)
+        self.source_axes.imshow(self.ima,interpolation='nearest')
+        self.dest_axes = self.dest_fig.add_subplot(111)
+
+    def _transform_but_fired(self):
+        F = fft.fftshift(np.abs(fft.fft2(self.ima)))
+        self.dest_axes.imshow(F,interpolation='nearest')
+        wx.CallAfter(self.dest_fig.canvas.draw)
+
+
+if __name__ == "__main__":
+    # Create a window to demo the editor
+
+    ima = plt.imread('../../test/data/cameraman.tif')
+
+
+    app = FourierInteractive(ima)
+    app.configure_traits()

dependencies/pymic/fourier/mpl_figure_editor.py

+import wx
+
+import matplotlib
+# We want matplotlib to use a wxPython backend
+matplotlib.use('WXAgg')
+from matplotlib.backends.backend_wxagg import FigureCanvasWxAgg as FigureCanvas
+from matplotlib.figure import Figure
+from matplotlib.backends.backend_wx import NavigationToolbar2Wx
+
+from traits.api import Any, Instance
+from traitsui.wx.editor import Editor
+from traitsui.wx.basic_editor_factory import BasicEditorFactory
+
+class _MPLFigureEditor(Editor):
+
+    scrollable  = True
+
+    def init(self, parent):
+        self.control = self._create_canvas(parent)
+        self.set_tooltip()
+
+    def update_editor(self):
+        pass
+
+    def _create_canvas(self, parent):
+        """ Create the MPL canvas. """
+        # The panel lets us add additional controls.
+        panel = wx.Panel(parent, -1, style=wx.CLIP_CHILDREN)
+        sizer = wx.BoxSizer(wx.VERTICAL)
+        panel.SetSizer(sizer)
+        # matplotlib commands to create a canvas
+        mpl_control = FigureCanvas(panel, -1, self.value)
+        sizer.Add(mpl_control, 1, wx.LEFT | wx.TOP | wx.GROW)
+        toolbar = NavigationToolbar2Wx(mpl_control)
+        sizer.Add(toolbar, 0, wx.EXPAND)
+        self.value.canvas.SetMinSize((10,10))
+        return panel
+
+class MPLFigureEditor(BasicEditorFactory):
+
+    klass = _MPLFigureEditor
+
+
+if __name__ == "__main__":
+    # Create a window to demo the editor
+    from traits.api import HasTraits
+    from traitsui.api import View, Item
+    from numpy import sin, cos, linspace, pi
+
+    class Test(HasTraits):
+
+        figure = Instance(Figure, ())
+
+        view = View(Item('figure', editor=MPLFigureEditor(),
+            show_label=False),
+            width=400,
+            height=300,
+            resizable=True)
+
+        def __init__(self):
+            super(Test, self).__init__()
+            axes = self.figure.add_subplot(111)
+            t = linspace(0, 2*pi, 200)
+            axes.plot(sin(t)*(1+0.5*cos(11*t)), cos(t)*(1+0.5*cos(11*t)))
+
+    Test().configure_traits()

dependencies/pymic/levelsets/__init__.py

+from chenveze import chenveze

dependencies/pymic/levelsets/checkstop.py

+import numpy as npy
+
+def checkstop(old,new,dt):
+    """indicate whether we should performance further iterations or stop
+    """
+
+    if new.ndim==2: #grayscale image
+        ind = npy.abs(new)<=.5
+        M = npy.sum(ind)
+        print 'M',M
+        Q = npy.sum(npy.abs(new[ind]-old[ind]))/M
+        if Q<dt*.18**2:
+            indicator = 1
+        else:
+            indicator = 0
+
+    else: #rgb image
+        new0 = new[:,:,0]
+        new1 = new[:,:,1]
+        old0 = old[:,:,0]
+        old1 = old[:,:,1]
+        ind0 = npy.abs(old0)<1
+        ind1 = npy.abs(old1)<1
+        M0 = npy.sum(ind0)
+        M1 = npy.sum(ind1)
+        Q0 = npy.sum(npy.abs(new0[ind0]-old0[ind0]))/M0
+        Q1 = npy.sum(npy.abs(new1[ind1]-old1[ind1]))/M1
+        if Q0<dt*.18**2 & Q1<dt*.18**2:
+            indicator = 1
+        else:
+            indicator = 0
+
+    return indicator
+
+def main():
+    print 'no test'
+
+if __name__ == '__main__':
+    main()
+

dependencies/pymic/levelsets/chenveze.py

+import numpy as npy
+import matplotlib.pyplot as plt
+import matplotlib.cm as cm
+
+from scipy.ndimage import distance_transform_cdt,gaussian_filter
+from scipy.ndimage.filters import convolve,convolve1d
+from maskcircle2 import maskcircle2
+from heaviside import heaviside
+from kappa import kappa
+from checkstop import checkstop
+
+
+def chenveze(I,mask,num_iter,mu=.2,method='chan'):
+    """Active contour with Chen-Vese Method
+        for image segementation
+
+        Implemented by Yue Wu (yue.wu@tufts.edu)
+        Tufts University
+        Feb 2009
+        http://sites.google.com/site/rexstribeofimageprocessing/
+
+        all rights reserved
+        Last update 02/26/2009
+        --------------------------------------------------------------------------
+        Usage of varibles:
+        input:
+           I           = any gray/double/RGB input image
+           mask        = initial mask, either customerlized or built-in
+           num_iter    = total number of iterations
+           mu          = chenvezeweight of length term
+           method      = submethods pick from ('chen','vector','multiphase')
+
+        Types of built-in mask functions
+           'small'     = create a small circular mask
+           'medium'    = create a medium circular mask
+           'large'     = create a large circular mask
+           'whole'     = create a mask with holes around
+           'whole+small' = create a two layer mask with one layer small
+                           circular mask and the other layer with holes around
+                           (only work for method 'multiphase')
+        Types of methods
+           'chan'      = general CV method
+           'vector'    = CV method for vector image
+           'multiphase'= CV method for multiphase (2 phases applied here)
+
+        output:
+           phiS0        = updated level set function
+    """
+
+    eps = 1.0e-10
+
+    #check mask
+    mask = maskcircle2(I,mask)
+
+    #initialisation
+    P = I.astype(float)
+    pmin = npy.min(P)
+    pmax = npy.max(P)
+    prange = pmax-pmin
+    P = (P-pmin)/prange
+    force = eps
+
+    #core function
+    if method == 'chan':
+        #SDF
+        #get the distance map of the initial mask
+#        phi0 = distance_transform_cdt(mask) - distance_transform_cdt(1-mask) + mask - .5
+        phi0 = distance_transform_cdt(mask) - distance_transform_cdt(1-mask) - mask
+
+        #initial force, set to eps to avoid division by zeros
+
+        #main loop
+        fig = plt.figure()
+        for n in range(num_iter):
+            inidx = phi0>=0
+            outidx = phi0<0
+
+            L = P
+            c1 = npy.sum(L*(heaviside(phi0)))/(npy.sum(inidx)+eps)
+            c2 = npy.sum(L*(1.0-heaviside(phi0)))/(npy.sum(outidx)+eps)
+            force_image = - (L-c1)**2 + (L-c2)**2
+            #sum on all the band
+            #todo
+
+            plt.clf()
+
+            plt.subplot(2,2,1)
+            plt.imshow(force_image)
+            plt.title('F')
+            plt.colorbar()
+            plt.subplot(2,2,2)
+            plt.imshow(phi0)
+            plt.title('phi')
+            plt.colorbar()
+            print n,c1,c2
+            # calculate the external force of the image
+            force = mu*kappa(phi0)/npy.max(npy.abs(kappa(phi0)))+force_image
+            plt.subplot(2,2,3)
+            plt.imshow(phi0<0)
+            plt.title('seg')
+            plt.colorbar()
+            #display graylevel distribution for in and out
+            hin = npy.histogram(P[inidx].flatten(),npy.arange(-1,1,.01))
+            hout = npy.histogram(P[outidx].flatten(),npy.arange(-1,1,.01))
+            plt.subplot(2,2,4)
+            plt.plot(hin[0],label='in')
+            plt.plot(hout[0],label='out')
+            plt.title('distribution')
+            plt.legend()
+            fig.savefig('../../test/temp/iter%04d.png'%n)
+
+            #normalize the force
+            force = 1.0 * force / npy.max(npy.abs(force))
+
+            #get step size dt
+            dt = 10.0
+
+            #stopping criterion
+            old = phi0
+            phi0 = phi0 + dt*force
+            new = phi0
+
+            indicator = checkstop(old,new,dt)
+
+#            if false:
+#                seg = phi0<0
+#                return seg
+
+    seg = phi0<0
+    mu_in = npy.sum(I*(heaviside(phi0)))/(npy.sum(inidx)+eps)
+    mu_out = npy.sum(I*(1.0-heaviside(phi0)))/(npy.sum(outidx)+eps)
+    return [seg,mu_in,mu_out]
+
+
+def Cx(ima):
+    """x' derivative of image"""
+    c = convolve1d(ima,npy.array([-1,0,1]),axis=1,cval=1)
+    return c/2.0
+
+def Cy(ima):
+    """y' derivative of image"""
+    c = convolve1d(ima,npy.array([-1,0,1]),axis=0,cval=1)
+    return c/2.0
+
+def border(I):
+    Ix = Cx(I)
+    Iy = Cy(I)
+    return npy.maximum(npy.abs(Ix),npy.abs(Iy))
+
+def local_orientation(I,s=1.5):
+#    I = gaussian_filter(I,s)
+    Ix = Cx(I)
+    Iy = Cy(I)
+    angle = npy.arctan2(Iy,Ix)
+    angle = gaussian_filter(angle,s)
+    return angle
+
+def demo_texture():
+#    I = plt.imread('../../test/data/texture.png')[:,:,0]
+    I = plt.imread('../../test/data/hlines_mix2.tif')
+    #add some noise
+    n = npy.random.random(I.shape)
+    n_th = .1
+    I[n>n_th] = 255-I[n>n_th]
+
+    A = local_orientation(I,3)
+
+    plt.figure()
+    plt.subplot(1,2,1)
+    plt.imshow(I)
+    plt.subplot(1,2,2)
+    plt.imshow(A)
+    plt.colorbar()
+
+    S,mu_in,mu_out = chenveze(A,mask='large',num_iter=20,mu=.2,method='chan')
+
+    print mu_in,mu_out
+
+    plt.figure()
+    plt.imshow(I,cmap=cm.gray)
+    plt.imshow(S*255,cmap=cm.jet,alpha=.8)
+    plt.xlabel('$\mu_{in}=%f ^{\circ}$  $\mu_{out}=%f ^{\circ}$'%(180*mu_in/npy.pi,180*mu_out/npy.pi))
+    plt.show()
+
+
+def main():
+
+#    I = plt.imread('../../test/data/testpng.png')[:,:,0]
+
+#    I = plt.imread('../../test/data/image4.png')[:,:,0]
+#    I = local_orientation(I,1)
+#    I = plt.imread('../../test/data/exp0012.jpg')
+#    I = plt.imread('../../test/data/dh_phase.png')
+#    I = plt.imread('../../test/data/brain.jpg')[:,:,1]
+#    I = plt.imread('../../test/data/test.png')*255.0
+#    I = plt.imread('../../test/data/flowers.jpg')[:,:,0]
+#    I = plt.imread('../../test/data/cameraman.tif')
+#    I = plt.imread('../../test/data/phase_s.png')
+    print I.shape
+
+    plt.imshow(I)
+    plt.colorbar()
+    plt.show()
+
+    S = chenveze(I,mask='whole',num_iter=100,mu=.2,method='chan')
+    plt.show()
+
+    plt.imshow(S)
+    plt.colorbar()
+    plt.show()
+if __name__ == '__main__':
+    demo_texture()
+

dependencies/pymic/levelsets/heaviside.py

+import numpy as npy
+
+def heaviside(z):
+    """Heaviside step function (smoothed version)"""
+    epsilon = 1e-1
+    H = npy.zeros(z.shape)
+    idx1 = z>epsilon
+    idx2 = npy.logical_and(z<epsilon,z>-epsilon)
+    ze = z[idx2]/epsilon
+
+    H[idx1] = 1
+    H[idx2] = .5 *(1 + ze + 1.0/npy.pi * npy.sin(npy.pi*ze))
+
+
+    return H
+
+def main():
+    import matplotlib.pyplot as plt
+
+    z = npy.linspace(-1,1,1000)
+    h = heaviside(z)
+    plt.plot(h)
+    plt.show()
+
+if __name__ == '__main__':
+    main()
+

dependencies/pymic/levelsets/kappa.py

+import numpy as npy
+
+def kappa(I):
+    """get curvature information of input image
+    input: 2D image I
+    output: curvature matrix KG (size+=1)
+    """
+    m,n = I.shape
+    P = npy.hstack((npy.zeros((m,1)),I,npy.zeros((m,1))))
+    m,n = P.shape
+    P = npy.vstack((npy.zeros((1,n)),P,npy.zeros((1,n))))
+
+    fy = P[2:,1:-1] - P[0:-2,1:-1]
+    fx = P[1:-1,2:] - P[1:-1,0:-2]
+    fyy = P[2:,1:-1] + P[0:-2,1:-1] - 2 * I
+    fxx = P[1:-1,2:] + P[1:-1,0:-2] - 2 * I
+    fxy = .25*(P[2:,2:]+P[:-2,:-2]-P[:-2,2:]-P[2:,:-2])
+    G = npy.sqrt(fx**2+fy**2)
+    K = (fxx*fy**2-2*fxy*fx*fy+fyy*fx**2)/(npy.power(fx**2+fy**2,1.5)+1e-10)
+    KG = K*G
+
+    return KG
+
+def main():
+    import matplotlib.pyplot as plt
+
+    I = plt.imread('brain.jpg')[:,:,1]
+
+    KG = kappa(I)
+
+    plt.imshow(KG)
+    plt.colorbar()
+    plt.show()
+
+if __name__ == '__main__':
+    main()
+

dependencies/pymic/levelsets/maskcircle2.py

+import numpy as npy
+from scipy.ndimage.filters import convolve
+
+def maskcircle2(I,type='small'):
+    """auto pick a circular mask for image I
+    built-in mask creation function
+    Input: I   : input image
+    type: mask shape keywords
+     Output: m  : mask image
+     """
+    if I.ndim!=2:
+        temp = I[:,:,0]
+    else:
+        temp = I
+
+    h = npy.asarray([[0,1,0],[1,-4,1],[0,1,0]])
+    T = convolve(temp,h)
+    T[0,:] = 0
+    T[-1,:] = 0
+    T[:,0] = 0
+    T[:,-1] = 0
+
+    thre = npy.max(T)*.5
+    idx = npy.abs(T)>thre
+    n_idx = npy.sum(idx)
+
+    x,y = npy.meshgrid(range(0,temp.shape[1]),range(0,temp.shape[0]))
+
+    print x.shape,y.shape,T.shape
+
+
+    cx = npy.sum(x[idx])/n_idx
+    cy = npy.sum(y[idx])/n_idx
+
+    m = npy.zeros(temp.shape)
+    p,q = temp.shape
+
+    if type=='small':
+        r = 10
+        n = npy.zeros(x.shape)
+        n[((x-cx)**2+(y-cy)**2)<r**2] = 1
+        m = n
+    elif type=='medium':
+        r = min(min(cx,p-cx),min(cy,q-cy))
+        r = max(2/3*r,25)
+        n = npy.zeros(x.shape)
+        n[((x-cx)**2+(y-cy)**2)<r**2] = 1
+        m = n
+    elif type=='large':
+        r = min(min(cx,p-cx),min(cy,q-cy))
+        r = max(2/3*r,60)
+        n = npy.zeros(x.shape)
+        n[((x-cx)**2+(y-cy)**2)<r**2] = 1
+        m = n
+    elif type=='whole':
+        n = npy.zeros(x.shape)
+        n[10:-10,10:-10] = 1
+        m = n
+
+
+
+
+
+    return m
+
+
+
+def main():
+    import matplotlib.pyplot as plt
+
+    I = plt.imread('brain.jpg')[:,:,1]
+
+    m = maskcircle2(I,type='small')
+    plt.imshow(m)
+    plt.figure()
+    m = maskcircle2(I,type='medium')
+    plt.imshow(m)
+    plt.figure()
+    m = maskcircle2(I,type='large')
+    plt.imshow(m)
+    m = maskcircle2(I,type='whole')
+    plt.imshow(m)
+    plt.show()
+
+
+if __name__ == '__main__':
+    main()
+

dependencies/pymic/levelsets/reinitialization.py

+import numpy as npy
+
+def reinitialization(D,dt):
+    """reinitialize the distance map for active contour
+    """
+    m,n = D.shape
+    T = npy.hstack((npy.zeros((m,1)),D,npy.zeros((m,1))))
+    m,n = T.shape
+    T = npy.vstack((npy.zeros((1,n)),T,npy.zeros((1,n))))
+
+    a = D - T[:-2,1:-1]
+    b = T[2:,1:-1] - D
+    c = D - T[2:,:-2]
+    d = T[1:-1,2:]-D
+
+    a_p = npy.maximum(a,0)
+    a_m = npy.minimum(a,0)
+    b_p = npy.maximum(b,0)
+    b_m = npy.minimum(b,0)
+    c_p = npy.maximum(c,0)
+    c_m = npy.minimum(c,0)
+    d_p = npy.maximum(d,0)
+    d_m = npy.minimum(d,0)
+
+    G = npy.zeros(D.shape)
+    ind_plus = D>0
+    ind_minus = D<0
+
+    G[ind_plus] = npy.sqrt(npy.maximum(a_p[ind_plus]**2,b_m[ind_plus]**2)
+    + npy.maximum(c_p[ind_plus]**2,d_m[ind_plus]**2))-1
+    G[ind_minus] = npy.sqrt(npy.maximum(a_m[ind_minus]**2,b_p[ind_minus]**2)
+    + npy.maximum(c_m[ind_minus]**2,d_p[ind_minus]**2))-1
+
+    sign_D = D/npy.sqrt(D**2+1)
+    D = D - dt * sign_D * G
+
+    return D
+
+def main():
+    import matplotlib.pyplot as plt
+
+    X,Y = npy.meshgrid(npy.linspace(-1,1,200),npy.linspace(-1,1,200))
+    Z = X**2+Y**2
+
+    D = reinitialization(Z,1)
+
+    plt.imshow(D)
+    plt.show()
+
+if __name__ == '__main__':
+    main()
+

dependencies/pymic/meanshift/__init__.py

+from meanshift import meanshift, LUT

dependencies/pymic/meanshift/integralshift.py

+# -*- coding: utf-8 -*-
+'''
+Mean shift implementation using integral images
+integral function is extracted from opencv library
+'''
+
+__author__ = 'olivier'
+
+import sys
+import cv
+import numpy as np
+
+import pylab
+
+import matplotlib.patches as mpatches
+import matplotlib.lines as mlines
+
+from scipy.ndimage import imread
+
+def test_image(s=0):
+    t = np.zeros((100,120),dtype='uint8')
+    if s==0:
+        t[30:50,20:40] = 255
+    if s==1:
+        x,y = np.meshgrid(range(20,40),range(30,50))
+        t[30:50,20:40] = x
+    if s==2:
+        x,y = np.meshgrid(range(20,40),range(30,50))
+        t[30:50,20:40] = 55 - np.sqrt((x-30)**2+(y-40)**2)
+    if s==3:
+        t = imread('../../test/data/exp0001crop.jpg').astype(float)
+        #normalization
+        t = (t-t.min())/(t.max()-t.min())
+        t = t**2
+    return t
+
+def get_sum(integral_image,x1,y1,x2,y2):
+    """returns sum from CVarr integral image
+    """
+    a = cv.Get2D(integral_image,y1,x1)[0]
+    b = cv.Get2D(integral_image,y1,x2)[0]
+    c = cv.Get2D(integral_image,y2,x2)[0]
+    d = cv.Get2D(integral_image,y2,x1)[0]
+    return c+a-b-d
+
+class IntegratedImage(object):
+    """Object keeping the integral image and the X,Y integral images (internal opencv)
+    """
+    def __init__(self,array):
+        self.m,self.n = array.shape
+        X,Y = np.meshgrid(range(1,self.n+1),range(1,self.m+1))
+
+        self.mat = cv.fromarray(array)
+        x_mat = cv.fromarray((array*X).astype(float))
+        y_mat = cv.fromarray((array*Y).astype(float))
+
+        self.sum = cv.CreateMat(self.m+1,self.n+1,cv.CV_64F) #!!!!! +1
+        self.xsum = cv.CreateMat(self.m+1,self.n+1,cv.CV_64F) #!!!!! +1
+        self.ysum = cv.CreateMat(self.m+1,self.n+1,cv.CV_64F) #!!!!! +1
+
+        cv.Integral(self.mat,self.sum)
+        cv.Integral(x_mat,self.xsum)
+        cv.Integral(y_mat,self.ysum)
+
+    def _valid_rectangle(self,rect):
+        """return a valid rectangle w.r.t. image size
+        """
+
+        x1,x2 = (min(rect[0],rect[2]),max(rect[0],rect[2]))
+        y1,y2 = (min(rect[1],rect[3]),max(rect[1],rect[3]))
+        x1 = int(max(0,x1))
+        x2 = int(min(self.n,max(x2,0)))
+        y1 = int(max(0,y1))
+        y2 = int(min(self.m,max(y2,0)))
+
+        return (x1,y1,x2,y2)
+    def _get_sum(self,integral_image,rect):
+        """returns sum from CVarr integral image
+        """
+        x1,y1,x2,y2 = self._valid_rectangle(rect)
+
+        a = cv.Get2D(integral_image,y1,x1)[0]
+        b = cv.Get2D(integral_image,y1,x2)[0]
+        c = cv.Get2D(integral_image,y2,x2)[0]
+        d = cv.Get2D(integral_image,y2,x1)[0]
+        
+        return c+a-b-d
+
+    def _surf_rect(self,rect):
+        return (rect[2]-rect[0])*(rect[3]-rect[1])
+
+    def get_value(self,x,y):
+        return cv.Get2D(self.mat,int(y),int(x))[0]
+    
+    def find_g(self,rect):
+        """returns centroid of the rectangle
+        """
+        sum = self._get_sum(self.sum,rect)
+        xsum = self._get_sum(self.xsum,rect)
+        ysum = self._get_sum(self.ysum,rect)
+
+        if sum>0:
+            xg = xsum/sum
+            yg = ysum/sum
+            mean = sum/self._surf_rect(rect)
+        else:
+            xg = (rect[0]+rect[2])/2.0
+            yg = (rect[1]+rect[3])/2.0
+            mean = 0
+
+        return (xg,yg,mean)
+
+
+def shift(target,x0,y0,w,adapt = 'none'):
+    """converge using meanshift
+    """
+    N = 40
+    path = np.zeros((N,2))
+    x1 = x0-w/2.0
+    y1 = y0-w/2.0
+    x2 = x0+w/2.0
+    y2 = y0+w/2.0
+    path[0,:] = [x0,y0]
+    w0 = w
+    xg = x0
+    yg = y0
+    for iter in range(1,N):
+        rect = (x1,y1,x2,y2)
+        xg,yg,mean = target.find_g(rect)
+        mean_c= target.get_value(xg,yg)
+        path[iter,:] = [xg,yg]
+        d2 = (x1+w/2.0-xg)**2+(y1+w/2.0-yg)**2
+        if d2<.1:
+            path[iter+1,:] = [xg,yg]
+            path = path[0:iter+1,:]
+            break
+        if adapt is 'none':
+            pass
+        if adapt is 'iter_dec':
+            w = w*.8
+        if adapt is 'step_dec':
+            w = .5*np.sqrt(d2)
+        if adapt is 'intensity':
+            w = w0*(1-mean_c)
+        if adapt is 'custom':
+            weight = [1.0]*5 + [.75]*5 + [.5]*5 + [.25]*(N-15)
+            w = w0*weight[iter]
+        x1 = xg-w/2.0
+        y1 = yg-w/2.0
+        x2 = xg+w/2.0
+        y2 = yg+w/2.0
+
+    return (xg,yg,path)
+
+def test1():
+    """randow test
+    """
+    im = test_image(3)
+    pylab.imshow(im, interpolation='nearest')
+
+    target = IntegratedImage(im)
+    m,n = im.shape
+    print m,n
+    for i in range(10):
+
+        x1 = int(np.random.random()*n)-25
+        y1 = int(np.random.random()*m)-25
+        x2 = x1+25
+        y2 = y1+25
+
+        rect = (x1,y1,x2,y2)
+        xg,yg,mean = target.find_g(rect)
+
+        rect = mpatches.Rectangle((x1,y1),x2-x1,y2-y1,fill=False)
+        start = mpatches.Circle((xg,yg), 0.5,ec=[1,1,1],fc='none')
+        stop = mpatches.Circle(((x2+x1)/2.0,(y1+y2)/2.0), 0.5,fc=[1,1,1],ec='none')
+
+        arrow = mpatches.Arrow((x2+x1)/2.0,(y1+y2)/2.0, xg-(x2+x1)/2.0, yg-(y2+y1)/2.0 , width=0.5, color = 'y')
+
+        pylab.gca().add_patch(rect)
+        pylab.gca().add_patch(start)
+        pylab.gca().add_patch(stop)
+        pylab.gca().add_patch(arrow)
+
+    pylab.colorbar()
+
+def test2():
+    """grid test
+    """
+    im = test_image(3)
+
+    target = IntegratedImage(im)
+    m,n = im.shape
+    print m,n
+    w = 30
+    step = 5
+    pylab.imshow(im, interpolation='nearest')
+
+    for x in range(0,n,step):
+        for y in range(0,m,step):
+
+            x1 = x-w/2
+            y1 = y-w/2
+            x2 = x+w/2
+            y2 = y+w/2
+
+            rect = (x1,y1,x2,y2)
+            xg,yg,mean = target.find_g(rect)
+
+
+            rect = mpatches.Rectangle((x1,y1),x2-x1,y2-y1,fill=False)
+            start = mpatches.Circle((xg,yg), 0.5,fc=[1,1,1],ec='none')
+            stop = mpatches.Circle(((x2+x1)/2.0,(y1+y2)/2.0), 0.5,ec=[1,1,1],fc='none')
+
+            arrow = mpatches.Arrow((x2+x1)/2.0,(y1+y2)/2.0, xg-(x2+x1)/2.0, yg-(y2+y1)/2.0 , width=0.5, color = 'y')
+
+#            pylab.gca().add_patch(rect)
+            pylab.gca().add_patch(start)
+            pylab.gca().add_patch(stop)
+            pylab.gca().add_patch(arrow)
+
+    pylab.colorbar()
+
+def test3():
+    """iterations
+    """
+    im = test_image(3)
+
+    target = IntegratedImage(im)
+    m,n = im.shape
+    print m,n
+    w = 30
+    step = 10
+
+    fig = pylab.figure()
+    ax = pylab.axes([0,0,1,1])
+
+    pylab.imshow(im, interpolation='nearest')
+
+    for x in range(0,n,step):
+        for y in range(0,m,step):
+
+            xg,yg,path = shift(target,x,y,w,adapt='custom')
+
+            start = mpatches.Circle((x,y), 0.5,fc=[1,1,1],ec='none')
+            stop = mpatches.Circle((xg,yg), 0.5,ec=[1,1,1],fc='none')
+
+            arrow = mpatches.Arrow(x,y, xg-x, yg-y , width=0.5, color = 'y')
+            line = mlines.Line2D(path[:,0], path[:,1], lw=path.shape[0]*.3, color = 'y')
+
+#            pylab.gca().add_patch(start)
+            pylab.gca().add_patch(stop)
+            ax.add_line(line)
+
+
+    pylab.colorbar()
+
+if __name__ == "__main__":
+    import doctest
+    print cv
+    test3()
+    pylab.show()
+
+
+

dependencies/pymic/meanshift/meanshift.c

+#define MIN(a, b)  (((a) < (b)) ? (a) : (b))
+#define MAX(a, b)  (((a) > (b)) ? (a) : (b))
+
+#include <stdio.h>
+#include <math.h>
+
+int compute_g(unsigned char * pdata,int sizeX,int sizeY,double off_x,double off_y,double *TRIANGLE,double *DATA,double *LUT) {
+	int bbx0,bby0,bbx1,bby1;
+	int xint,yint;
+	long int pos;
+	unsigned char *pimage;
+
+	double a0,b0,c0,a1,b1,c1,a2,b2,c2;
+	double s0,s1,s2;
+	double x0,y0,x1,y1,x2,y2,x,y;
+
+	double mx0,my0,mx1,my1,mx2,my2;
+
+	double xcentre,ycentre,xr,yr;
+	double alpha,surf,surfalpha,value,totalvalue,gmin,gmax;
+	double n;
+
+	double sumX,sumY,sumXw,sumYw;
+
+	pimage = (unsigned char*)pdata;
+
+	mx0 = TRIANGLE[0]+off_x;
+	my0 = TRIANGLE[1]+off_y;
+	mx1 = TRIANGLE[2]+off_x;
+	my1 = TRIANGLE[3]+off_y;
+	mx2 = TRIANGLE[4]+off_x;
+	my2 = TRIANGLE[5]+off_y;
+
+	x0 = mx1;
+	y0 = my1;
+	x1 = mx0;
+	y1 = my0;
+	x2 = mx2;
+	y2 = my2;
+
+	//bounding box du triangle
+	bbx0 = (int)(MIN(MIN(x0,x1),x2))-1;
+	bbx1 = (int)(MAX(MAX(x0,x1),x2))+1;
+	bby0 = (int)(MIN(MIN(y0,y1),y2))-1;
+	bby1 = (int)(MAX(MAX(y0,y1),y2))+1;
+
+	//centre du triangle
+	xcentre = (x0+x1+x2)/3.0;
+	ycentre = (y0+y1+y2)/3.0;
+
+	//calcul des droites
+	if(x0==x1){
+		a0 = 1.0;
+		b0 = 0.0;
+		c0 = -x0;
+	}
+	else{
+		if(y0==y1){
+			a0 = 0.0;
+			b0 = 1.0;