Python implementation of Prasad DK, et al.'s parameter-independent line-fitting method

prasad_linefit.py

# a more or less line-for-line transcription of the source available here:
# https://docs.google.com/open?id=0B10RxHxW3I92dG9SU0pNMV84alk
# which is described in this paper:
# http://ieeexplore.ieee.org/document/6166585/
# D. K. Prasad, C. Quek, M. K. H. Leung and S. Y. Cho, "A parameter independent line fitting method," The First Asian Conference on Pattern Recognition, Beijing, 2011, pp. 441-445.
# doi: 10.1109/ACPR.2011.6166585
#
#
# Operation is pretty simple, just pass an **ordered** Nx2 array of x/y coordinates, get line segments.
# A convenience ordering function is provided free of charge ;)
import numpy as np

# needed imports if you want to use the order_edges function :)
#from skimage import morphology
#from scipy.spatial import distance
#from collections import deque as dq

def prasad_lines(edge):
    # edge should be a list of ordered coordinates
    # all credit to http://ieeexplore.ieee.org/document/6166585/
    # adapted from MATLAB scripts here: https://docs.google.com/open?id=0B10RxHxW3I92dG9SU0pNMV84alk
    # don't expect a lot of commenting from me here,
    # I don't claim to *understand* it, I just transcribed
    '''
    % This function takes each array of edgepoints in edgelist, finds the
    % size and position of the optimal deviation due to digitization
    % (Prasad, D.K., et. al. 2011, ACPR) from the line that joins the
    % endpoints, if the maximum deviation exceeds the optimal deviation due to
    % digitization then the
    % edge is shortened to the point of maximum deviation and the test is
    % repeated.  In this manner each edge is broken down to line segments,
    % each of which adhere to the original data with the specified tolerance.
    %
    % Copyright (c) 2012 Dilip K. Prasad
    % School of Computer Engineering
    % Nanyang Technological University, Singapore
    % http://www.ntu.edu.sg/
    %
    % Permission is hereby granted, free of charge, to any person obtaining a copy
    % of this software and associated documentation files (the "Software"), to deal
    % in the Software with restriction for its use for research purpose only, subject to the following conditions:
    %
    % The above copyright notice and this permission notice shall be included in
    % all copies or substantial portions of the Software.
    %
    % The Software is provided "as is", without warranty of any kind.
    %
    % Please cite the following work if this program is used
    % [1] Dilip K. Prasad, Chai Quek, Maylor K.H. Leung, and Siu-Yeung Cho
    %       “A parameter independent line fitting method,” 1st Asian Conference
    %       on Pattern Recognition (ACPR 2011), Beijing, China, 28-30 Nov, 2011.
    % [2] Dilip K. Prasad and Maylor K. H. Leung, "Polygonal Representation of
    %       Digital Curves," Digital Image Processing, Stefan G. Stanciu (Ed.),
    %       InTech, 2012.

    % April 2012 (Dilip K. Prasad) - Changes are made to adapt the line fitting method based on optimal deviation due to digitization
    '''

    x = edge[:,0]
    y = edge[:,1]

    first = 0
    last = len(edge)-1

    seglist = []
    seglist.append([x[0], y[0]])

    D = []
    precision = []
    reliability = []
    sum_dev = []
    D_temp = []


    while first<last:

        mdev_results = prasad_maxlinedev(x[first:last], y[first:last])
        print(mdev_results['d_max'])
        print(mdev_results['del_tol_max'])

        while mdev_results['d_max'] > mdev_results['del_tol_max']:
            print(last)
            last = mdev_results['index_d_max']+first
            print(last)
            mdev_results = prasad_maxlinedev(x[first:last], y[first:last])

        D.append(mdev_results['S_max'])
        seglist.append([x[last], y[last]])
        precision.append(mdev_results['precision'])
        reliability.append(mdev_results['reliability'])
        sum_dev.append(mdev_results['sum_dev'])

        first = last
        last = len(x)-1

    precision_edge = np.mean(precision)
    reliability_edge = np.sum(sum_dev)/np.sum(D)

    return seglist, precision_edge, reliability_edge



def prasad_maxlinedev(x, y):
    # all credit to http://ieeexplore.ieee.org/document/6166585/
    # adapted from MATLAB scripts here: https://docs.google.com/open?id=0B10RxHxW3I92dG9SU0pNMV84alk
    
    '''
    % This function takes each array of edgepoints in edgelist, finds the
    % size and position of the optimal deviation due to digitization (Prasad, D.K., et. al. 2011, ACPR) from the line that joins the
    % endpoints, if the maximum deviation exceeds the optimal deviation due to digitization then the
    % edge is shortened to the point of maximum deviation and the test is
    % repeated.  In this manner each edge is broken down to line segments,
    % each of which adhere to the original data with the specified tolerance.
    '''

    x = x.astype(np.float)
    y = y.astype(np.float)

    results = {}

    first = 0
    last = len(x)-1

    X = np.array([[x[0], y[0]], [x[last], y[last]]])
    A = np.array([
        [(y[0]-y[last]) / (y[0]*x[last] - y[last]*x[0])],
        [(x[0]-x[last]) / (x[0]*y[last] - x[last]*y[0])]
    ])

    if np.isnan(A[0]) and np.isnan(A[1]):
        devmat = np.column_stack((x-x[first], y-y[first])) ** 2
        dev = np.abs(np.sqrt(np.sum(devmat, axis=1)))
    elif np.isinf(A[0]) and np.isinf(A[1]):
        c = x[0]/y[0]
        devmat = np.column_stack((
            x[:]/np.sqrt(1+c**2),
            -c*y[:]/np.sqrt(1+c**2)
        ))
        dev = np.abs(np.sum(devmat, axis=1))
    else:
        devmat = np.column_stack((x, y))
        dev = np.abs(np.matmul(devmat, A)-1.)/np.sqrt(np.sum(A**2))

    results['d_max'] = np.max(dev)
    results['index_d_max'] = np.argmax(dev)
    results['precision'] = np.linalg.norm(dev, ord=2)/np.sqrt(float(last))
    s_mat = np.column_stack((x-x[first], y-y[first])) ** 2
    results['S_max'] = np.max(np.sqrt(np.sum(s_mat, axis=1)))
    results['reliability'] = np.sum(dev)/results['S_max']
    results['sum_dev'] = np.sum(dev)
    results['del_phi_max'] = prasad_digital_error(results['S_max'])
    results['del_tol_max'] = np.tan((results['del_phi_max']*results['S_max']))
    return results

