Initial commit

circle_fit.py

@@ -0,0 +1,52 @@
+""" circle_fit -- Ben Gamari (2012)
+    Based on,
+      L. Maisonobe. "Finding the circle that best fits a set of points" (2007).
+      <http://www.spaceroots.org/documents/circle/circle-fitting.pdf>
+"""
+
+import numpy as np
+from numpy import sqrt, mean, sum
+import scipy.optimize
+
+def est_center(points):
+    sx = 0
+    sy = 0
+    q = 0
+    n = points.shape[0]
+    for i in xrange(0, n-2):
+        for j in xrange(i+1, n-1):
+            for k in xrange(j+1, n):
+                (xi,yi) = points[i]
+                (xj,yj) = points[j]
+                (xk,yk) = points[k]
+                d = (xk-xj)*(yj-yi) - (xj-xi)*(yk-yj)
+                if abs(d) < 1e-6: continue
+                xc = (yk-yj) * (xi**2+yi**2) \
+                   + (yi-yk) * (xj**2+yj**2) \
+                   + (yj-yi) * (xk**2+yk**2)
+
+                yc = (xk-xj) * (xi**2+yi**2) \
+                   + (xi-xk) * (xj**2+yj**2) \
+                   + (xj-xi) * (xk**2+yk**2)
+                sx  += +xc / 2 / d
+                sy  += -yc / 2 / d
+                q += 1
+    if q == 0:
+        raise RuntimeError("Not enough non-aligned points")
+    else:
+        return sx/q, sy/q
+
+def est_radius(points, center):
+    return mean(sqrt(sum((points - center[newaxis,:])**2, axis=0)))
+
+def fit_circle(points, c0=None):
+    ps = points
+    if c0 is None:
+        c0 = est_center(points)
+    print 'Estimated center', c0
+
+    d = lambda (cx, cy, r): sqrt((points[:,0]-cx)**2 + (points[:,1]-cy)**2) - r
+    out = scipy.optimize.leastsq(d, (c0[0], c0[1], 1))
+    xopt = out[0]
+    return ((xopt[0], xopt[1]), xopt[2])
+

gaussian_smooth.py

@@ -0,0 +1,22 @@
+import numpy as np
+import scipy.signal as signal
+
+def gauss_kern(size, sizey=None):
+    """ Returns a normalized 2D gauss kernel array for convolutions """
+    size = int(size)
+    if not sizey:
+        sizey = size
+    else:
+        sizey = int(sizey)
+        x, y = np.mgrid[-size:size+1, -sizey:sizey+1]
+        g = exp(-(x**2/float(size)+y**2/float(sizey)))
+    return g / g.sum()
+
+def blur_image(im, n, ny=None) :
+    """ blurs the image by convolving with a gaussian kernel of typical
+        size n. The optional keyword argument ny allows for a different
+        size in the y direction.
+    """
+    g = gauss_kern(n, sizey=ny)
+    improc = signal.convolve(im,g, mode='valid')
+    return improc

project.py

@@ -0,0 +1,59 @@
+import sys
+from PIL import Image
+import numpy as np
+from numpy import cos, sin, pi, r_, c_
+from scipy.interpolate import griddata
+from matplotlib import pyplot as pl
+from circle_fit import fit_circle
+
+def prompt_center(d):
+    points = []
+    fig = pl.figure()
+
+    def onclick(event):
+        points.append((event.xdata, event.ydata))
+        pl.plot([event.xdata], [event.ydata], 'w+')
+        fig.canvas.draw()
+
+    def key_press(event):
+        if event.key == "enter":
+            pl.close()
+
+    fig.canvas.mpl_connect('button_press_event', onclick)
+    fig.canvas.mpl_connect('key_press_event', key_press)
+    pl.imshow(np.log1p(d))
+    pl.autoscale(enable=False)
+    pl.show()
+    ((cx,cy), r) = fit_circle(np.array(points))
+    return (cx,cy)
+
+def cartesian_projection(d, center, r_min=None, r_max=None, nphi=300, nr=100):
+    r_max = min([center[0], center[1], d.shape[0]-center[0], d.shape[1]-center[1]] + ([r_max] if r_max is not None else []))
+    if r_min is None: r_min = 0
+    xs,ys = np.indices(d.shape)
+    rs, phis = np.meshgrid(np.linspace(r_min, r_max, nr), np.linspace(0, 2*pi, nphi))
+    sample_pts = (rs * cos(phis) + center[0], rs * sin(phis) + center[1])
+    samples = griddata((xs.flatten(),ys.flatten()), d.flatten(), sample_pts, method='linear')
+    return np.rec.fromarrays([rs, phis, samples], names='r,phi,i')
+
+img = Image.open(sys.argv[1]).convert('F')
+print "Value extrema", img.getextrema()
+img = np.array(img, dtype='f')
+print "Image size", img.shape
+center = np.array(img.shape) / 2
+
+center = prompt_center(img)
+print "Image center", center
+pl.figure()
+pl.imshow(img)
+pl.savefig('proj1.png')
+
+p = cartesian_projection(img, center)
+pl.figure()
+pl.imshow(p['i'], aspect='auto')
+pl.savefig('proj.png')
+
+pl.figure()
+p = p[np.logical_and(100 < p['r'], p['r'] < 400)]
+pl.errorbar(p['r'][0,:], np.mean(p['i'], axis=0), yerr=np.std(p['i'], axis=0))
+pl.savefig('i.png')