public class AlgorithmHoughCardioid
extends AlgorithmBase

This Hough transform uses (xi, yi) points in the original image space to generate theta0, a0 points in the Hough transform. xc and yc of the cusp are obtained from finding the point of maximum curvature. Hough space is used to check the (xcdim * ycdim) possibilities of xc-xchalf to xc + xchalf, yc - ychalf to yc + ychalf, where xchalf = (xcdim - 1)/2 and ychalf = (ycdim - 1)/2. This Hough transform module only works with binary images. Before it is used the user must compute the gradient of an image and threshold it to obtain a binary image. Noise removal and thinning should also be performed, if necessary, before this program is run. The user is asked for the number of theta0 bins, a0 bins, side points for curvature, and number of cardioids. The default size for theta0Num is 720 and the default size for a0Num is min(512, 2*max(image.getExtents()[0], image.getExtents()[1]). The default number of cardioids is 1. The program generates a Hough transform of the source image using the basic equations: theta = atan2(y - yc, x - xc) where xc and yc are the cusp points. Calculate d2 = sqrt((x - d1)**2 + (y - d2)**2)/(1 - cos(theta + theta0)) if theta != -theta0 In general: sqrt((x - xc)**2 + (y - yc)**2) = a*(1 - cos(theta + theta0)) ((x-xc)**2 + (y-yc)**2 - a*((x-xc)*cos*(theta0) - (y-yc)*sin(theta0)) = a*sqrt((x-xc)**2 + (y-yc)**2) dy/dx = (-2*(x-xc) + a*cos(theta0) + a*(x-xc)/sqrt((x-xc)**2 + (y-yc)**2))/ (2*(y-yc) - a*sin(theta0) - a*(y-yc)/sqrt((x-xc)**2 + (y-yc)**2)) x = xc + a*cos(theta)*(1 - cos(theta + theta0)) = xc + a*(-0.5*cos(theta0) + cos(theta) - 0.5*cos(2*theta + theta0)) y = yc + a*sin(theta)*(1 - cos(theta + theta0)) = yc + a*(0.5*sin(theta0) + sin(theta) -0.5*sin(2*theta + theta0)) dy/dx = dy/dtheta/dx/dtheta = (-cos(2*theta + theta0) + cos(theta))/(sin(2*theta + theta0) - sin(theta)) = tan((1/2)*(3*theta + theta0)) Every slope value occurs 3 times so we must use the second derivative to find which of the 3 theta values is correct. dy'/dtheta = -a*sin(theta) + 2a*sin(2*theta + theta0) d2y/dx2 = dy'/dtheta/dx/dtheta = (-sin(theta) + 2*sin(2*theta + theta0))/(-sin(theta) + sin(2*theta + theta0)) curvature(theta) = ((dx/dtheta)*(d2y/dtheta2) - (dy/dtheta)*(d2x/dtheta2))/((dx/dtheta)**2 + (dy/dtheta)**2)**1.5 d2x/dtheta2 = -a*cos(theta) + 2*a*cos(2*theta + theta0) d2y/dtheta2 = -a*sin(theta) + 2*a*sin(2*theta + theta0) curvature(theta) = 3/(2*sqrt(2)*a*sqrt(1 - cos(theta+theta0)) so the curvature is infinite at theta = -theta0. All cusp chords are of length 2 * a. The tangents to the endpoints of a cusp chord are perpindicular. If 3 points have parallel tangents, the lines from the cusp to these 3 points make equal angles of 2*PI/3 at the cusp. For cusp on the left: sqrt((x - xc)**2 + (y - yc)**2) = a*(1 + cos(theta)). x = xc + (a/2)*(1 + 2*cos(theta) + cos(2*theta))) = xc + a*cos(theta)*(1 + cos(theta)) y = yc + (a/2)*(2*sin(theta) + sin(2*theta)) = yc + a*sin(theta)*(1 + cos(theta)) For cusp on the right: sqrt((x - xc)**2 + (y - yc)**2) = a*(1 - cos(theta)). x = xc + (a/2)*(-1 + 2*cos(theta) - cos(2*theta))) = xc + a*cos(theta)*(1 - cos(theta)) y = yc + (a/2)*(2*sin(theta) - sin(2*theta)) = yc + a*sin(theta)*(1 - cos(theta)) For cusp on top: sqrt((x - xc)**2 + (y - yc)**2) = a*(1 + sin(theta)). x = xc + (a/2)*(2*cos(theta) + sin(2*theta))) = xc + a*cos(theta)*(1 + sin(theta)) y = yc + (a/2)*(1 + 2*sin(theta) - cos(2*theta)) = yc + a*sin(theta)*(1 + sin(theta)) For cusp on bottom: sqrt((x - xc)**2 + (y - yc)**2) = a*(1 - sin(theta)). x = xc + (a/2)*(2*cos(theta) - sin(2*theta))) = xc + a*cos(theta)*(1 - sin(theta)) y = yc + (a/2)*(-1 + 2*sin(theta) + cos(2*theta)) = yc + a*sin(theta)*(1 - sin(theta)) The program finds the cardioids containing the largest number of points. The program produces a dialog which allows the user to select which cardioids should be drawn. The Hough transform for the entire image is generated a separate time to find each cardioid. Calculate theta = atan2(y - yc, x - xc). Calculate d3 = sqrt((x - xc)**2 + (y - yc)**2)/(1 - cos(theta + theta0)) if theta != -theta0 Don't calculate d2 if theta = -theta0. d2 goes from 0 to maxA = sqrt((xDim-1)**2 + (yDim-1)**2)/2.0. s2 is the dimension 2 scaling factor. s2 * (a0 - 1) = maxA. s2 = maxA/(a0 - 1) m = d2*(a0 - 1)/maxA. Only calculate the Hough transform for d2