def prasad_digital_error(ss):
    # all credit to http://ieeexplore.ieee.org/document/6166585/
    # adapted from MATLAB scripts here: https://docs.google.com/open?id=0B10RxHxW3I92dG9SU0pNMV84alk
    
    '''
    % This function computes the analytical error bound of the slope of a
    % digital line. See [1] for the derivation.
    %'''

    phii = np.arange(0, np.pi*2, np.pi / 360)

    s, phi = np.meshgrid(ss, phii)

    term1 = []

    term1.append(np.abs(np.cos(phi)))
    term1.append(np.abs(np.sin(phi)))
    term1.append(np.abs(np.sin(phi) + np.cos(phi)))
    term1.append(np.abs(np.sin(phi) - np.cos(phi)))

    term1.append(np.abs(np.cos(phi)))
    term1.append(np.abs(np.sin(phi)))
    term1.append(np.abs(np.sin(phi) + np.cos(phi)))
    term1.append(np.abs(np.sin(phi) - np.cos(phi)))

    tt2 = []
    tt2.append((np.sin(phi))/ s)
    tt2.append((np.cos(phi))/ s)
    tt2.append((np.sin(phi) - np.cos(phi))/ s)
    tt2.append((np.sin(phi) + np.cos(phi))/ s)

    tt2.append(-(np.sin(phi))/ s)
    tt2.append(-(np.cos(phi))/ s)
    tt2.append(-(np.sin(phi) - np.cos(phi))/ s)
    tt2.append(-(np.sin(phi) + np.cos(phi))/ s)


    term2 = []
    term2.append(s* (1 - tt2[0] + tt2[0]**2))
    term2.append(s* (1 - tt2[1] + tt2[1]**2))
    term2.append(s* (1 - tt2[2] + tt2[2]**2))
    term2.append(s* (1 - tt2[3] + tt2[3]**2))

    term2.append(s* (1 - tt2[4] + tt2[4]**2))
    term2.append(s* (1 - tt2[5] + tt2[5]**2))
    term2.append(s* (1 - tt2[6] + tt2[6]**2))
    term2.append(s* (1 - tt2[7] + tt2[7]**2))


    case_value = []
    for c_i in range(8):
        ss = s[:,0]
        t1 = term1[c_i]
        t2 = term2[c_i]
        case_value.append((1/ ss ** 2) * t1 * t2)

    return np.max(case_value)

  
def order_points(edge_points):
    # convert edge masks to ordered x-y coords

    if isinstance(edge_points, list):
        edge_points = np.array(edge_points)

    if edge_points.shape[1] > 2:
        # starting w/ imagelike thing, get the points
        edge_points = np.column_stack(np.where(edge_points))

    dists = distance.squareform(distance.pdist(edge_points))
    # make tri...nary connectedness graph, max dist is ~1.4 for diagonal pix
    # convert to 3 and 2 so singly-connected points are always <4
    dists[dists > 1.5] = 0
    dists[dists >1]    = 3
    dists[dists == 1]  = 2

    # check if we have easy edges
    # any point w/ sum of connectivity weights <4 can only have single connection
    dists_sum = np.sum(dists, axis=1)

    ends = np.where(dists_sum<4)[0]
    if len(ends)>0:
        pt_i = ends[0]
        first_i = ends[0]
        got_end = True
    else:
        # but if the edge has forex. a horiz & diag neighbor that'll fail
        # so just pick zero and get going.
        pt_i = 0
        first_i = 0
        got_end = False

    # walk through our dist graph, gathering points as we go
    inds = range(len(edge_points))
    new_pts = dq()
    forwards = True
    
    # this tight lil bundle will get reused a bit...
    # we are making a new list of points, and don't want to double-count points
    # but we can't pop from edge_points directly, because then the indices from the
    # dist mat will be wrong. Instead we make a list of all the indices and pop
    # from that. But since popping will change the index/value parity, we
    # have to double index inds.pop(inds.index(etc.))

    new_pts.append(edge_points[inds.pop(inds.index(pt_i))])
    
    while True:
        # get dict of connected points and distances
        # filtered by whether the index hasn't been added yet
        connected_pts = {k: dists[pt_i,k] for k in np.where(dists[pt_i,:])[0] if k in inds}

        # if we get nothing, we're either done or we have to go back to the first pt
        if len(connected_pts) == 0:
            if got_end:
                # still have points left, go back to first and go backwards
                # not conditioning on len(inds) because we might drop 
                # some degenerate diag points along the way
                # if we accidentally started at an end it won't hurt much
                pt_i = first_i
                forwards = False
                got_end = False
                continue
            else:
                # got to the end lets get outta here
                break

        # find point with min distance (take horiz/vert points before diags)
        pt_i = min(connected_pts, key=connected_pts.get)
        if forwards:
            new_pts.append(edge_points[inds.pop(inds.index(pt_i))])
        else:
            new_pts.appendleft(edge_points[inds.pop(inds.index(pt_i))])

    return np.array(new_pts)