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