Initial commit

1 files changed, 1023 insertions, 0 deletions

diff --git a/pkg/images/lib/rgmerge.x b/pkg/images/lib/rgmerge.x
new file mode 100644
index 00000000..5e218bd0
--- /dev/null
+++ b/pkg/images/lib/rgmerge.x
@@ -0,0 +1,1023 @@
+include <mach.h>
+include <plset.h>
+include "xyxymatch.h"
+
+# RG_MATCH -- Compute the intersection of two lists using a pattern matching
+# algorithm. This algorithm is based on one developed by Edward Groth
+# 1986 A.J. 91, 1244. The algorithm matches pairs of coordinates from
+# two lists based on the triangles that can be formed from triplets of
+# points in each list. The algorithm is insensitive to coordinate translation,
+# rotation, magnification, or inversion and can tolerate distortions and
+# random errors.
+
+int procedure rg_match (xref, yref, nref, xin, yin, nin, reftri, reftrirat,
+	nreftri, nrmaxtri, nrefstars, intri, intrirat, nintri, ninmaxtri,
+	nliststars, tolerance, ptolerance, ratio, nreject)
+
+real	xref[ARB]		#I the reference x coordinates
+real	yref[ARB]		#I the reference y coordinates
+int	nref			#I the number of reference coordinates
+real	xin[ARB]		#I the input x coordinates
+real	yin[ARB]		#I the input y coordinates
+int	nin			#I the number of input coordinates
+int	reftri[nrmaxtri,ARB]	#U list of reference triangles
+real	reftrirat[nrmaxtri,ARB]	#U list of reference triangle parameters
+int	nreftri			#U number of reference triangles
+int	nrmaxtri		#I maximum number of reference triangles
+int	nrefstars		#I the number of reference stars
+int	intri[ninmaxtri,ARB]	#U list of input triangles
+real	intrirat[ninmaxtri,ARB]	#U list of input triangle parameters
+int	nintri			#U number of input triangles
+int	ninmaxtri		#I maximum number of input triangles
+int	nliststars		#I the number of input stars
+real	tolerance		#I the reference triangles matching tolerance
+real	ptolerance		#I the input triangles matching tolerance
+real	ratio			#I the maximum ratio of triangle sides
+int	nreject			#I maximum number of rejection iterations
+
+int	i, nmerge, nkeep, nmatch, ncheck
+pointer	sp, rindex, lindex
+int	rg_tmerge(), rg_treject(), rg_tvote(), rg_triangle
+
+begin
+	# Match the triangles in the input list to those in the reference list.
+	if (nreftri < nintri)
+	    nmerge = rg_tmerge (reftri, reftrirat, nreftri, nrmaxtri, intri,
+	        intrirat, nintri, ninmaxtri)
+	else
+	    nmerge = rg_tmerge (intri, intrirat, nintri, ninmaxtri, reftri,
+	        reftrirat, nreftri, nrmaxtri)
+	if (nmerge <= 0)
+	    return (0)
+
+	# Perform the rejection cycle.
+	nkeep = rg_treject (reftri, reftrirat, nreftri, nrmaxtri,
+	    intri, intrirat, nintri, ninmaxtri, nmerge, nreject)
+	if (nkeep <= 0)
+	    return (0)
+
+	# Match the coordinates.
+	nmatch = rg_tvote (reftri, nrmaxtri, nrefstars, intri, ninmaxtri,
+	    nliststars, nkeep)
+	if (nmatch <= 0)
+	    return (0)
+	else if (nmatch <= 3 && nkeep < nmerge)
+	    return (0)
+
+	# If all the coordinates were not matched then make another pass
+	# through the triangles matching algorithm. If the number of
+	# matches decreases as a result of this then all the matches were
+	# not true matches and declare the list unmatched.
+	if (nmatch < min (nref, nin) && nmatch > 2) {
+
+	    # Find the indices of the matched points.
+	    call smark (sp)
+	    call salloc (rindex, nmatch, TY_INT)
+	    call salloc (lindex, nmatch, TY_INT)
+	    do i = 1, nmatch {
+		Memi[rindex+i-1] = reftri[i,RG_MATCH]
+		Memi[lindex+i-1] = intri[i,RG_MATCH]
+	    }
+
+	    # Recompute the triangles.
+	    nreftri = rg_triangle (xref, yref, Memi[rindex], nmatch, reftri,
+		reftrirat, nrmaxtri, nrefstars, tolerance, ratio)
+	    nintri = rg_triangle (xin, yin, Memi[lindex], nmatch, intri,
+		intrirat, ninmaxtri, nliststars, ptolerance, ratio)
+
+	    # Rematch the triangles.
+	    if (nreftri < nintri)
+	        nmerge = rg_tmerge (reftri, reftrirat, nreftri, nrmaxtri, intri,
+	            intrirat, nintri, ninmaxtri)
+	    else
+	        nmerge = rg_tmerge (intri, intrirat, nintri, ninmaxtri, reftri,
+	            reftrirat, nreftri, nrmaxtri)
+
+	    # Reperform the rejection cycle.
+	    if (nmerge > 0)
+	        nkeep = rg_treject (reftri, reftrirat, nreftri, nrmaxtri,
+	            intri, intrirat, nintri, ninmaxtri, nmerge, nreject)
+
+	    # Reperform the vote.
+	    if (nkeep > 0) {
+	        ncheck = rg_tvote (reftri, nrmaxtri, nrefstars, intri,
+	            ninmaxtri, nliststars, nkeep)
+		if (ncheck <= 3 && nkeep < nmerge)
+		    ncheck = 0
+	    } else
+		ncheck = 0
+
+	    if (ncheck < nmatch)
+		nmatch = 0
+	    else
+		nmatch = ncheck
+
+	    call sfree (sp)
+	}
+
+	return (nmatch)
+end
+ 
+
+# RG_TRIANGLE -- Construct all the the possible triangles from
+# an input coordinate list. The triangles are constructed in such a way
+# that the shortest side of the triangle lies between vertices 1 and 2 and the
+# longest side between vertices 1 and 3. The parameters of each triangle
+# including the log of the perimeter, the ratio of the longest to shortest
+# side, the cosine of the angle at vertex 1, the tolerances in the ratio
+# and cosine and the sense of the triangle (clockwise or anti-clockwise)
+# are also computed. Triangles with a ratio greater than maxratio are
+# rejected as are triangles with vertices closer together than tolerance.
+
+int procedure rg_triangle (xref, yref, refindex, nrefstars, reftri, tripar,
+	nmaxtri, maxnpts, tolerance, maxratio)
+
+real	xref[ARB]		#I x reference coordinates
+real	yref[ARB]		#I y reference coordinates
+int	refindex[ARB]		#I the reference list sort index
+int	nrefstars		#I number of reference stars
+int	reftri[nmaxtri,ARB]	#O reference triangles
+real	tripar[nmaxtri,ARB]	#O triangle parameters
+int	nmaxtri			#I maximum number of triangles
+int	maxnpts			#I the maximum number of points
+real	tolerance		#I matching tolerance
+real	maxratio		#I maximum ratio of triangle sides
+
+int	i, j, k, nsample, npts, ntri
+real	rij, rjk, rki, dx1, dy1, dx2, dy2, dx3, dy3, r1, r2sq, r2, r3sq, r3
+real	ratio, cosc, cosc2, sinc2, tol2, tol
+
+begin
+	# Create the tolerance.
+	tol2 = tolerance ** 2
+	nsample = max (1, nrefstars / maxnpts) 
+	npts = min (nrefstars, nsample * maxnpts)
+
+	# Construct the triangles.
+	ntri = 1
+	do i = 1, npts - 2 * nsample, nsample {
+	    do j = i + nsample, npts - nsample, nsample {
+		do k = j + nsample, npts, nsample {
+
+		    # Compute the lengths of the three sides of the triangle,
+		    # eliminating triangles with sides that are less than
+		    # tolerance.
+		    rij = (xref[refindex[i]] - xref[refindex[j]]) ** 2 +
+		        (yref[refindex[i]] - yref[refindex[j]]) ** 2
+		    if (rij <= tol2)
+		        next
+		    rjk = (xref[refindex[j]] - xref[refindex[k]]) ** 2 +
+		        (yref[refindex[j]] - yref[refindex[k]]) ** 2
+		    if (rjk <= tol2)
+			next
+		    rki = (xref[refindex[k]] - xref[refindex[i]]) ** 2 +
+		        (yref[refindex[k]] - yref[refindex[i]]) ** 2
+		    if (rki <= tol2)
+			next
+
+		    # Order the vertices with the shortest side of the triangle
+		    # between vertices 1 and 2 and the intermediate side between
+		    # vertices 2 and 3.
+		    reftri[ntri,RG_INDEX] = ntri
+		    if (rij <= rjk) {
+	                if (rki <= rij) {
+			    reftri[ntri,RG_X1] = refindex[k]
+			    reftri[ntri,RG_X2] = refindex[i]
+			    reftri[ntri,RG_X3] = refindex[j]
+	    		} else if (rki >= rjk) {
+			    reftri[ntri,RG_X1] = refindex[i]
+			    reftri[ntri,RG_X2] = refindex[j]
+			    reftri[ntri,RG_X3] = refindex[k]
+	    		} else {
+			    reftri[ntri,RG_X1] = refindex[j]
+			    reftri[ntri,RG_X2] = refindex[i]
+			    reftri[ntri,RG_X3] = refindex[k]
+	    	        }
+		    } else {
+	    	        if (rki <= rjk) {
+			    reftri[ntri,RG_X1] = refindex[i]
+			    reftri[ntri,RG_X2] = refindex[k]
+			    reftri[ntri,RG_X3] = refindex[j]
+	    	        } else if (rki >= rij) {
+			    reftri[ntri,RG_X1] = refindex[k]
+			    reftri[ntri,RG_X2] = refindex[j]
+			    reftri[ntri,RG_X3] = refindex[i]
+	    	        } else {
+			    reftri[ntri,RG_X1] = refindex[j]
+			    reftri[ntri,RG_X2] = refindex[k]
+			    reftri[ntri,RG_X3] = refindex[i]
+	    	        }
+		    }
+
+		    # Compute the lengths of the sides.
+		    dx1 = xref[reftri[ntri,RG_X3]] - xref[reftri[ntri,RG_X2]]
+		    dy1 = yref[reftri[ntri,RG_X3]] - yref[reftri[ntri,RG_X2]]
+		    dx2 = xref[reftri[ntri,RG_X2]] - xref[reftri[ntri,RG_X1]]
+		    dy2 = yref[reftri[ntri,RG_X2]] - yref[reftri[ntri,RG_X1]]
+		    dx3 = xref[reftri[ntri,RG_X3]] - xref[reftri[ntri,RG_X1]]
+		    dy3 = yref[reftri[ntri,RG_X3]] - yref[reftri[ntri,RG_X1]]
+
+		    # Compute the ratio of the longest side of the triangle
+		    # to the shortest side.
+		    r1 = sqrt (dx1 ** 2 + dy1 ** 2)
+		    r2sq = dx2 ** 2 + dy2 ** 2
+		    r2 = sqrt (r2sq)
+		    r3sq = dx3 ** 2 + dy3 ** 2
+		    r3 = sqrt (r3sq)
+		    if (r2 <= 0.)
+	    	        next
+		    ratio = r3 / r2
+		    if (ratio > maxratio)
+	    		next
+
+		    # Compute the cos, cos ** 2 and sin ** 2 of the angle at 
+		    # vertex 1.
+		    cosc = (dx3 * dx2 + dy3 * dy2) / (r3 * r2)
+		    cosc2 = max (0.0, min (1.0, cosc * cosc))
+		    sinc2 = max (0.0, min (1.0, 1.0 - cosc2))
+
+		    # Determine whether the triangles vertices are arranged
+		    # clockwise of anticlockwise.
+		    if ((dx2 * dy1 - dy2 * dx1) > 0.0)
+			reftri[ntri,RG_CC] = YES
+		    else
+			reftri[ntri,RG_CC] = NO
+
+		    # Compute the tolerances.
+		    tol = (1.0 / r3sq - cosc / (r3 * r2) + 1.0 / r2sq)
+		    tripar[ntri,RG_TOLR] = 2.0 * ratio ** 2 * tol2 * tol 
+		    tripar[ntri,RG_TOLC] = 2.0 * sinc2 * tol2 * tol + 3.0 *
+		        cosc2 * tol2 ** 2 * tol * tol
+
+		    # Compute the perimeter.
+		    tripar[ntri,RG_LOGP] = log (r1 + r2 + r3)
+		    tripar[ntri,RG_RATIO] = ratio
+		    tripar[ntri,RG_COS1] = cosc
+
+		    ntri = ntri + 1
+		}
+	    }
+	}
+
+	ntri = ntri - 1
+
+	# Sort the triangles in increasing order of ratio.
+	call rg_qsortr (tripar[1,RG_RATIO], reftri[1,RG_INDEX],
+	    reftri[1,RG_INDEX], ntri)
+
+	return (ntri)
+end
+
+
+# RG_TMERGE -- Compute the intersection of two sorted files of triangles
+# using the tolerance parameter.
+
+int procedure rg_tmerge (reftri, rtripar, nrtri, nmrtri, listri, ltripar,
+	nltri, nmltri)
+
+int	reftri[nmrtri,ARB]	#U list of reference triangles
+real	rtripar[nmrtri,ARB]	#I reference triangle parameters
+int	nrtri			#I number of reference triangles
+int	nmrtri			#I maximum number of reference triangles
+int	listri[nmltri,ARB]	#U list of reference triangles
+real	ltripar[nmltri,ARB]	#I reference triangle parameters
+int	nltri			#I number of reference triangles
+int	nmltri			#I maximum number of reference triangles
+
+int	rp, blp, lp, ninter, rindex, lindex, mindex
+real	rmaxtol, lmaxtol, maxtol, dr, dr2, mdr2, dcos2, mdcos2, dtolr, dtolc 
+
+begin
+	# Find the maximum tolerance for each list.
+	call alimr (rtripar[1,RG_TOLR], nrtri, maxtol, rmaxtol)
+	call alimr (ltripar[1,RG_TOLR], nltri, maxtol, lmaxtol)
+	maxtol = sqrt (rmaxtol + lmaxtol)
+
+	# Define the beginning of the search range for each triangle.
+	blp = 1
+
+	# Loop over all the triangles in the reference list.
+	ninter = 0
+	for (rp = 1; rp <= nrtri; rp = rp + 1) {
+
+	    # Get the index for the reference triangle.
+	    rindex = reftri[rp,RG_INDEX]
+
+	    # Move to the first triangle in the input list that satisfies the
+	    # ratio tolerance requirement.
+	    for  (; blp <= nltri; blp = blp + 1) {
+		lindex = listri[blp,RG_INDEX]
+		dr = rtripar[rindex,RG_RATIO] - ltripar[lindex,RG_RATIO]
+		if (dr <= maxtol)
+		    break
+	    }
+
+	    # If the beginning of the search range becomes greater than
+	    # the length of the list then there is no match.
+	    if (blp > nltri)
+		break
+
+	    # If the first triangle in the list is past the tolerance
+	    # limit skip to the next reference triangle
+	    if (dr < -maxtol)
+		next
+
+	    # Search through the appropriate range of triangles for the 
+	    # closest fit.
+
+	    # Initialize the tolerances.
+	    mindex = 0
+	    mdr2 = 0.5 * MAX_REAL
+	    mdcos2 = 0.5 * MAX_REAL
+
+	    for (lp = blp; lp <= nltri; lp = lp + 1) {
+
+		# Quit the loop if the next triangle is out of match range.
+		lindex = listri[lp,RG_INDEX]
+		dr = rtripar[rindex,RG_RATIO] - ltripar[lindex,RG_RATIO]
+		if (dr < -maxtol)
+		    break
+
+		# Compute the tolerances for the two triangles.
+		dr2 = dr * dr
+		dcos2 = (rtripar[rindex,RG_COS1] - ltripar[lindex,RG_COS1]) ** 2
+		dtolr = rtripar[rindex,RG_TOLR] + ltripar[lindex,RG_TOLR]
+		dtolc = rtripar[rindex,RG_TOLC] + ltripar[lindex,RG_TOLC] 
+
+		# Find the best of all possible matches.
+		if (dr2 <= dtolr && dcos2 <= dtolc) {
+		    if ((dr2 + dcos2) < (mdr2 + mdcos2)) {
+			mindex = lindex
+			mdr2 = dr2
+			mdcos2 = dcos2
+		    }
+		}
+
+	    }
+
+	    # Add the match to the list.
+	    if (mindex > 0) {
+		ninter = ninter + 1
+		reftri[ninter,RG_MATCH] = rindex
+		listri[ninter,RG_MATCH] = mindex
+	    }
+	}
+
+	return (ninter)
+end
+
+
+# RG_TREJECT -- Remove false matches from the list of matched triangles.
+
+int procedure rg_treject (reftri, rtripar, nrtri, nmrtri, listri, ltripar,
+	nltri, nmltri, nmatch, maxiter)
+
+int	reftri[nmrtri,ARB]	#U list of reference triangles
+real	rtripar[nmrtri,ARB]	#I reference triangle parameters
+int	nrtri			#I number of reference triangles
+int	nmrtri			#I maximum number of reference triangles
+int	listri[nmltri,ARB]	#U list of reference triangles
+real	ltripar[nmltri,ARB]	#I reference triangle parameters
+int	nltri			#I number of reference triangles
+int	nmltri			#I maximum number of reference triangles
+int	nmatch			#I initial number of matches
+int	maxiter			#I maximum number of rejection iterations
+
+double	dif, mode, sum, sumsq
+int	i, nrej, nplus, nminus, ntrue, nfalse, npts, ncount, niter, rindex
+int	lindex
+pointer	sp, adif
+real	sigma, factor, locut, hicut
+double	rg_moded()
+
+begin
+	call smark (sp)
+	call salloc (adif, nmatch, TY_DOUBLE)
+
+	# Accumulate the number of same sense and number of opposite sense
+	# matches as well as the log perimeter statistics.
+	sum = 0.0d0
+	sumsq = 0.0d0
+	nplus = 0
+	do i = 1, nmatch {
+	    rindex = reftri[i,RG_MATCH]
+	    lindex = listri[i,RG_MATCH]
+	    dif = (rtripar[rindex,RG_LOGP] - ltripar[lindex,RG_LOGP])
+	    Memd[adif+i-1] = dif
+	    sum = sum + dif
+	    sumsq = sumsq + dif * dif
+	    if (reftri[rindex,RG_CC] == listri[lindex,RG_CC])
+		nplus = nplus + 1
+	}
+	nminus = nmatch - nplus
+
+	# Compute the mean, mode, and sigma of the logP distribution,
+	ntrue = abs (nplus - nminus)
+	nfalse = nplus + nminus - ntrue
+	#mean = sum / nmatch
+	if (nmatch <= 1)
+	    sigma = 0.0
+	else
+	    sigma = (sumsq - (sum / nmatch) * sum) / (nmatch - 1)
+	if (sigma <= 0.0) {
+	    call sfree (sp)
+	    return (nmatch)
+	} else
+	    sigma = sqrt (sigma)
+	call asrtd (Memd[adif], Memd[adif], nmatch)
+        #if (mod (nmatch,2) == 1)
+            #median = Memd[adif+nmatch/2]
+        #else
+	    #median = (Memd[adif+nmatch/2] + Memd[adif+(nmatch-1)/2]) / 2.0d0
+	mode = rg_moded (Memd[adif], nmatch, 10, 1.0d0, 0.1d0 * sigma,
+	    0.01d0 * sigma)
+	if (nfalse > ntrue)
+	    factor = 1.0
+	else if ((0.1 * ntrue) > nfalse)
+	    factor = 3.0
+	else
+	    factor = 2.0
+
+	# Begin the rejection cycle.
+	npts = nmatch
+	niter = 0
+	repeat {
+
+	    ncount = 0
+	    nrej = 0
+	    locut = mode - factor * sigma
+	    hicut = mode + factor * sigma
+
+	    # Reject matched triangles which are too far from the mean logP.
+	    do i = 1, npts {
+		rindex = reftri[i,RG_MATCH]
+		lindex = listri[i,RG_MATCH]
+		dif = rtripar[rindex,RG_LOGP] - ltripar[lindex,RG_LOGP]
+		if (dif < locut || dif > hicut) {
+		    sum = sum - dif
+		    sumsq = sumsq - dif * dif
+	    	    if (reftri[rindex,RG_CC] == listri[lindex,RG_CC])
+			nplus = nplus - 1
+		    else
+			nminus = nminus - 1
+		    nrej = nrej + 1
+		} else {
+		    Memd[adif+ncount] = dif
+		    ncount = ncount + 1
+		    reftri[ncount,RG_MATCH] = rindex
+		    listri[ncount,RG_MATCH] = lindex
+		}
+	    }
+
+	    # No more points were rejected.
+	    npts = ncount
+	    if (nrej <= 0)
+		break
+
+	    # All the points were rejected.
+	    if (npts <= 0)
+		break
+
+	    # The rejection iteration limit was reached.
+	    niter = niter + 1
+	    if (niter >= maxiter)
+		break
+
+	    # Compute the new mean and sigma of the logP distribution.
+	    #mean = sum / npts
+	    if (npts <= 1)
+	        sigma = 0.0
+	    else
+	        sigma = (sumsq - (sum / npts) * sum) / (npts - 1)
+	    if (sigma <= 0.0)
+		break
+	    sigma = sqrt (sigma)
+	    call asrtd (Memd[adif], Memd[adif], npts)
+            #if (mod (npts,2) == 1)
+                #median = Memd[adif+npts/2]
+            #else
+		#median = (Memd[adif+npts/2] + Memd[adif+(npts-1)/2]) / 2.0d0
+	    mode = rg_moded (Memd[adif], npts, 10, 1.0d0, 0.10d0 * sigma,
+	        0.01d0 * sigma)
+
+	    # Recompute the ksigma rejection criterion based on the number of
+	    # same and opposite sense matches.
+	    ntrue = abs (nplus - nminus)
+	    nfalse = nplus + nminus - ntrue
+	    if (nfalse > ntrue)
+	        factor = 1.0
+	    else if ((0.1 * ntrue) > nfalse)
+	        factor = 3.0
+	    else
+	        factor = 2.0
+	}
+
+	# One last iteration to get rid of opposite sense of matches.
+	if (npts <= 0)
+	    npts = 0
+	else if (nplus > nminus) {
+	    ncount = 0
+	    do i = 1, npts {
+	        rindex = reftri[i,RG_MATCH]
+	        lindex = listri[i,RG_MATCH]
+	    	if (reftri[rindex,RG_CC] == listri[lindex,RG_CC]) {
+		    ncount = ncount + 1
+		    reftri[ncount,RG_MATCH] = rindex
+		    listri[ncount,RG_MATCH] = lindex
+		}
+	    }
+	    npts = ncount
+	} else {
+	    ncount = 0
+	    do i = 1, npts {
+	        rindex = reftri[i,RG_MATCH]
+	        lindex = listri[i,RG_MATCH]
+	    	if (reftri[rindex,RG_CC] != listri[lindex,RG_CC]) {
+		    ncount = ncount + 1
+		    reftri[ncount,RG_MATCH] = rindex
+		    listri[ncount,RG_MATCH] = lindex
+		}
+	    }
+	    npts = ncount
+	}
+
+	call sfree (sp)
+	return (npts)
+end
+
+
+# RG_TVOTE -- Count the number a times a particular pair of
+# coordinates is matched in the set of matched triangles. If a particular
+# pair of points occurs in many triangles it is much more likely to be
+# a true match than if it occurs in very few. Since this vote array
+# may be quite sparsely occupied, use the PLIO package to store and
+# maintain the list.
+
+int procedure rg_tvote (reftri, nmrtri, nrefstars, listri, nmltri, nliststars,
+	nmatch)
+
+int	reftri[nmrtri,ARB]		#U reference triangles
+int	nmrtri				#I maximum number of reference triangles
+int	nrefstars			#I number of reference stars
+int	listri[nmltri,ARB]		#U input list triangles
+int	nmltri				#I maximum number of list triangles
+int	nliststars			#I number of list stars
+int	nmatch				#I number of match triangles
+
+int	i, j, rp, lp, vp, pixval, tminvote, tmaxvote, minvote, maxvote, hmaxvote
+int	ninter, axes[2], laxes[2], pvp
+pointer	sp, vote, vindex, pl, lmatch, rmatch
+bool	pl_linenotempty()
+pointer	pl_create()
+
+begin
+	# Open the pixel list.
+	axes[1] = nliststars
+	axes[2] = nrefstars
+	pl = pl_create (2, axes, 16)
+
+	# Acumulate the votes.
+	do i = 1, nmatch {
+	    rp = reftri[i,RG_MATCH]
+	    lp = listri[i,RG_MATCH]
+	    laxes[1] = listri[lp,RG_X1]
+	    laxes[2] = reftri[rp,RG_X1]
+	    if (! pl_linenotempty (pl, laxes))
+	        call pl_point (pl, laxes[1], laxes[2], PIX_SET + PIX_VALUE(1))
+	    else {
+		call pl_glpi (pl, laxes, pixval, 16, 1, PIX_SRC)
+		pixval = pixval + 1
+	        call pl_point (pl, laxes[1], laxes[2], PIX_SET +
+		    PIX_VALUE(pixval))
+	    }
+	    laxes[1] = listri[lp,RG_X2]
+	    laxes[2] = reftri[rp,RG_X2]
+	    if (! pl_linenotempty (pl, laxes))
+	        call pl_point (pl, laxes[1], laxes[2], PIX_SET + PIX_VALUE(1))
+	    else {
+		call pl_glpi (pl, laxes, pixval, 16, 1, PIX_SRC)
+		pixval = pixval + 1
+	        call pl_point (pl, laxes[1], laxes[2], PIX_SET +
+		    PIX_VALUE(pixval))
+	    }
+	    laxes[1] = listri[lp,RG_X3]
+	    laxes[2] = reftri[rp,RG_X3]
+	    if (! pl_linenotempty (pl, laxes))
+	        call pl_point (pl, laxes[1], laxes[2], PIX_SET + PIX_VALUE(1))
+	    else {
+		call pl_glpi (pl, laxes, pixval, 16, 1, PIX_SRC)
+		pixval = pixval + 1
+	        call pl_point (pl, laxes[1], laxes[2], PIX_SET +
+		    PIX_VALUE(pixval))
+	    }
+	}
+
+	# Allocate temporary working space.
+	call smark (sp)
+	call salloc (vote, axes[1], TY_INT)
+	call salloc (vindex, axes[1], TY_INT)
+	call salloc (lmatch, axes[1], TY_INT)
+	call salloc (rmatch, axes[2], TY_INT)
+	call amovki (NO, Memi[lmatch], axes[1])
+	call amovki (NO, Memi[rmatch], axes[2])
+
+	# Find the maximum value in the mask.
+	minvote = MAX_INT
+	maxvote = -MAX_INT
+	do i = 1, axes[2] {
+	    laxes[1] = 1
+	    laxes[2] = i
+	    if (! pl_linenotempty (pl, laxes))
+		next
+	    call pl_glpi (pl, laxes, Memi[vote], 16, axes[1], PIX_SRC)
+	    call alimi (Memi[vote], axes[1], tminvote, tmaxvote)
+	    minvote = min (minvote, tminvote)
+	    maxvote = max (maxvote, tmaxvote)
+	}
+	if (maxvote < 0) {
+	    maxvote = 0
+	    hmaxvote = 0
+	} else
+	    hmaxvote = maxvote / 2
+
+	# Vote on the matched pairs.
+	ninter = 0
+	if (maxvote > 0) {
+	    do j = 1, axes[2] {
+
+	        # Sort the vote array.
+	        do i = 1, axes[1]
+		    Memi[vindex+i-1] = i
+	        laxes[1] = 1
+	        laxes[2] = j
+	        call pl_glpi (pl, laxes, Memi[vote], 16, axes[1], PIX_SRC)
+	        call rg_qsorti (Memi[vote], Memi[vindex], Memi[vindex],
+		    axes[1])
+
+	        # Reject points which have no votes, which have only a
+		# single vote if the maximum number of votest is > 1,
+		# less or equal to half the  maximum number of votes,
+		# the same number of votes as the next largest index,
+		# or which have already been matched.
+
+		vp = Memi[vindex+axes[1]-1]
+		pvp = Memi[vindex+axes[1]-2]
+	        if (Memi[vote+vp-1] <= 0)
+		    next
+		if (Memi[vote+vp-1] == Memi[vote+pvp-1])
+		    next
+		if (Memi[vote+vp-1] <= hmaxvote)
+		    next
+		if (Memi[lmatch+vp-1] == YES || Memi[rmatch+j-1] == YES)
+		    next
+	        if (Memi[vote+vp-1] == 1 && (maxvote > 1 || nmatch > 1))
+		    next
+
+		ninter = ninter + 1
+		reftri[ninter, RG_MATCH] = j
+		listri[ninter,RG_MATCH] = vp
+		Memi[rmatch+j-1] = YES
+		Memi[lmatch+vp-1] = YES
+	    }
+	} else if (maxvote > 0) {
+	}
+
+	call sfree (sp)
+	call pl_close (pl)
+
+	return (ninter)
+end
+
+
+# RG_MLINCOEFF -- Compute the coefficients of a new linear transformation
+# using the first one to three matched stars as input.
+
+int procedure rg_mlincoeff (xref, yref, xlist, ylist, reftri, nmrtri, listri,
+	nmltri, nmatch, coeff, ncoeff)
+
+real	xref[ARB]		#I the x reference coordinates
+real	yref[ARB]		#I the y reference coordinates
+real	xlist[ARB]		#I the x list coordinates
+real	ylist[ARB]		#I the y list coordinates
+int	reftri[nmrtri,ARB]	#I list of reference triangles
+int	nmrtri			#I maximum number of reference triangles
+int	listri[nmltri,ARB]	#I list of reference triangles
+int	nmltri			#I maximum number of list triangles
+int	nmatch			#I number of matches
+real	coeff[ARB]		#O the new computed coefficients
+int	ncoeff			#I the number of coefficients
+
+int	i, rindex, lindex, stat
+pointer	sp, xr, yr, xin, yin
+int	rg_lincoeff()
+
+begin
+	if (nmatch <= 0)
+	    return (ERR)
+
+	call smark (sp)
+	call salloc (xr, nmatch, TY_REAL)
+	call salloc (yr, nmatch, TY_REAL)
+	call salloc (xin, nmatch, TY_REAL)
+	call salloc (yin, nmatch, TY_REAL)
+
+	# Load the points to be fit.
+	do i = 1, nmatch {
+	    rindex = reftri[i,RG_MATCH]
+	    lindex = listri[i,RG_MATCH]
+	    Memr[xr+i-1] = xref[rindex] 
+	    Memr[yr+i-1] = yref[rindex] 
+	    Memr[xin+i-1] = xlist[lindex]
+	    Memr[yin+i-1] = ylist[lindex]
+	}
+
+	# Compute the new coefficients.
+	stat = rg_lincoeff (Memr[xr], Memr[yr], Memr[xin], Memr[yin],
+	    nmatch, coeff, ncoeff)
+
+	call sfree (sp)
+
+	return (stat)
+end
+
+
+# RG_MWRITE -- Write out the intersection of the two matched pixel lists to the
+# output file.
+
+procedure rg_mwrite (ofd, xref, yref, rlineno, xlist, ylist, ilineno,
+	reftri, nmrtri, listri, nmltri, nmatch, xformat, yformat)
+
+int	ofd			#I the output file descriptor
+real	xref[ARB]		#I the x reference coordinates
+real	yref[ARB]		#I the y reference coordinates
+int	rlineno[ARB]		#I the reference coordinate line numbers
+real	xlist[ARB]		#I the x list coordinates
+real	ylist[ARB]		#I the y list coordinates
+int	ilineno[ARB]		#I the input list line numbers
+int	reftri[nmrtri,ARB]	#I list of reference triangles
+int	nmrtri			#I maximum number of reference triangles
+int	listri[nmltri,ARB]	#I list of reference triangles
+int	nmltri			#I maximum number of list triangles
+int	nmatch			#I number of matches
+char	xformat[ARB]		#I the output x column format
+char	yformat[ARB]		#I the output y column format
+
+int	i, lindex, rindex
+pointer	sp, fmtstr
+
+begin
+	call smark (sp)
+	call salloc (fmtstr, SZ_LINE, TY_CHAR)
+
+	# Construct the format string.
+	call sprintf (Memc[fmtstr], SZ_LINE, "%s %s  %s %s  %%5d %%5d\n")
+	if (xformat[1] == EOS)
+	    call pargstr ("%13.7g")
+	else
+	    call pargstr (xformat)
+	if (yformat[1] == EOS)
+	    call pargstr ("%13.7g")
+	else
+	    call pargstr (yformat)
+	if (xformat[1] == EOS)
+	    call pargstr ("%13.7g")
+	else
+	    call pargstr (xformat)
+	if (yformat[1] == EOS)
+	    call pargstr ("%13.7g")
+	else
+	    call pargstr (yformat)
+
+	do i = 1, nmatch {
+	    rindex = reftri[i,RG_MATCH]
+	    lindex = listri[i,RG_MATCH]
+	    call fprintf (ofd, Memc[fmtstr])
+		call pargr (xref[rindex])
+		call pargr (yref[rindex])
+		call pargr (xlist[lindex])
+		call pargr (ylist[lindex])
+		call pargi (rlineno[rindex])
+		call pargi (ilineno[lindex])
+	}
+
+	call sfree (sp)
+end
+
+
+# RG_LMWRITE -- Write out the intersection of the matched celestial coordinate
+# and pixel lists to the output file.
+
+procedure rg_lmwrite (ofd, lngref, latref, rlineno, xlist, ylist, ilineno,
+	reftri, nmrtri, listri, nmltri, nmatch, lngformat, latformat,
+	xformat, yformat)
+
+int	ofd			#I the output file descriptor
+double	lngref[ARB]		#I the x reference coordinates
+double	latref[ARB]		#I the y reference coordinates
+int	rlineno[ARB]		#I the reference coordinate line numbers
+real	xlist[ARB]		#I the x list coordinates
+real	ylist[ARB]		#I the y list coordinates
+int	ilineno[ARB]		#I the input list line numbers
+int	reftri[nmrtri,ARB]	#I list of reference triangles
+int	nmrtri			#I maximum number of reference triangles
+int	listri[nmltri,ARB]	#I list of reference triangles
+int	nmltri			#I maximum number of list triangles
+int	nmatch			#I number of matches
+char	lngformat[ARB]		#I the output longitude column format
+char	latformat[ARB]		#I the output latitude column format
+char	xformat[ARB]		#I the output x column format
+char	yformat[ARB]		#I the output y column format
+
+int	i, lindex, rindex
+pointer	sp, fmtstr
+
+begin
+	call smark (sp)
+	call salloc (fmtstr, SZ_LINE, TY_CHAR)
+
+	# Construct the format string.
+	call sprintf (Memc[fmtstr], SZ_LINE, "%s %s  %s %s  %%5d %%5d\n")
+	if (lngformat[1] == EOS)
+	    call pargstr ("%13.7g")
+	else
+	    call pargstr (lngformat)
+	if (latformat[1] == EOS)
+	    call pargstr ("%13.7g")
+	else
+	    call pargstr (latformat)
+	if (xformat[1] == EOS)
+	    call pargstr ("%13.7g")
+	else
+	    call pargstr (xformat)
+	if (yformat[1] == EOS)
+	    call pargstr ("%13.7g")
+	else
+	    call pargstr (yformat)
+
+	do i = 1, nmatch {
+	    rindex = reftri[i,RG_MATCH]
+	    lindex = listri[i,RG_MATCH]
+	    call fprintf (ofd, Memc[fmtstr])
+		call pargd (lngref[rindex])
+		call pargd (latref[rindex])
+		call pargr (xlist[lindex])
+		call pargr (ylist[lindex])
+		call pargi (rlineno[rindex])
+		call pargi (ilineno[lindex])
+	}
+
+	call sfree (sp)
+end
+
+
+# RG_FACTORIAL -- Compute the combinatorial function which is defined as
+# n! / ((n - ngroup)! * ngroup!).
+
+int procedure rg_factorial (n, ngroup)
+
+int	n		#I input argument
+int	ngroup		#I combinatorial factor
+
+int	i, fac, grfac
+
+begin
+	if (n <= 0)
+	    return (1)
+
+	fac = n
+	do i = n - 1, n - 3 + 1, -1
+	    fac = fac * i
+
+	grfac = ngroup
+	do i = ngroup - 1, 2, -1
+	    grfac = grfac * i
+
+	return (fac / grfac)
+end
+
+
+# RG_MODED -- Compute mode of an array.  The mode is found by binning
+# with a bin size based on the data range over a fraction of the
+# pixels about the median and a bin step which may be smaller than the
+# bin size.  If there are too few points the median is returned.
+# The input array must be sorted.
+
+double procedure rg_moded (a, npts, nmin, zrange, fzbin, fzstep)
+
+double  a[npts]                 #I the sorted input data array
+int     npts                    #I the number of points
+int     nmin                    #I the minimum number of points
+double  zrange                  #I fraction of pixels around median to use
+double  fzbin                   #I the bin size for the mode search
+double  fzstep                  #I the step size for the mode search
+
+int     x1, x2, x3, nmax
+double  zstep, zbin, y1, y2, mode
+bool    fp_equald()
+
+begin
+        # If there are too few points return the median.
+        if (npts < nmin) {
+            if (mod (npts,2) == 1)
+                return (a[1+npts/2])
+            else
+                return ((a[npts/2] + a[1+npts/2]) / 2.0d0)
+        }
+
+        # Compute the data range that will be used to do the mode search.
+        # If the data has no range then the constant value will be returned.
+        x1 = max (1, int (1.0d0 + npts * (1.0d0 - zrange) / 2.0d0))
+        x3 = min (npts, int (1.0d0 + npts * (1.0d0 + zrange) / 2.0d0))
+        if (fp_equald (a[x1], a[x3]))
+            return (a[x1])
+
+
+        # Compute the bin and step size. The bin size is based on the
+        # data range over a fraction of the pixels around the median
+        # and a bin step which may be smaller than the bin size.
+
+        zstep = fzstep #* (a[x3] - a[x1])
+        zbin = fzbin #* (a[x3] - a[x1])
+
+        nmax = 0
+        x2 = x1
+        for (y1 = a[x1]; x2 < x3; y1 = y1 + zstep) {
+            for (; a[x1] < y1; x1 = x1 + 1)
+                ;
+            y2 = y1 + zbin
+            for (; (x2 < x3) && (a[x2] < y2); x2 = x2 + 1)
+                ;
+            if (x2 - x1 > nmax) {
+                nmax = x2 - x1
+                if (mod (x2+x1,2) == 0)
+                    mode = a[(x2+x1)/2]
+                else
+                    mode = (a[(x2+x1)/2] + a[(x2+x1)/2+1]) / 2.0d0
+            }
+        }
+
+        return (mode)
+end
+
+
+#define  NMIN    10        # Minimum number of pixels for mode calculation
+#define  ZRANGE  0.8d0     # Fraction of pixels about median to use
+#define  ZSTEP   0.01d0    # Step size for search for mode
+#define  ZBIN    0.1d0     # Bin size for mode.
+#
+## RG_MODED -- Compute mode of an array.  The mode is found by binning
+## with a bin size based on the data range over a fraction of the
+## pixels about the median and a bin step which may be smaller than the
+## bin size.  If there are too few points the median is returned.
+## The input array must be sorted.
+#
+#double	procedure rg_moded (a, n)
+#
+#double  a[n]                    # Data array
+#int     n                       # Number of points
+#
+#int     i, j, k, nmax
+#real    z1, z2, zstep, zbin
+#double  mode
+#bool    fp_equald()
+#
+#begin
+#        if (n < NMIN)
+#            return (a[n/2])
+#
+#        # Compute the mode.  The array must be sorted.  Consider a
+#        # range of values about the median point.  Use a bin size which
+#        # is ZBIN of the range.  Step the bin limits in ZSTEP fraction of
+#        # the bin size.
+#
+#	i = 1 + n * (1. - ZRANGE) / 2.0d0
+#        j = 1 + n * (1. + ZRANGE) / 2.0d0
+#        z1 = a[i]
+#        z2 = a[j]
+#        if (fp_equald (z1, z2)) {
+#            mode = z1
+#            return (mode)
+#        }
+#
+#        zstep = ZSTEP * (z2 - z1)
+#        zbin = ZBIN * (z2 - z1)
+#
+#        z1 = z1 - zstep
+#        k = i
+#        nmax = 0
+#        repeat {
+#            z1 = z1 + zstep
+#            z2 = z1 + zbin
+#            for (; i < j && a[i] < z1; i=i+1)
+#                ;
+#            for (; k < j && a[k] < z2; k=k+1)
+#                ;
+#            if (k - i > nmax) {
+#                nmax = k - i
+#                mode = a[(i+k)/2]
+#            }
+#        } until (k >= j)
+#
+#        return (mode)
+#end
+#