PRO findmatch, x1, y1, x2, y2, mins, maxs, critical_unitarity, nmatch, subs1, subs2, lpar, ldet, QUIET = quiet
; Description: This module searches for matching points between two coordinate lists
; of two-dimensional points where each list corresponds to a different
; coordinate space. The coordinates of the points in the first coordinate
; space S1 are specified by "(x1,y1)", and the coordinates of the points
; in the second coordinate space S2 are specified by "(x2,y2)". The
; algorithm used for point matching is based on the similar triangles
; technique, and specifically the algorithm used here is based on the
; methodology and ideas in the paper by Pal & Bakos (Pal A. & Bakos G.A.,
; 2006, PASP, 118, 1474).
; To summarise, (extended) Delaunay triangulations of each coordinate
; list are generated and triangles are matched in "continuous" triangle
; space using symmetric point matching. Each triangle match for which the
; ratio of triangle perimeters does not satisfy the minimum and maximum
; transformation scale factors ("mins" and "maxs" respectively) is rejected.
; A vote algorithm is used to identify candidate point pair matches from
; which an initial linear transformation is derived. If the unitarity of
; the linear transformation satisfies the "critical_unitarity" parameter,
; then the module reports the number of point pair matches "nmatch" along
; with the lists of subscripts "subs1" and "subs2" of these points in each
; coordinate list. The module also returns the fitted initial linear
; transformation using the parameters "lpar" and "ldet".
; In the case that a particular (extended) Delaunay triangulation fails
; to produce a set of reliable point pair matches, then either the extension
; level of the Delaunay triangulation for each coordinate list is increased
; or the chirality of the triangle space coordinates for the second
; coordinate list is reversed, and the above algorithm is repeated. A
; maximum extension level of two for the extended Delaunay triangulation has
; been chosen, because our implementation of the generation of an extended
; Delaunay triangulation produces more unique triangles than the method
; implemented by Pal & Bakos 2006 for their maximum extension level of four
; (see "triangulate_extend.pro").
;
; Input Parameters:
;
; x1 - FLOAT/DOUBLE VECTOR - A one-dimensional vector of x coordinates for the points
; in the first coordinate list.
; y1 - FLOAT/DOUBLE VECTOR - A one-dimensional vector of y coordinates for the points
; in the first coordinate list.
; x2 - FLOAT/DOUBLE VECTOR - A one-dimensional vector of x coordinates for the points
; in the second coordinate list.
; y2 - FLOAT/DOUBLE VECTOR - A one-dimensional vector of y coordinates for the points
; in the second coordinate list.
; mins - FLOAT/DOUBLE - The minimum possible transformation scale factor (magnification)
; for the transformation mapping S1 to S2. If "mins" is negative,
; then this restriction to the transformation scale factor is
; ignored.
; maxs - FLOAT/DOUBLE - The maximum possible transformation scale factor (magnification)
; for the transformation mapping S1 to S2. If "maxs" is negative
; or zero, then this restriction to the transformation scale
; factor is ignored.
; critical_unitarity - FLOAT/DOUBLE - The maximum acceptable unitarity of the initial
; linear transformation fitted to the candidate point
; pair matches. This should be set to some factor (~4)
; times the distortion factor f_d due to the
; difference between the orthographic and gnomonic
; projections given by the approximation:
;
; f_d ~ 1 - cos(d)
;
; where d is the field-of-view (FOV) in degrees. For
; any FOV less than 1.25 degrees, then setting the
; value of "critical_unitarity" to 0.001 is recommended
; (and sufficient). A negative or zero value for this
; parameter will be reset to a value of "0.001".
;
; N.B: Both the first coordinate list and the second coordinate list must not contain
; duplicate entries, otherwise the triangle matching algorithm will fail.
;
; Output Parameters:
;
; nmatch - LONG - The number of matching point pairs between the coordinate lists.
; subs1 - LONG VECTOR - A one-dimensional vector of subscripts from the first
; coordinate list corresponding to the matched points. This
; vector has "nmatch" elements, and each element in "subs1"
; corresponds in order with the subscripts of the matching
; points from the second coordinate list "subs2".
; subs2 - LONG VECTOR - A one-dimensional vector of subscripts from the second
; coordinate list corresponding to the matched points. This
; vector has "nmatch" elements, and each element in "subs2"
; corresponds in order with the subscripts of the matching
; points from the first coordinate list "subs1".
; lpar - DOUBLE VECTOR - A vector of six DOUBLES representing the fitted parameters
; A, B, C, D, E and F, in this order, of the initial linear
; transformation mapping S1 to S2. The linear transformation
; is defined by:
;
; x2 = A*x1 + B*y1 + E
; y2 = C*x1 + D*y1 + F
;
; ldet - DOUBLE - The determinant of the fitted initial linear transformation defined
; by:
;
; ldet = (A*D) - (B*C)
;
; Keywords:
;
; If the keyword QUIET is set (as "/QUIET"), then the program will not print any
; progress messages on screen.
;
;
; Author: Dan Bramich (dan.bramich@hotmail.co.uk)
;
; History:
;
; 24/04/2009 - Further module optimisation performed (dmb).
; 30/07/2008 - Module created (dmb).
;
; This module does the following:
;
; 1) Checks that the parameters "x1" and "y1" are one-dimensional vectors of the correct
; number type.
;
; 2) Checks that the first coordinate list has at least three points and that the number
; of x coordinates is equal to the number of y coordinates for this list.
;
; 3) Checks that the parameters "x2" and "y2" are one-dimensional vectors of the correct
; number type.
;
; 4) Checks that the second coordinate list has at least three points and that the number
; of x coordinates is equal to the number of y coordinates for this list.
;
; 5) Sets up a vote array to contain the votes for each possible point pair combination.
;
; 6) While the extension level of the Delaunay triangulation is two or less, and while
; the number of successfully matched point pairs is zero, the program will repeat
; steps (7)-(13). The initial extension level of the Delaunay triangulation is zero.
;
; 7) Generates an (extended) Delaunay triangulation (see Pal & Bakos 2006 and
; "triangulate_extend.pro" from DanIDL) using the current extension level for each
; coordinate list.
;
; 8) Generates the "continuous" triangle space coordinates (see Pal & Bakos 2006 and
; "triangle_space.pro" from DanIDL) for each triangulation.
;
; 9) Matches the triangles between the two triangle lists using symmetric point matching
; (see Pal & Bakos 2006 and "sympoint_match_2d.pro" from DanIDL). Each triangle match
; for which the ratio of triangle perimeters does not satisfy the minimum and maximum
; transformation scale factors ("mins" and "maxs" respectively) is rejected.
;
; 10) If there is at least one matching triangle between the two triangle lists, then
; the module executes point voting using the list of matching triangles. Each matching
; triangle yields three votes, one for each pair of matching vertices between the
; two coordinate lists. This is in contrast to the Pal & Bakos 2006 algorithm which
; gives Nt votes to each of the three point pairs from the closest matching pair of
; triangles in triangle space coordinates, where Nt is the number of matching triangles,
; and gives Nt-1 votes to each of the three point pairs from the next closest matching
; triangle pair, and so on.
; A set of candidate point pair matches are chosen such that each candidate point
; pair match has at least one vote and also at least half of the number of votes of the
; highest vote getting point pair.
; If there are no matching triangles between the two triangle lists, then the number
; of candidate point pair matches is set to zero.
;
; 11) If there are at least three candidate point pair matches, then the program fits an
; initial linear transformation to the candidate point pair matches that transforms the
; coordinate space S1 to the coordinate space S2. If there are less than three candidate
; point pair matches, then the program moves on to step (13).
;
; 12) If the unitarity of the fitted linear transformation is less than the critical
; unitarity "critical_unitarity", then the program declares the matching process a
; success, and reports the number of matching point pairs "nmatch" along with the lists
; of subscripts "subs1" and "subs2" of these points in each coordinate list. The program
; also returns the fitted initial linear transformation using the parameters "lpar" and
; "ldet". Otherwise, the program moves on to step (13).
; Note that, in the case that the matching process is a success, "nmatch" will always
; be returned with a value of at least three, due to the necessity of fitting an initial
; linear transformation in step (11).
;
; 13) In reaching this step, the program declares the matching process for the current
; triangulations a failure. If the extension level of the Delaunay triangulation is
; less than two, and the x coordinates of the points in the second coordinate list have
; not yet been reversed, then the extension level of the Delaunay triangulation is
; incremented by one, and the program returns to step (7).
; If the extension level of the Delaunay triangulation is equal to two, and the x
; coordinates of the points in the second coordinate list have not yet been reversed,
; then the x coordinates of the points in the second coordinate list are reversed and
; the extension level of the Delaunay triangulation is reset to zero, and the program
; returns to step (7).
; If the extension level of the Delaunay triangulation is less than two, and the x
; coordinates of the points in the second coordinate list have already been reversed,
; then the extension level of the Delaunay triangulation is incremented by one, and
; the program returns to step (7).
; Finally, if the extension level of the Delaunay triangulation is equal to two, and
; the x coordinates of the points in the second coordinate list have already been reversed,
; then the program ends having failed to find a reliable set of point pair matches, and
; returns with "nmatch" set to zero.
;If the keyword QUIET is set, then do not print progress messages
;Check that "x1" is a vector of the correct number type
;Check that "y1" is a vector of the correct number type
;If there are less than 3 points in the first coordinate list, or the number of x coordinates
;is not equal to the number of y coordinates for this list, then return with zero matches
;Check that "x2" is a vector of the correct number type
;Check that "y2" is a vector of the correct number type
;If there are less than 3 points in the second coordinate list, or the number of x coordinates
;is not equal to the number of y coordinates for this list, then return with zero matches
;Check that "mins", "maxs" and "critical_unitarity" are all numbers
;Set up a vote array
;Attempt to find a reliable set of point pair candidate matches starting with the Delaunay
;triangulation. If this fails to produce a reliable initial linear transformation, then use
;successive levels of extended Delaunay triangulations. If this also fails, then switch the
;chirality of the triangle space coordinates for the second coordinate list and retry,
;starting with the Delaunay triangulation and continuing with the extended Delaunay
;triangulations.
;Generate the relevant triangulation for both the first and second coordinate lists
;Generate the triangle space coordinates for both triangulations
;Match the triangles between the two triangle lists
;Continue if there is at least one matching triangle
;Execute the point voting using the list of matching triangles in order to select candidate
;point pair matches from the two coordinate lists
;Sort the vote array
;Choose the set of candidate point pair matches with at least half of the votes of the
;highest vote getting point pair
;If there are no matching triangles
;If there are at least three candidate point pair matches, then derive an initial linear
;transformation
;Fit a linear transformation to the candidate point pair matches that transforms the
;coordinate space of the points from the first coordinate list to the coordinate space
;of the points from the second coordinate list
;If the fit of the initial linear transformation was successful
;Calculate the unitarity of the fitted linear transformation
;If the unitarity of the fitted linear transformation is less than the critical
;unitarity, then declare the matching process a success and save the details
;of the fitted initial linear transformation
;If the unitarity of the fitted linear transformation is equal to or greater than
;the critical unitarity
;If the fit of the initial linear transformation was unsuccessful
;If there are not enough candidate point pair matches to derive an initial linear transformation
;If the point matching has been unsuccessful for the current triangulations and triangle
;space chirality, then move on to the next attempt
;Finish