Adapative iteration code that is running v slowly. The main function is run_all(), which then calls adaptive_iteration() which then calls the other functions.

gistfile1.txt

def make_ellipse(x, x0, y, y0, theta, a, b):
    c = np.cos(theta)
    s = np.sin(theta)
    a2 = a**2
    b2 = b**2
    xnew = x - x0
    ynew = y - y0
    ellipse = (xnew * c + ynew * s)**2/a2 + (xnew * s - ynew * c)**2/b2 <= 1

    return ellipse

def next_location(image, x, y):
    neighbours = image[y-1:y+2, x-1:x+2]
    if neighbours.shape != (3, 3):
        print "Edge? ",
        return (x, y)
    neighbours[1,1] = 0
    #print neighbours
    next_loc = np.where(neighbours == 1)
    if len(next_loc[0]) == 0:
        # If there are no new points then we've reached
        # the end of the line, so continue giving this point
        y_offset = 1
        x_offset = 1
        print "Error?"
    else:
        y_offset = next_loc[0][0]
        x_offset = next_loc[1][0]
    
    #print "Offsets: %d, %d" % (x_offset, y_offset)
    
    next_loc_x = (x-1)+x_offset
    next_loc_y = (y-1)+y_offset
    return np.array([next_loc_x, next_loc_y])


def walk_line(input_image, start_x, start_y, n):
    x = start_x
    y = start_y
    im = input_image.copy()
    for i in range(n):
        res = next_location(im, x, y)
        temp_x = res[0]
        temp_y = res[1]
        
        if temp_y >= input_image.shape[1] or temp_x >= input_image.shape[0]:
            print "Skipping because outside of array"
            return (x, y)
        else:
            x = temp_x
            y = temp_y
            im[x, y] = 0
        
    return (x, y)

def get_edges(image):
    edges, mag = mycanny(image, sigma=2.0, low_threshold=0.2, high_threshold=0.95)
    print "Found %d edges." % (np.sum(edges))
    edges = skeletonize(edges)
    labelled_edges, n = ndimage.measurements.label(edges, np.ones((3,3)))
    return edges, labelled_edges

def get_endpoints(image):
    footprint = np.array([[1,1,1],
                          [1,0,1],
                          [1,1,1]])
    endpoints =  np.bitwise_and(convolve2d(image, footprint, mode='same') == 1, image)
    return endpoints, np.where(endpoints == 1)

def adaptive_iteration(input_edges, lab_edges, s, a, b, prune_threshold):
    edges = input_edges.copy()
    endpoints_image, (endpoints_x, endpoints_y) = get_endpoints(edges)

    x, y = np.meshgrid(np.arange(-1*input_edges.shape[1]/2,input_edges.shape[1]/2), np.arange(-1*input_edges.shape[0]/2,input_edges.shape[0]/2))

    for i in xrange(len(endpoints_x)):
        print "%.2f" % (i/float(len(endpoints_x)) * 100),
        endpoint_x = endpoints_x[i]
        endpoint_y = endpoints_y[i]

        walked_y, walked_x = walk_line(edges, endpoint_y, endpoint_x, s)

        if walked_x < endpoint_x:
            left_x = walked_x
            left_y = walked_y
            right_x = endpoint_x
            right_y = endpoint_y
        else:
            left_x = endpoint_x
            left_y = endpoint_y
            right_x = walked_x
            right_y = walked_y
            
        slope = (right_x - left_x) / float(right_y - left_y)
        theta = np.arctan(slope)

        x0 = (endpoint_y - input_edges.shape[1]/2)
        y0 = (endpoint_x - input_edges.shape[0]/2)
        ellipse = make_ellipse(x, x0, y, y0, theta, a, b)

        edges = np.bitwise_or(edges, ellipse)
    
    skel = skeletonize(edges)
    return skel

def run_all():
  padhot = np.pad(hot, 1, mode='constant')

  edges, lab_edges = get_edges(padhot)
  orig_edges = edges.copy()
  
  result = edges.copy()

  result = adaptive_iteration(result, lab_edges, 3, 5, 3, 3)