include include include # 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) 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 #