Development snapshot.

Author: Stefan van der Walt

.bzrignore

@@ -1,3 +1,5 @@
 supreme/ext/.sconsign.dblite
 supreme/ext/*.os
-build
+build
+*.os
+.sconsign.dblite

supreme/doc/examples/match.py

@@ -0,0 +1,105 @@
+import os.path
+import glob
+
+from numpy.testing import set_local_path, restore_path
+import numpy as N
+import pylab as P
+
+set_local_path('../../..')
+import supreme as sr
+from supreme.config import data_path
+restore_path()
+
+print "Reading images and features..."
+features = []
+images = []
+
+dataset = 'reflectometer'
+basename = 'D'
+imagetype = 'png'
+featuretype = 'sift'
+T = 0.6
+image_files = sorted(glob.glob(os.path.join(data_path,'%s/%s*.%s' % (dataset,basename,imagetype))))
+feature_files = sorted(glob.glob(os.path.join(data_path,'%s/%s*.%s' % (dataset,basename,featuretype))))
+
+images = [sr.imread(fn,flatten=True) for fn in image_files]
+features = [sr.feature.SIFT.fromfile(fn,mode=featuretype.upper()) for fn in feature_files]
+
+for pair in zip(image_files,feature_files):
+    print pair[0], '->', pair[1]
+
+image_plane = N.hstack((images[0],images[1]))
+
+print "Matching features..."
+ref = features[0]
+tf_matrices = [N.eye(3)]
+valid_matrices = [True]
+for frame in features[1:]:
+    match,dist,valid = sr.feature.match(frame['data'],ref['data'],threshold=T)
+
+    valid_ref = match[valid]
+    xf,yf = frame['column'][valid],frame['row'][valid]
+    xr,yr = ref['column'][valid_ref],ref['row'][valid_ref]
+
+    M,ier = sr.register.sparse(yr,xr,yf,xf)
+    valid_matrices.append(ier == 1)
+    tf_matrices.append(M)
+
+    print "Found %d matches." % valid.sum()
+
+# Adjust coordinates for image plane
+rows,cols = images[0].shape
+xf += cols
+yf = rows - yf
+yr = rows - yr
+
+# Display matches
+#P.imshow(image_plane,cmap=P.cm.gray)
+#P.gca().set_autoscale_on(False)
+
+#P.plot(xf,yf,'.b')
+#P.plot(xr,yr,'.r')
+#xzip = zip(xf,xr)
+#yzip = zip(yf,yr)
+#for i in range(len(xzip)):
+#    P.plot(xzip[i:i+1],yzip[i:i+1],'-w')
+#P.show()
+
+images = [i for i,v in zip(images,valid_matrices) if v]
+tf_matrices = [t for t,v in zip(tf_matrices,valid_matrices) if v]
+
+# Scale for super-resolution
+scale = 3
+for M in tf_matrices:
+    M[:2,:] *= scale
+oshape = N.ceil(N.array(images[0].shape)*scale).astype(int)
+
+for img in images:
+    img -= img.min()
+print "Reconstructing using interpolation..."
+out1 = sr.register.stack.with_transform(images,tf_matrices,oshape=oshape)
+
+# Astronomy
+#out[out > 10] = 10
+#out /= out.max()
+
+
+print "Stacking using polygon overlap..."
+out2 = N.zeros(oshape,float)
+for i,(img,M) in enumerate(zip(images,tf_matrices)):
+    print "Stacking frame %d" % i
+    out2 += sr.ext.interp_transf_polygon(img,N.linalg.inv(M),oshape)
+out2 /= len(images)
+#out2[out2 > 500] = 500
+
+import scipy as S
+imsave = S.misc.pilutil.imsave
+imsave('original.png',images[0])
+imsave('_linear.png',out1)
+imsave('_polygon.png',out2)
+
+P.subplot(121)
+P.imshow(out1,interpolation='nearest',cmap=P.cm.gray)
+P.subplot(122)
+P.imshow(out2,interpolation='nearest',cmap=P.cm.gray)
+P.show()

supreme/ext/SConstruct

@@ -1,5 +1,5 @@
 from glob import glob
 
 env = Environment()
-env.Replace(CCFLAGS=['-O2','-Wall','-ansi','-pedantic'])
+env.Replace(CCFLAGS=['-O2','-ggdb','-Wall','-ansi','-pedantic'])
 env.SharedLibrary('supreme_', glob('*.c'))

supreme/ext/libsupreme.py

@@ -34,7 +34,7 @@ class POI(Structure):
    'variance_map' : (None,
                      [c_int, c_int, array_2d_uchar, array_2d_double,
                       c_int, c_int]),
-                     
+
    'line_intersect' : (None,
                        [c_double, c_double, c_double, c_double,
                         c_double, c_double, c_double, c_double,
@@ -52,11 +52,13 @@ class POI(Structure):
                        [c_int, c_int, array_2d_uchar,
                         c_int, c_int, array_2d_double, array_2d_double,
                         c_char, c_uint8, array_2d_uchar]),
+   'interp_transf_polygon': (None,
+                             [c_int, c_int, array_2d_uchar,
+                              c_int, c_int, array_2d_double,
+                              array_2d_double]),
     }
 
 def register_api(lib,api):
-    import inspect
-    parent_frame = inspect.currentframe().f_back
     for f, (restype, argtypes) in api.iteritems():
         func = getattr(lib, f)
         func.restype = restype
@@ -66,13 +68,13 @@ def register_api(lib,api):
 
 # Python wrappers for libsupreme functions
 
-def _atype(arrays, types):
+def atype(arrays, types):
     """Return contiguous arrays of given types.
 
     arrays - list of input arrays
     types - list of corresponding types
-    
-    """        
+
+    """
     out = ()
     try:
         za = zip(arrays,types)
@@ -81,28 +83,25 @@ def _atype(arrays, types):
 
     out = ()
     for A,T in za:
-        try:
-            out += N.ascontiguousarray(A,T),
-        except:
-            out += N.ascontiguousarray(A).astype(T),
+        out += N.ascontiguousarray(A,T),
 
     return out
-        
+
 def npnpoly(x_vertices, y_vertices, x_points, y_points):
     """Calculate whether points are in a given polygon.
 
     Returns a boolean array of length len(x_points).
-    
+
     """
-    xi,yi,x,y = _atype([x_vertices,y_vertices,
+    xi,yi,x,y = atype([x_vertices,y_vertices,
                         x_points,y_points],[N.double]*4)
 
     if xi[0] != xi[-1] or yi[0] != yi[-1]:
         xi = N.append(xi,x[0])
         yi = N.append(yi,y[0])
 
     out = N.empty(len(x),dtype=N.uint8)
-    
+
     _lib.npnpoly(len(xi), xi, yi,
                  len(x), x, y,
                  out)
@@ -118,7 +117,7 @@ def variance_map(image, shape=(3,3)):
     the variance of all the pixels in the window is stored.
 
     """
-    image, = _atype(image,N.uint8)
+    image, = atype(image,N.uint8)
     assert image.ndim == 2, "Image must be 2-dimensional"
     window_size_rows, window_size_columns = shape
     rows, columns = image.shape
@@ -151,16 +150,16 @@ def poly_clip(x, y, xleft, xright, ytop, ybottom):
     x and y are 1D arrays describing the coordinates of the vertices.
     xleft, xright, ytop and ybottom specify the borders of the
     bounding box.  Note that a cartesian axis system is used such that
-    the following must hold:
+    the following must hold true:
 
     x_left < x_right
     y_bottom < y_top
 
     The x and y coordinates of the vertices of the resulting polygon
     are returned.
-    
+
     """
-    x,y = _atype([x,y],[N.double,N.double])
+    x,y = atype([x,y],[N.double,N.double])
 
     assert len(x) == len(y), "Equal number of x and y coordinates required"
     assert ytop > ybottom
@@ -170,11 +169,11 @@ def poly_clip(x, y, xleft, xright, ytop, ybottom):
     if x[0] != x[-1] or y[0] != y[-1]:
         x = N.append(x,x[0])
         y = N.append(y,y[0])
-        
+
     xleft,xright,ytop,ybottom = map(N.double,[xleft,xright,ytop,ybottom])
 
-    workx = N.empty(2*len(x-1),dtype=N.double)
-    worky = N.empty(2*len(x-1),dtype=N.double)
+    workx = N.empty(2*len(x)-1,dtype=N.double)
+    worky = N.empty_like(workx)
     M = _lib.poly_clip(len(x),x,y,xleft,xright,ytop,ybottom,workx,worky)
     workx[M] = workx[0]
     worky[M] = worky[0]
@@ -193,10 +192,10 @@ def correlate(A,B,
     mode_column -- values outside boundaries, 'zero' or 'mirror'
 
     Returns a rows-by-columns array of correlation values.
-    
+
     """
 
-    A,B = _atype([A,B],[N.double,N.double])
+    A,B = atype([A,B],[N.double,N.double])
     assert A.ndim == 2 and B.ndim == 2, "Input arrays must be two dimensional"
 
     A_r,A_c = A.shape
@@ -210,7 +209,7 @@ def correlate(A,B,
 
     modes = {'zero': 0,
              'mirror': 1}
-        
+
     output = N.empty((rows,columns),dtype=N.double)
     _lib.correlate(A_r, A_c, A,
                    B_r, B_c, B,
@@ -220,7 +219,8 @@ def correlate(A,B,
 
     return output
 
-def interp_bilinear(grey_image,transform_coords_r,transform_coords_c,
+def interp_bilinear(grey_image,
+                    transform_coords_r=None,transform_coords_c=None,
                     mode='N',cval=0,output=None):
     """Calculate values at given coordinates using bi-linear interpolation.
 
@@ -249,23 +249,36 @@ def interp_bilinear(grey_image,transform_coords_r,transform_coords_c,
     An image of shape transform_coords_r.shape and type N.uint8.
 
     """
-    grey_image, = _atype(grey_image,N.uint8)
-    transform_coords_r,transform_coords_c = _atype([transform_coords_r,transform_coords_c],
+    grey_image, = atype(grey_image,N.uint8)
+    transform_coords_r,transform_coords_c = atype([transform_coords_r,transform_coords_c],
                                                    [N.double,N.double])
     assert grey_image.ndim == 2, "Input image must be 2-dimensional"
     assert transform_coords_r.ndim == 2 and transform_coords_c.ndim == 2, \
            "Transform coordinates must be 2-dimensional"
     if output is None:
-        output = N.empty(transform_coords_r.shape,dtype=N.uint8)        
+        output = N.empty(transform_coords_r.shape,dtype=N.uint8)
     else:
-        output, = _atype(output,N.uint8)
-        output.shape = transform_coords_r.shape        
+        output, = atype(output,N.uint8)
+        output.shape = transform_coords_r.shape
 
     rows,columns = grey_image.shape
     tf_rows,tf_columns = transform_coords_r.shape
     _lib.interp_bilinear(rows,columns,grey_image,
                          tf_rows,tf_columns,
                          transform_coords_r,transform_coords_c,
                          mode[0],cval,output)
-    
+
     return output
+
+def interp_transf_polygon(grey_image,transform,oshape=None):
+    grey_image,transform = atype([grey_image,transform],[N.uint8,N.double])
+    ishape = grey_image.shape
+    if oshape is None:
+        oshape = ishape
+
+    out = N.empty(oshape,dtype=SC.ftype)
+    _lib.interp_transf_polygon(ishape[0],ishape[1],grey_image,
+                               oshape[0],oshape[1],out,
+                               transform)
+    return out
+