LSST Applications  21.0.0-172-gfb10e10a+18fedfabac,22.0.0+297cba6710,22.0.0+80564b0ff1,22.0.0+8d77f4f51a,22.0.0+a28f4c53b1,22.0.0+dcf3732eb2,22.0.1-1-g7d6de66+2a20fdde0d,22.0.1-1-g8e32f31+297cba6710,22.0.1-1-geca5380+7fa3b7d9b6,22.0.1-12-g44dc1dc+2a20fdde0d,22.0.1-15-g6a90155+515f58c32b,22.0.1-16-g9282f48+790f5f2caa,22.0.1-2-g92698f7+dcf3732eb2,22.0.1-2-ga9b0f51+7fa3b7d9b6,22.0.1-2-gd1925c9+bf4f0e694f,22.0.1-24-g1ad7a390+a9625a72a8,22.0.1-25-g5bf6245+3ad8ecd50b,22.0.1-25-gb120d7b+8b5510f75f,22.0.1-27-g97737f7+2a20fdde0d,22.0.1-32-gf62ce7b1+aa4237961e,22.0.1-4-g0b3f228+2a20fdde0d,22.0.1-4-g243d05b+871c1b8305,22.0.1-4-g3a563be+32dcf1063f,22.0.1-4-g44f2e3d+9e4ab0f4fa,22.0.1-42-gca6935d93+ba5e5ca3eb,22.0.1-5-g15c806e+85460ae5f3,22.0.1-5-g58711c4+611d128589,22.0.1-5-g75bb458+99c117b92f,22.0.1-6-g1c63a23+7fa3b7d9b6,22.0.1-6-g50866e6+84ff5a128b,22.0.1-6-g8d3140d+720564cf76,22.0.1-6-gd805d02+cc5644f571,22.0.1-8-ge5750ce+85460ae5f3,master-g6e05de7fdc+babf819c66,master-g99da0e417a+8d77f4f51a,w.2021.48
LSST Data Management Base Package
pessimistic_pattern_matcher_b_3D.py
Go to the documentation of this file.
1 
2 import numpy as np
3 from scipy.optimize import least_squares
4 from scipy.spatial import cKDTree
5 from scipy.stats import sigmaclip
6 
7 from .pessimisticPatternMatcherUtils import (find_candidate_reference_pair_range,
8  create_pattern_spokes,)
9 import lsst.pipe.base as pipeBase
10 
11 
12 def _rotation_matrix_chi_sq(flattened_rot_matrix,
13  pattern_a,
14  pattern_b,
15  max_dist_rad):
16  """Compute the squared differences for least squares fitting.
17 
18  Given a flattened rotation matrix, one N point pattern and another N point
19  pattern to transform into to, compute the squared differences between the
20  points in the two patterns after the rotation.
21 
22  Parameters
23  ----------
24  flattened_rot_matrix : `numpy.ndarray`, (9, )
25  A flattened array representing a 3x3 rotation matrix. The array is
26  flattened to comply with the API of scipy.optimize.least_squares.
27  Flattened elements are [[0, 0], [0, 1], [0, 2], [1, 0]...]
28  pattern_a : `numpy.ndarray`, (N, 3)
29  A array containing N, 3 vectors representing the objects we would like
30  to transform into the frame of pattern_b.
31  pattern_b : `numpy.ndarray`, (N, 3)
32  A array containing N, 3 vectors representing the reference frame we
33  would like to transform pattern_a into.
34  max_dist_rad : `float`
35  The maximum distance allowed from the pattern matching. This value is
36  used as the standard error for the resultant chi values.
37 
38  Returns
39  -------
40  noralized_diff : `numpy.ndarray`, (9,)
41  Array of differences between the vectors representing of the source
42  pattern rotated into the reference frame and the converse. This is
43  used to minimize in a least squares fitter.
44  """
45  # Unflatten the rotation matrix
46  rot_matrix = flattened_rot_matrix.reshape((3, 3))
47  # Compare the rotated source pattern to the references.
48  rot_pattern_a = np.dot(rot_matrix, pattern_a.transpose()).transpose()
49  diff_pattern_a_to_b = rot_pattern_a - pattern_b
50  # Return the flattened differences and length tolerances for use in a least
51  # squares fitter.
52  return diff_pattern_a_to_b.flatten() / max_dist_rad
53 
54 
56  """Class implementing a pessimistic version of Optimistic Pattern Matcher
57  B (OPMb) from Tabur 2007. See `DMTN-031 <http://ls.st/DMTN-031`_
58 
59  Parameters
60  ----------
61  reference_array : `numpy.ndarray`, (N, 3)
62  spherical points x, y, z of to use as reference objects for
63  pattern matching.
64  log : `lsst.log.Log` or `logging.Logger`
65  Logger for outputting debug info.
66 
67  Notes
68  -----
69  The class loads and stores the reference object
70  in a convenient data structure for matching any set of source objects that
71  are assumed to contain each other. The pessimistic nature of the algorithm
72  comes from requiring that it discovers at least two patterns that agree on
73  the correct shift and rotation for matching before exiting. The original
74  behavior of OPMb can be recovered simply. Patterns matched between the
75  input datasets are n-spoked pinwheels created from n+1 points. Refer to
76  DMTN #031 for more details. http://github.com/lsst-dm/dmtn-031
77  """
78 
79  def __init__(self, reference_array, log):
80  self._reference_array_reference_array = reference_array
81  self._n_reference_n_reference = len(self._reference_array_reference_array)
82  self.loglog = log
83 
84  if self._n_reference_n_reference > 0:
85  self._build_distances_and_angles_build_distances_and_angles()
86  else:
87  raise ValueError("No reference objects supplied")
88 
89  def _build_distances_and_angles(self):
90  """Create the data structures we will use to search for our pattern
91  match in.
92 
93  Throughout this function and the rest of the class we use id to
94  reference the position in the input reference catalog and index to
95  'index' into the arrays sorted on distance.
96  """
97  # Create empty arrays to store our pair information per
98  # reference object.
99  self._dist_array_dist_array = np.empty(
100  int(self._n_reference_n_reference * (self._n_reference_n_reference - 1) / 2),
101  dtype="float32")
102  self._id_array_id_array = np.empty(
103  (int(self._n_reference_n_reference * (self._n_reference_n_reference - 1) / 2), 2),
104  dtype="uint16")
105 
106  startIdx = 0
107  # Loop over reference objects storing pair distances and ids.
108  for ref_id, ref_obj in enumerate(self._reference_array_reference_array):
109  # Set the ending slicing index to the correct length for the
110  # pairs we are creating.
111  endIdx = startIdx + self._n_reference_n_reference - 1 - ref_id
112 
113  # Reserve and fill the ids of each reference object pair.
114  # 16 bit is safe for the id array as the catalog input from
115  # MatchPessimisticB is limited to a max length of 2 ** 16.
116  self._id_array_id_array[startIdx:endIdx, 0] = ref_id
117  self._id_array_id_array[startIdx:endIdx, 1] = np.arange(ref_id + 1,
118  self._n_reference_n_reference,
119  dtype="uint16")
120 
121  # Compute the vector deltas for each pair of reference objects.
122  # Compute and store the distances.
123  self._dist_array_dist_array[startIdx:endIdx] = np.sqrt(
124  ((self._reference_array_reference_array[ref_id + 1:, :]
125  - ref_obj) ** 2).sum(axis=1))
126  # Set startIdx of the slice to the end of the previous slice.
127  startIdx = endIdx
128 
129  # Sort each array on the pair distances for the initial
130  # optimistic pattern matcher lookup.
131  sorted_dist_args = self._dist_array_dist_array.argsort()
132  self._dist_array_dist_array = self._dist_array_dist_array[sorted_dist_args]
133  self._id_array_id_array = self._id_array_id_array[sorted_dist_args]
134 
135  def match(self, source_array, n_check, n_match, n_agree,
136  max_n_patterns, max_shift, max_rotation, max_dist,
137  min_matches, pattern_skip_array=None):
138  """Match a given source catalog into the loaded reference catalog.
139 
140  Given array of points on the unit sphere and tolerances, we
141  attempt to match a pinwheel like pattern between these input sources
142  and the reference objects this class was created with. This pattern
143  informs of the shift and rotation needed to align the input source
144  objects into the frame of the references.
145 
146  Parameters
147  ----------
148  source_array : `numpy.ndarray`, (N, 3)
149  An array of spherical x,y,z coordinates and a magnitude in units
150  of objects having a lower value for sorting. The array should be
151  of shape (N, 4).
152  n_check : `int`
153  Number of sources to create a pattern from. Not all objects may be
154  checked if n_match criteria is before looping through all n_check
155  objects.
156  n_match : `int`
157  Number of objects to use in constructing a pattern to match.
158  n_agree : `int`
159  Number of found patterns that must agree on their shift and
160  rotation before exiting. Set this value to 1 to recover the
161  expected behavior of Optimistic Pattern Matcher B.
162  max_n_patters : `int`
163  Number of patterns to create from the input source objects to
164  attempt to match into the reference objects.
165  max_shift : `float`
166  Maximum allowed shift to match patterns in arcseconds.
167  max_rotation : `float`
168  Maximum allowed rotation between patterns in degrees.
169  max_dist : `float`
170  Maximum distance in arcseconds allowed between candidate spokes in
171  the source and reference objects. Also sets that maximum distance
172  in the intermediate verify, pattern shift/rotation agreement, and
173  final verify steps.
174  pattern_skip_array : `int`
175  Patterns we would like to skip. This could be due to the pattern
176  being matched on a previous iteration that we now consider invalid.
177  This assumes the ordering of the source objects is the same
178  between different runs of the matcher which, assuming no object
179  has been inserted or the magnitudes have changed, it should be.
180 
181  Returns
182  -------
183  output_struct : `lsst.pipe.base.Struct`
184  Result struct with components
185 
186  - ``matches`` : (N, 2) array of matched ids for pairs. Empty list if no
187  match found (`numpy.ndarray`, (N, 2) or `list`)
188  - ``distances_rad`` : Radian distances between the matched objects.
189  Empty list if no match found (`numpy.ndarray`, (N,))
190  - ``pattern_idx``: Index of matched pattern. None if no match found
191  (`int`).
192  - ``shift`` : Magnitude for the shift between the source and reference
193  objects in arcseconds. None if no match found (`float`).
194  """
195 
196  # Given our input source_array we sort on magnitude.
197  sorted_source_array = source_array[source_array[:, -1].argsort(), :3]
198  n_source = len(sorted_source_array)
199 
200  # Initialize output struct.
201  output_match_struct = pipeBase.Struct(
202  match_ids=[],
203  distances_rad=[],
204  pattern_idx=None,
205  shift=None,
206  max_dist_rad=None,)
207 
208  if n_source <= 0:
209  self.loglog.warning("Source object array is empty. Unable to match. Exiting matcher.")
210  return None
211 
212  # To test if the shifts and rotations we find agree with each other,
213  # we first create two test points situated at the top and bottom of
214  # where the z axis on the sphere bisects the source catalog.
215  test_vectors = self._compute_test_vectors_compute_test_vectors(source_array[:, :3])
216 
217  # We now create an empty list of our resultant rotated vectors to
218  # compare the different rotations we find.
219  rot_vect_list = []
220 
221  # Convert the tolerances to values we will use in the code.
222  max_cos_shift = np.cos(np.radians(max_shift / 3600.))
223  max_cos_rot_sq = np.cos(np.radians(max_rotation)) ** 2
224  max_dist_rad = np.radians(max_dist / 3600.)
225 
226  # Loop through the sources from brightest to faintest, grabbing a
227  # chunk of n_check each time.
228  for pattern_idx in range(np.min((max_n_patterns,
229  n_source - n_match))):
230 
231  # If this pattern is one that we matched on the past but we
232  # now want to skip, we do so here.
233  if pattern_skip_array is not None and \
234  np.any(pattern_skip_array == pattern_idx):
235  self.loglog.debug(
236  "Skipping previously matched bad pattern %i...",
237  pattern_idx)
238  continue
239  # Grab the sources to attempt to create this pattern.
240  pattern = sorted_source_array[
241  pattern_idx: np.min((pattern_idx + n_check, n_source)), :3]
242 
243  # Construct a pattern given the number of points defining the
244  # pattern complexity. This is the start of the primary tests to
245  # match our source pattern into the reference objects.
246  construct_return_struct = \
247  self._construct_pattern_and_shift_rot_matrix_construct_pattern_and_shift_rot_matrix(
248  pattern, n_match, max_cos_shift, max_cos_rot_sq,
249  max_dist_rad)
250 
251  # Our struct is None if we could not match the pattern.
252  if construct_return_struct.ref_candidates is None or \
253  construct_return_struct.shift_rot_matrix is None or \
254  construct_return_struct.cos_shift is None or \
255  construct_return_struct.sin_rot is None:
256  continue
257 
258  # Grab the output data from the Struct object.
259  ref_candidates = construct_return_struct.ref_candidates
260  shift_rot_matrix = construct_return_struct.shift_rot_matrix
261  cos_shift = construct_return_struct.cos_shift
262  sin_rot = construct_return_struct.sin_rot
263 
264  # If we didn't match enough candidates we continue to the next
265  # pattern.
266  if len(ref_candidates) < n_match:
267  continue
268 
269  # Now that we know our pattern and shift/rotation are valid we
270  # store the the rotated versions of our test points for later
271  # use.
272  tmp_rot_vect_list = []
273  for test_vect in test_vectors:
274  tmp_rot_vect_list.append(np.dot(shift_rot_matrix, test_vect))
275  # Since our test point are in the source frame, we can test if
276  # their lengths are mostly preserved under the transform.
277  if not self._test_pattern_lengths_test_pattern_lengths(np.array(tmp_rot_vect_list),
278  max_dist_rad):
279  continue
280 
281  tmp_rot_vect_list.append(pattern_idx)
282  rot_vect_list.append(tmp_rot_vect_list)
283 
284  # Test if we have enough rotations, which agree, or if we
285  # are in optimistic mode.
286  if self._test_rotation_agreement_test_rotation_agreement(rot_vect_list, max_dist_rad) < \
287  n_agree - 1:
288  continue
289 
290  # Run the final verify step.
291  match_struct = self._final_verify_final_verify(source_array[:, :3],
292  shift_rot_matrix,
293  max_dist_rad,
294  min_matches)
295  if match_struct.match_ids is None or \
296  match_struct.distances_rad is None or \
297  match_struct.max_dist_rad is None:
298  continue
299 
300  # Convert the observed shift to arcseconds
301  shift = np.degrees(np.arccos(cos_shift)) * 3600.
302  # Print information to the logger.
303  self.loglog.debug("Succeeded after %i patterns.", pattern_idx)
304  self.loglog.debug("\tShift %.4f arcsec", shift)
305  self.loglog.debug("\tRotation: %.4f deg",
306  np.degrees(np.arcsin(sin_rot)))
307 
308  # Fill the struct and return.
309  output_match_struct.match_ids = \
310  match_struct.match_ids
311  output_match_struct.distances_rad = \
312  match_struct.distances_rad
313  output_match_struct.pattern_idx = pattern_idx
314  output_match_struct.shift = shift
315  output_match_struct.max_dist_rad = match_struct.max_dist_rad
316  return output_match_struct
317 
318  self.loglog.debug("Failed after %i patterns.", pattern_idx + 1)
319  return output_match_struct
320 
321  def _compute_test_vectors(self, source_array):
322  """Compute spherical 3 vectors at the edges of the x, y, z extent
323  of the input source catalog.
324 
325  Parameters
326  ----------
327  source_array : `numpy.ndarray`, (N, 3)
328  array of 3 vectors representing positions on the unit
329  sphere.
330 
331  Returns
332  -------
333  test_vectors : `numpy.ndarray`, (N, 3)
334  Array of vectors representing the maximum extents in x, y, z
335  of the input source array. These are used with the rotations
336  the code finds to test for agreement from different patterns
337  when the code is running in pessimistic mode.
338  """
339 
340  # Get the center of source_array.
341  if np.any(np.logical_not(np.isfinite(source_array))):
342  self.loglog.warning("Input source objects contain non-finite values. "
343  "This could end badly.")
344  center_vect = np.nanmean(source_array, axis=0)
345 
346  # So that our rotation test works over the full sky we compute
347  # the max extent in each Cartesian direction x,y,z.
348  xbtm_vect = np.array([np.min(source_array[:, 0]), center_vect[1],
349  center_vect[2]], dtype=np.float64)
350  xtop_vect = np.array([np.max(source_array[:, 0]), center_vect[1],
351  center_vect[2]], dtype=np.float64)
352  xbtm_vect /= np.sqrt(np.dot(xbtm_vect, xbtm_vect))
353  xtop_vect /= np.sqrt(np.dot(xtop_vect, xtop_vect))
354 
355  ybtm_vect = np.array([center_vect[0], np.min(source_array[:, 1]),
356  center_vect[2]], dtype=np.float64)
357  ytop_vect = np.array([center_vect[0], np.max(source_array[:, 1]),
358  center_vect[2]], dtype=np.float64)
359  ybtm_vect /= np.sqrt(np.dot(ybtm_vect, ybtm_vect))
360  ytop_vect /= np.sqrt(np.dot(ytop_vect, ytop_vect))
361 
362  zbtm_vect = np.array([center_vect[0], center_vect[1],
363  np.min(source_array[:, 2])], dtype=np.float64)
364  ztop_vect = np.array([center_vect[0], center_vect[1],
365  np.max(source_array[:, 2])], dtype=np.float64)
366  zbtm_vect /= np.sqrt(np.dot(zbtm_vect, zbtm_vect))
367  ztop_vect /= np.sqrt(np.dot(ztop_vect, ztop_vect))
368 
369  # Return our list of vectors for later rotation testing.
370  return np.array([xbtm_vect, xtop_vect, ybtm_vect, ytop_vect,
371  zbtm_vect, ztop_vect])
372 
373  def _construct_pattern_and_shift_rot_matrix(self, src_pattern_array,
374  n_match, max_cos_theta_shift,
375  max_cos_rot_sq, max_dist_rad):
376  """Test an input source pattern against the reference catalog.
377 
378  Returns the candidate matched patterns and their
379  implied rotation matrices or None.
380 
381  Parameters
382  ----------
383  src_pattern_array : `numpy.ndarray`, (N, 3)
384  Sub selection of source 3 vectors to create a pattern from
385  n_match : `int`
386  Number of points to attempt to create a pattern from. Must be
387  >= len(src_pattern_array)
388  max_cos_theta_shift : `float`
389  Maximum shift allowed between two patterns' centers.
390  max_cos_rot_sq : `float`
391  Maximum rotation between two patterns that have been shifted
392  to have their centers on top of each other.
393  max_dist_rad : `float`
394  Maximum delta distance allowed between the source and reference
395  pair distances to consider the reference pair a candidate for
396  the source pair. Also sets the tolerance between the opening
397  angles of the spokes when compared to the reference.
398 
399  Return
400  -------
401  output_matched_pattern : `lsst.pipe.base.Struct`
402  Result struct with components:
403 
404  - ``ref_candidates`` : ids of the matched pattern in the internal
405  reference_array object (`list` of `int`).
406  - ``src_candidates`` : Pattern ids of the sources matched
407  (`list` of `int`).
408  - ``shift_rot_matrix_src_to_ref`` : 3x3 matrix specifying the full
409  shift and rotation between the reference and source objects.
410  Rotates source into reference frame. `None` if match is not
411  found. (`numpy.ndarray`, (3, 3))
412  - ``shift_rot_matrix_ref_to_src`` : 3x3 matrix specifying the full
413  shift and rotation of the reference and source objects. Rotates
414  reference into source frame. None if match is not found
415  (`numpy.ndarray`, (3, 3)).
416  - ``cos_shift`` : Magnitude of the shift found between the two
417  patten centers. `None` if match is not found (`float`).
418  - ``sin_rot`` : float value of the rotation to align the already
419  shifted source pattern to the reference pattern. `None` if no match
420  found (`float`).
421  """
422 
423  # Create our place holder variables for the matched sources and
424  # references. The source list starts with the 0th and first indexed
425  # objects as we are guaranteed to use those and these define both
426  # the shift and rotation of the final pattern.
427  output_matched_pattern = pipeBase.Struct(
428  ref_candidates=[],
429  src_candidates=[],
430  shift_rot_matrix=None,
431  cos_shift=None,
432  sin_rot=None)
433 
434  # Create the delta vectors and distances we will need to assemble the
435  # spokes of the pattern.
436  src_delta_array = np.empty((len(src_pattern_array) - 1, 3))
437  src_delta_array[:, 0] = (src_pattern_array[1:, 0]
438  - src_pattern_array[0, 0])
439  src_delta_array[:, 1] = (src_pattern_array[1:, 1]
440  - src_pattern_array[0, 1])
441  src_delta_array[:, 2] = (src_pattern_array[1:, 2]
442  - src_pattern_array[0, 2])
443  src_dist_array = np.sqrt(src_delta_array[:, 0]**2
444  + src_delta_array[:, 1]**2
445  + src_delta_array[:, 2]**2)
446 
447  # Our first test. We search the reference dataset for pairs
448  # that have the same length as our first source pairs to with
449  # plus/minus the max_dist tolerance.
450 
451  # TODO: DM-32299 Convert all code below and subsiquent functions to
452  # C++.
453  candidate_range = find_candidate_reference_pair_range(
454  src_dist_array[0], self._dist_array_dist_array, max_dist_rad)
455 
456  def generate_ref_dist_indexes(low, high):
457  """Generator to loop outward from the midpoint between two values.
458 
459  Parameters
460  ----------
461  low : `int`
462  Minimum of index range.
463  high : `int`
464  Maximum of index range.
465 
466  Yields
467  ------
468  index : `int`
469  Current index.
470  """
471  mid = (high + low) // 2
472  for idx in range(high - low):
473  if idx%2 == 0:
474  mid += idx
475  yield mid
476  else:
477  mid -= idx
478  yield mid
479 
480  # Start our loop over the candidate reference objects.
481  for ref_dist_idx in generate_ref_dist_indexes(candidate_range[0],
482  candidate_range[1]):
483  # We have two candidates for which reference object corresponds
484  # with the source at the center of our pattern. As such we loop
485  # over and test both possibilities.
486  tmp_ref_pair_list = self._id_array_id_array[ref_dist_idx]
487  for pair_idx, ref_id in enumerate(tmp_ref_pair_list):
488  src_candidates = [0, 1]
489  ref_candidates = []
490  shift_rot_matrix = None
491  cos_shift = None
492  sin_rot = None
493  # Test the angle between our candidate ref center and the
494  # source center of our pattern. This angular distance also
495  # defines the shift we will later use.
496  ref_center = self._reference_array_reference_array[ref_id]
497  cos_shift = np.dot(src_pattern_array[0], ref_center)
498  if cos_shift < max_cos_theta_shift:
499  continue
500 
501  # We can now append this one as a candidate.
502  ref_candidates.append(ref_id)
503  # Test to see which reference object to use in the pair.
504  if pair_idx == 0:
505  ref_candidates.append(
506  tmp_ref_pair_list[1])
507  ref_delta = (self._reference_array_reference_array[tmp_ref_pair_list[1]]
508  - ref_center)
509  else:
510  ref_candidates.append(
511  tmp_ref_pair_list[0])
512  ref_delta = (self._reference_array_reference_array[tmp_ref_pair_list[0]]
513  - ref_center)
514 
515  # For dense fields it will be faster to compute the absolute
516  # rotation this pair suggests first rather than saving it
517  # after all the spokes are found. We then compute the cos^2
518  # of the rotation and first part of the rotation matrix from
519  # source to reference frame.
520  test_rot_struct = self._test_rotation_test_rotation(
521  src_pattern_array[0], ref_center, src_delta_array[0],
522  ref_delta, cos_shift, max_cos_rot_sq)
523  if test_rot_struct.cos_rot_sq is None or \
524  test_rot_struct.proj_ref_ctr_delta is None or \
525  test_rot_struct.shift_matrix is None:
526  continue
527 
528  # Get the data from the return struct.
529  cos_rot_sq = test_rot_struct.cos_rot_sq
530  proj_ref_ctr_delta = test_rot_struct.proj_ref_ctr_delta
531  shift_matrix = test_rot_struct.shift_matrix
532 
533  # Now that we have a candidate first spoke and reference
534  # pattern center, we mask our future search to only those
535  # pairs that contain our candidate reference center.
536  tmp_ref_id_array = np.arange(len(self._reference_array_reference_array),
537  dtype="uint16")
538  tmp_ref_dist_array = np.sqrt(
539  ((self._reference_array_reference_array
540  - self._reference_array_reference_array[ref_id])
541  ** 2).sum(axis=1)).astype("float32")
542  tmp_sorted_args = np.argsort(tmp_ref_dist_array)
543  tmp_ref_id_array = tmp_ref_id_array[tmp_sorted_args]
544  tmp_ref_dist_array = tmp_ref_dist_array[tmp_sorted_args]
545 
546  # Now we feed this sub data to match the spokes of
547  # our pattern.
548  pattern_spokes = create_pattern_spokes(
549  src_pattern_array[0], src_delta_array, src_dist_array,
550  self._reference_array_reference_array[ref_id], proj_ref_ctr_delta,
551  tmp_ref_dist_array, tmp_ref_id_array, self._reference_array_reference_array,
552  max_dist_rad, n_match)
553 
554  # If we don't find enough candidates we can continue to the
555  # next reference center pair.
556  if len(pattern_spokes) < n_match - 2:
557  continue
558 
559  # If we have the right number of matched ids we store these.
560  ref_candidates.extend([cand[0] for cand in pattern_spokes])
561  src_candidates.extend([cand[1] for cand in pattern_spokes])
562 
563  # We can now create our full rotation matrix for both the
564  # shift and rotation. Reminder shift, aligns the pattern
565  # centers, rotation rotates the spokes on top of each other.
566  shift_rot_struct = self._create_shift_rot_matrix_create_shift_rot_matrix(
567  cos_rot_sq, shift_matrix, src_delta_array[0],
568  self._reference_array_reference_array[ref_id], ref_delta)
569  # If we fail to create the rotation matrix, continue to the
570  # next objects.
571  if shift_rot_struct.sin_rot is None or \
572  shift_rot_struct.shift_rot_matrix is None:
573  continue
574 
575  # Get the data from the return struct.
576  sin_rot = shift_rot_struct.sin_rot
577  shift_rot_matrix = shift_rot_struct.shift_rot_matrix
578 
579  # Now that we have enough candidates we test to see if it
580  # passes intermediate verify. This shifts and rotates the
581  # source pattern into the reference frame and then tests that
582  # each source/reference object pair is within max_dist. It also
583  # tests the opposite rotation that is reference to source
584  # frame.
585  fit_shift_rot_matrix = self._intermediate_verify_intermediate_verify(
586  src_pattern_array[src_candidates],
587  self._reference_array_reference_array[ref_candidates],
588  shift_rot_matrix, max_dist_rad)
589 
590  if fit_shift_rot_matrix is not None:
591  # Fill the struct and return.
592  output_matched_pattern.ref_candidates = ref_candidates
593  output_matched_pattern.src_candidates = src_candidates
594  output_matched_pattern.shift_rot_matrix = \
595  fit_shift_rot_matrix
596  output_matched_pattern.cos_shift = cos_shift
597  output_matched_pattern.sin_rot = sin_rot
598  return output_matched_pattern
599 
600  return output_matched_pattern
601 
602  def _test_rotation(self, src_center, ref_center, src_delta, ref_delta,
603  cos_shift, max_cos_rot_sq):
604  """ Test if the rotation implied between the source
605  pattern and reference pattern is within tolerance. To test this
606  we need to create the first part of our spherical rotation matrix
607  which we also return for use later.
608 
609  Parameters
610  ----------
611  src_center : `numpy.ndarray`, (N, 3)
612  pattern.
613  ref_center : `numpy.ndarray`, (N, 3)
614  3 vector defining the center of the candidate reference pinwheel
615  pattern.
616  src_delta : `numpy.ndarray`, (N, 3)
617  3 vector delta between the source pattern center and the end of
618  the pinwheel spoke.
619  ref_delta : `numpy.ndarray`, (N, 3)
620  3 vector delta of the candidate matched reference pair
621  cos_shift : `float`
622  Cosine of the angle between the source and reference candidate
623  centers.
624  max_cos_rot_sq : `float`
625  candidate reference pair after shifting the centers on top of each
626  other. The function will return None if the rotation implied is
627  greater than max_cos_rot_sq.
628 
629  Returns
630  -------
631  result : `lsst.pipe.base.Struct`
632  Result struct with components:
633 
634  - ``cos_rot_sq`` : magnitude of the rotation needed to align the
635  two patterns after their centers are shifted on top of each
636  other. `None` if rotation test fails (`float`).
637  - ``shift_matrix`` : 3x3 rotation matrix describing the shift needed to
638  align the source and candidate reference center. `None` if rotation
639  test fails (`numpy.ndarray`, (N, 3)).
640  """
641 
642  # Make sure the sine is a real number.
643  if cos_shift > 1.0:
644  cos_shift = 1.
645  elif cos_shift < -1.0:
646  cos_shift = -1.
647  sin_shift = np.sqrt(1 - cos_shift ** 2)
648 
649  # If the sine of our shift is zero we only need to use the identity
650  # matrix for the shift. Else we construct the rotation matrix for
651  # shift.
652  if sin_shift > 0:
653  rot_axis = np.cross(src_center, ref_center)
654  rot_axis /= sin_shift
655  shift_matrix = self._create_spherical_rotation_matrix_create_spherical_rotation_matrix(
656  rot_axis, cos_shift, sin_shift)
657  else:
658  shift_matrix = np.identity(3)
659 
660  # Now that we have our shift we apply it to the src delta vector
661  # and check the rotation.
662  rot_src_delta = np.dot(shift_matrix, src_delta)
663  proj_src_delta = (rot_src_delta
664  - np.dot(rot_src_delta, ref_center) * ref_center)
665  proj_ref_delta = (ref_delta
666  - np.dot(ref_delta, ref_center) * ref_center)
667  cos_rot_sq = (np.dot(proj_src_delta, proj_ref_delta) ** 2
668  / (np.dot(proj_src_delta, proj_src_delta)
669  * np.dot(proj_ref_delta, proj_ref_delta)))
670  # If the rotation isn't in tolerance return None.
671  if cos_rot_sq < max_cos_rot_sq:
672  return pipeBase.Struct(
673  cos_rot_sq=None,
674  proj_ref_ctr_delta=None,
675  shift_matrix=None,)
676  # Return the rotation angle, the plane projected reference vector,
677  # and the first half of the full shift and rotation matrix.
678  return pipeBase.Struct(
679  cos_rot_sq=cos_rot_sq,
680  proj_ref_ctr_delta=proj_ref_delta,
681  shift_matrix=shift_matrix,)
682 
683  def _create_spherical_rotation_matrix(self, rot_axis, cos_rotation,
684  sin_rotion):
685  """Construct a generalized 3D rotation matrix about a given
686  axis.
687 
688  Parameters
689  ----------
690  rot_axis : `numpy.ndarray`, (3,)
691  3 vector defining the axis to rotate about.
692  cos_rotation : `float`
693  cosine of the rotation angle.
694  sin_rotation : `float`
695  sine of the rotation angle.
696 
697  Return
698  ------
699  shift_matrix : `numpy.ndarray`, (3, 3)
700  3x3 spherical, rotation matrix.
701  """
702 
703  rot_cross_matrix = np.array(
704  [[0., -rot_axis[2], rot_axis[1]],
705  [rot_axis[2], 0., -rot_axis[0]],
706  [-rot_axis[1], rot_axis[0], 0.]], dtype=np.float64)
707  shift_matrix = (cos_rotation*np.identity(3)
708  + sin_rotion*rot_cross_matrix
709  + (1. - cos_rotation)*np.outer(rot_axis, rot_axis))
710 
711  return shift_matrix
712 
713  def _create_shift_rot_matrix(self, cos_rot_sq, shift_matrix, src_delta,
714  ref_ctr, ref_delta):
715  """ Create the final part of our spherical rotation matrix.
716 
717  Parameters
718  ----------
719  cos_rot_sq : `float`
720  cosine of the rotation needed to align our source and reference
721  candidate patterns.
722  shift_matrix : `numpy.ndarray`, (3, 3)
723  3x3 rotation matrix for shifting the source pattern center on top
724  of the candidate reference pattern center.
725  src_delta : `numpy.ndarray`, (3,)
726  3 vector delta of representing the first spoke of the source
727  pattern
728  ref_ctr : `numpy.ndarray`, (3,)
729  3 vector on the unit-sphere representing the center of our
730  reference pattern.
731  ref_delta : `numpy.ndarray`, (3,)
732  3 vector delta made by the first pair of the reference pattern.
733 
734  Returns
735  -------
736  result : `lsst.pipe.base.Struct`
737  Result struct with components:
738 
739  - ``sin_rot`` : float sine of the amount of rotation between the
740  source and reference pattern. We use sine here as it is
741  signed and tells us the chirality of the rotation (`float`).
742  - ``shift_rot_matrix`` : float array representing the 3x3 rotation
743  matrix that takes the source pattern and shifts and rotates
744  it to align with the reference pattern (`numpy.ndarray`, (3,3)).
745  """
746  cos_rot = np.sqrt(cos_rot_sq)
747  rot_src_delta = np.dot(shift_matrix, src_delta)
748  delta_dot_cross = np.dot(np.cross(rot_src_delta, ref_delta), ref_ctr)
749 
750  sin_rot = np.sign(delta_dot_cross) * np.sqrt(1 - cos_rot_sq)
751  rot_matrix = self._create_spherical_rotation_matrix_create_spherical_rotation_matrix(
752  ref_ctr, cos_rot, sin_rot)
753 
754  shift_rot_matrix = np.dot(rot_matrix, shift_matrix)
755 
756  return pipeBase.Struct(
757  sin_rot=sin_rot,
758  shift_rot_matrix=shift_rot_matrix,)
759 
760  def _intermediate_verify(self, src_pattern, ref_pattern, shift_rot_matrix,
761  max_dist_rad):
762  """ Perform an intermediate verify step.
763 
764  Rotate the matches references into the source frame and test their
765  distances against tolerance. Only return true if all points are within
766  tolerance.
767 
768  Parameters
769  ----------
770  src_pattern : `numpy.ndarray`, (N,3)
771  Array of 3 vectors representing the points that make up our source
772  pinwheel pattern.
773  ref_pattern : `numpy.ndarray`, (N,3)
774  Array of 3 vectors representing our candidate reference pinwheel
775  pattern.
776  shift_rot_matrix : `numpy.ndarray`, (3,3)
777  3x3 rotation matrix that takes the source objects and rotates them
778  onto the frame of the reference objects
779  max_dist_rad : `float`
780  Maximum distance allowed to consider two objects the same.
781 
782  Returns
783  -------
784  fit_shift_rot_matrix : `numpy.ndarray`, (3,3)
785  Return the fitted shift/rotation matrix if all of the points in our
786  source pattern are within max_dist_rad of their matched reference
787  objects. Returns None if this criteria is not satisfied.
788  """
789  if len(src_pattern) != len(ref_pattern):
790  raise ValueError(
791  "Source pattern length does not match ref pattern.\n"
792  "\t source pattern len=%i, reference pattern len=%i" %
793  (len(src_pattern), len(ref_pattern)))
794 
795  if self._intermediate_verify_comparison_intermediate_verify_comparison(
796  src_pattern, ref_pattern, shift_rot_matrix, max_dist_rad):
797  # Now that we know our initial shift and rot matrix is valid we
798  # want to fit the implied transform using all points from
799  # our pattern. This is a more robust rotation matrix as our
800  # initial matrix only used the first 2 points from the source
801  # pattern to estimate the shift and rotation. The matrix we fit
802  # are allowed to be non unitary but need to preserve the length of
803  # the rotated vectors to within the match tolerance.
804  fit_shift_rot_matrix = least_squares(
805  _rotation_matrix_chi_sq,
806  x0=shift_rot_matrix.flatten(),
807  args=(src_pattern, ref_pattern, max_dist_rad)
808  ).x.reshape((3, 3))
809  # Do another verify in case the fits have wondered off.
810  if self._intermediate_verify_comparison_intermediate_verify_comparison(
811  src_pattern, ref_pattern, fit_shift_rot_matrix,
812  max_dist_rad):
813  return fit_shift_rot_matrix
814 
815  return None
816 
817  def _intermediate_verify_comparison(self, pattern_a, pattern_b,
818  shift_rot_matrix, max_dist_rad):
819  """Test the input rotation matrix against one input pattern and
820  a second one.
821 
822  If every point in the pattern after rotation is within a distance of
823  max_dist_rad to its candidate point in the other pattern, we return
824  True.
825 
826  Parameters
827  ----------
828  pattern_a : `numpy.ndarray`, (N,3)
829  Array of 3 vectors representing the points that make up our source
830  pinwheel pattern.
831  pattern_b : `numpy.ndarray`, (N,3)
832  Array of 3 vectors representing our candidate reference pinwheel
833  pattern.
834  shift_rot_matrix : `numpy.ndarray`, (3,3)
835  3x3 rotation matrix that takes the source objects and rotates them
836  onto the frame of the reference objects
837  max_dist_rad : `float`
838  Maximum distance allowed to consider two objects the same.
839 
840 
841  Returns
842  -------
843  bool
844  True if all rotated source points are within max_dist_rad of
845  the candidate references matches.
846  """
847  shifted_pattern_a = np.dot(shift_rot_matrix,
848  pattern_a.transpose()).transpose()
849  tmp_delta_array = shifted_pattern_a - pattern_b
850  tmp_dist_array = (tmp_delta_array[:, 0] ** 2
851  + tmp_delta_array[:, 1] ** 2
852  + tmp_delta_array[:, 2] ** 2)
853  return np.all(tmp_dist_array < max_dist_rad ** 2)
854 
855  def _test_pattern_lengths(self, test_pattern, max_dist_rad):
856  """ Test that the all vectors in a pattern are unit length within
857  tolerance.
858 
859  This is useful for assuring the non unitary transforms do not contain
860  too much distortion.
861 
862  Parameters
863  ----------
864  test_pattern : `numpy.ndarray`, (N, 3)
865  Test vectors at the maximum and minimum x, y, z extents.
866  max_dist_rad : `float`
867  maximum distance in radians to consider two points "agreeing" on
868  a rotation
869 
870  Returns
871  -------
872  test : `bool`
873  Length tests pass.
874  """
875  dists = (test_pattern[:, 0] ** 2
876  + test_pattern[:, 1] ** 2
877  + test_pattern[:, 2] ** 2)
878  return np.all(
879  np.logical_and((1 - max_dist_rad) ** 2 < dists,
880  dists < (1 + max_dist_rad) ** 2))
881 
882  def _test_rotation_agreement(self, rot_vects, max_dist_rad):
883  """ Test this rotation against the previous N found and return
884  the number that a agree within tolerance to where our test
885  points are.
886 
887  Parameters
888  ----------
889  rot_vects : `numpy.ndarray`, (N, 3)
890  Arrays of rotated 3 vectors representing the maximum x, y,
891  z extent on the unit sphere of the input source objects rotated by
892  the candidate rotations into the reference frame.
893  max_dist_rad : `float`
894  maximum distance in radians to consider two points "agreeing" on
895  a rotation
896 
897  Returns
898  -------
899  tot_consent : `int`
900  Number of candidate rotations that agree for all of the rotated
901  test 3 vectors.
902  """
903 
904  self.loglog.debug("Comparing pattern %i to previous %i rotations...",
905  rot_vects[-1][-1], len(rot_vects) - 1)
906 
907  tot_consent = 0
908  for rot_idx in range(max((len(rot_vects) - 1), 0)):
909  tmp_dist_list = []
910  for vect_idx in range(len(rot_vects[rot_idx]) - 1):
911  tmp_delta_vect = (rot_vects[rot_idx][vect_idx]
912  - rot_vects[-1][vect_idx])
913  tmp_dist_list.append(
914  np.dot(tmp_delta_vect, tmp_delta_vect))
915  if np.all(np.array(tmp_dist_list) < max_dist_rad ** 2):
916  tot_consent += 1
917  return tot_consent
918 
919  def _final_verify(self,
920  source_array,
921  shift_rot_matrix,
922  max_dist_rad,
923  min_matches):
924  """Match the all sources into the reference catalog using the shift/rot
925  matrix.
926 
927  After the initial shift/rot matrix is found, we refit the shift/rot
928  matrix using the matches the initial matrix produces to find a more
929  stable solution.
930 
931  Parameters
932  ----------
933  source_array : `numpy.ndarray` (N, 3)
934  3-vector positions on the unit-sphere representing the sources to
935  match
936  shift_rot_matrix : `numpy.ndarray` (3, 3)
937  Rotation matrix representing inferred shift/rotation of the
938  sources onto the reference catalog. Matrix need not be unitary.
939  max_dist_rad : `float`
940  Maximum distance allowed for a match.
941  min_matches : `int`
942  Minimum number of matched objects required to consider the
943  match good.
944 
945  Returns
946  -------
947  output_struct : `lsst.pipe.base.Struct`
948  Result struct with components:
949 
950  - ``match_ids`` : Pairs of indexes into the source and reference
951  data respectively defining a match (`numpy.ndarray`, (N, 2)).
952  - ``distances_rad`` : distances to between the matched objects in
953  the shift/rotated frame. (`numpy.ndarray`, (N,)).
954  - ``max_dist_rad`` : Value of the max matched distance. Either
955  returning the input value of the 2 sigma clipped value of the
956  shift/rotated distances. (`float`).
957  """
958  output_struct = pipeBase.Struct(
959  match_ids=None,
960  distances_rad=None,
961  max_dist_rad=None,
962  )
963 
964  # Perform an iterative final verify.
965  match_sources_struct = self._match_sources_match_sources(source_array,
966  shift_rot_matrix)
967  cut_ids = match_sources_struct.match_ids[
968  match_sources_struct.distances_rad < max_dist_rad]
969 
970  n_matched = len(cut_ids)
971  clipped_struct = self._clip_distances_clip_distances(
972  match_sources_struct.distances_rad)
973  n_matched_clipped = clipped_struct.n_matched_clipped
974 
975  if n_matched < min_matches or n_matched_clipped < min_matches:
976  return output_struct
977 
978  # The shift/rotation matrix returned by
979  # ``_construct_pattern_and_shift_rot_matrix``, above, was
980  # based on only six points. Here, we refine that result by
981  # using all of the good matches from the “final verification”
982  # step, above. This will produce a more consistent result.
983  fit_shift_rot_matrix = least_squares(
984  _rotation_matrix_chi_sq,
985  x0=shift_rot_matrix.flatten(),
986  args=(source_array[cut_ids[:, 0], :3],
987  self._reference_array_reference_array[cut_ids[:, 1], :3],
988  max_dist_rad)
989  ).x.reshape((3, 3))
990 
991  # Redo the matching using the newly fit shift/rotation matrix.
992  match_sources_struct = self._match_sources_match_sources(
993  source_array, fit_shift_rot_matrix)
994 
995  # Double check the match distances to make sure enough matches
996  # survive still. We'll just overwrite the previously used variables.
997  n_matched = np.sum(
998  match_sources_struct.distances_rad < max_dist_rad)
999  clipped_struct = self._clip_distances_clip_distances(
1000  match_sources_struct.distances_rad)
1001  n_matched_clipped = clipped_struct.n_matched_clipped
1002  clipped_max_dist = clipped_struct.clipped_max_dist
1003 
1004  if n_matched < min_matches or n_matched_clipped < min_matches:
1005  return output_struct
1006 
1007  # Store our matches in the output struct. Decide between the clipped
1008  # distance and the input max dist based on which is smaller.
1009  output_struct.match_ids = match_sources_struct.match_ids
1010  output_struct.distances_rad = match_sources_struct.distances_rad
1011  if clipped_max_dist < max_dist_rad:
1012  output_struct.max_dist_rad = clipped_max_dist
1013  else:
1014  output_struct.max_dist_rad = max_dist_rad
1015 
1016  return output_struct
1017 
1018  def _clip_distances(self, distances_rad):
1019  """Compute a clipped max distance and calculate the number of pairs
1020  that pass the clipped dist.
1021 
1022  Parameters
1023  ----------
1024  distances_rad : `numpy.ndarray`, (N,)
1025  Distances between pairs.
1026 
1027  Returns
1028  -------
1029  output_struct : `lsst.pipe.base.Struct`
1030  Result struct with components:
1031 
1032  - ``n_matched_clipped`` : Number of pairs that survive the
1033  clipping on distance. (`float`)
1034  - ``clipped_max_dist`` : Maximum distance after clipping.
1035  (`float`).
1036  """
1037  clipped_dists, _, clipped_max_dist = sigmaclip(
1038  distances_rad,
1039  low=100,
1040  high=2)
1041  # Check clipped distances. The minimum value here
1042  # prevents over convergence on perfect test data.
1043  if clipped_max_dist < 1e-16:
1044  clipped_max_dist = 1e-16
1045  n_matched_clipped = np.sum(distances_rad < clipped_max_dist)
1046  else:
1047  n_matched_clipped = len(clipped_dists)
1048 
1049  return pipeBase.Struct(n_matched_clipped=n_matched_clipped,
1050  clipped_max_dist=clipped_max_dist)
1051 
1052  def _match_sources(self,
1053  source_array,
1054  shift_rot_matrix):
1055  """ Shift both the reference and source catalog to the the respective
1056  frames and find their nearest neighbor using a kdTree.
1057 
1058  Removes all matches who do not agree when either the reference or
1059  source catalog is rotated. Cuts on a maximum distance are left to an
1060  external function.
1061 
1062  Parameters
1063  ----------
1064  source_array : `numpy.ndarray`, (N, 3)
1065  array of 3 vectors representing the source objects we are trying
1066  to match into the source catalog.
1067  shift_rot_matrix : `numpy.ndarray`, (3, 3)
1068  3x3 rotation matrix that performs the spherical rotation from the
1069  source frame into the reference frame.
1070 
1071  Returns
1072  -------
1073  results : `lsst.pipe.base.Struct`
1074  Result struct with components:
1075 
1076  - ``matches`` : array of integer ids into the source and
1077  reference arrays. Matches are only returned for those that
1078  satisfy the distance and handshake criteria
1079  (`numpy.ndarray`, (N, 2)).
1080  - ``distances`` : Distances between each match in radians after
1081  the shift and rotation is applied (`numpy.ndarray`, (N)).
1082  """
1083  shifted_references = np.dot(
1084  np.linalg.inv(shift_rot_matrix),
1085  self._reference_array_reference_array.transpose()).transpose()
1086  shifted_sources = np.dot(
1087  shift_rot_matrix,
1088  source_array.transpose()).transpose()
1089 
1090  ref_matches = np.empty((len(shifted_references), 2),
1091  dtype="uint16")
1092  src_matches = np.empty((len(shifted_sources), 2),
1093  dtype="uint16")
1094 
1095  ref_matches[:, 1] = np.arange(len(shifted_references),
1096  dtype="uint16")
1097  src_matches[:, 0] = np.arange(len(shifted_sources),
1098  dtype="uint16")
1099 
1100  ref_kdtree = cKDTree(self._reference_array_reference_array)
1101  src_kdtree = cKDTree(source_array)
1102 
1103  ref_to_src_dist, tmp_ref_to_src_idx = \
1104  src_kdtree.query(shifted_references)
1105  src_to_ref_dist, tmp_src_to_ref_idx = \
1106  ref_kdtree.query(shifted_sources)
1107 
1108  ref_matches[:, 0] = tmp_ref_to_src_idx
1109  src_matches[:, 1] = tmp_src_to_ref_idx
1110 
1111  handshake_mask = self._handshake_match_handshake_match(src_matches, ref_matches)
1112  return pipeBase.Struct(
1113  match_ids=src_matches[handshake_mask],
1114  distances_rad=src_to_ref_dist[handshake_mask],)
1115 
1116  def _handshake_match(self, matches_src, matches_ref):
1117  """Return only those matches where both the source
1118  and reference objects agree they they are each others'
1119  nearest neighbor.
1120 
1121  Parameters
1122  ----------
1123  matches_src : `numpy.ndarray`, (N, 2)
1124  int array of nearest neighbor matches between shifted and
1125  rotated reference objects matched into the sources.
1126  matches_ref : `numpy.ndarray`, (N, 2)
1127  int array of nearest neighbor matches between shifted and
1128  rotated source objects matched into the references.
1129  Return
1130  ------
1131  handshake_mask_array : `numpy.ndarray`, (N,)
1132  Return the array positions where the two match catalogs agree.
1133  """
1134  handshake_mask_array = np.zeros(len(matches_src), dtype=bool)
1135 
1136  for src_match_idx, match in enumerate(matches_src):
1137  ref_match_idx = np.searchsorted(matches_ref[:, 1], match[1])
1138  if match[0] == matches_ref[ref_match_idx, 0]:
1139  handshake_mask_array[src_match_idx] = True
1140  return handshake_mask_array
int max
def _intermediate_verify(self, src_pattern, ref_pattern, shift_rot_matrix, max_dist_rad)
def _test_rotation(self, src_center, ref_center, src_delta, ref_delta, cos_shift, max_cos_rot_sq)
def _construct_pattern_and_shift_rot_matrix(self, src_pattern_array, n_match, max_cos_theta_shift, max_cos_rot_sq, max_dist_rad)
def _final_verify(self, source_array, shift_rot_matrix, max_dist_rad, min_matches)
def match(self, source_array, n_check, n_match, n_agree, max_n_patterns, max_shift, max_rotation, max_dist, min_matches, pattern_skip_array=None)
def _create_shift_rot_matrix(self, cos_rot_sq, shift_matrix, src_delta, ref_ctr, ref_delta)
def _intermediate_verify_comparison(self, pattern_a, pattern_b, shift_rot_matrix, max_dist_rad)
std::pair< size_t, size_t > find_candidate_reference_pair_range(float src_dist, ndarray::Array< float, 1, 1 > const &ref_dist_array, float max_dist_rad)
Find the range of reference spokes within a spoke distance tolerance of our source spoke.
std::vector< std::pair< size_t, size_t > > create_pattern_spokes(ndarray::Array< double, 1, 1 > const &src_ctr, ndarray::Array< double, 2, 1 > const &src_delta_array, ndarray::Array< double, 1, 1 > const &src_dist_array, ndarray::Array< double, 1, 1 > const &ref_ctr, ndarray::Array< double, 1, 1 > const &proj_ref_ctr_delta, ndarray::Array< float, 1, 1 > const &ref_dist_array, ndarray::Array< uint16_t, 1, 1 > const &ref_id_array, ndarray::Array< double, 2, 1 > const &reference_array, double max_dist_rad, size_t n_match)
Create the individual spokes that make up the pattern now that the shift and rotation are within tole...