3 from scipy.optimize
import least_squares
4 from scipy.spatial
import cKDTree
5 from scipy.stats
import sigmaclip
10 def _rotation_matrix_chi_sq(flattened_rot_matrix,
14 """Compute the squared differences for least squares fitting.
16 Given a flattened rotation matrix, one N point pattern and another N point
17 pattern to transform into to, compute the squared differences between the
18 points in the two patterns after the rotation.
22 flattened_rot_matrix : `numpy.ndarray`, (9, )
23 A flattened array representing a 3x3 rotation matrix. The array is
24 flattened to comply with the API of scipy.optimize.least_squares.
25 Flattened elements are [[0, 0], [0, 1], [0, 2], [1, 0]...]
26 pattern_a : `numpy.ndarray`, (N, 3)
27 A array containing N, 3 vectors representing the objects we would like
28 to transform into the frame of pattern_b.
29 pattern_b : `numpy.ndarray`, (N, 3)
30 A array containing N, 3 vectors representing the reference frame we
31 would like to transform pattern_a into.
32 max_dist_rad : `float`
33 The maximum distance allowed from the pattern matching. This value is
34 used as the standard error for the resultant chi values.
38 noralized_diff : `numpy.ndarray`, (9,)
39 Array of differences between the vectors representing of the source
40 pattern rotated into the reference frame and the converse. This is
41 used to minimize in a least squares fitter.
44 rot_matrix = flattened_rot_matrix.reshape((3, 3))
46 rot_pattern_a = np.dot(rot_matrix, pattern_a.transpose()).transpose()
47 diff_pattern_a_to_b = rot_pattern_a - pattern_b
50 return diff_pattern_a_to_b.flatten() / max_dist_rad
54 """Class implementing a pessimistic version of Optimistic Pattern Matcher
55 B (OPMb) from Tabur 2007. See `DMTN-031 <http://ls.st/DMTN-031`_
59 reference_array : `numpy.ndarray`, (N, 3)
60 spherical points x, y, z of to use as reference objects for
63 Logger for outputting debug info.
67 The class loads and stores the reference object
68 in a convenient data structure for matching any set of source objects that
69 are assumed to contain each other. The pessimistic nature of the algorithm
70 comes from requiring that it discovers at least two patterns that agree on
71 the correct shift and rotation for matching before exiting. The original
72 behavior of OPMb can be recovered simply. Patterns matched between the
73 input datasets are n-spoked pinwheels created from n+1 points. Refer to
74 DMTN #031 for more details. http://github.com/lsst-dm/dmtn-031
85 raise ValueError(
"No reference objects supplied")
87 def _build_distances_and_angles(self):
88 """Create the data structures we will use to search for our pattern
91 Throughout this function and the rest of the class we use id to
92 reference the position in the input reference catalog and index to
93 'index' into the arrays sorted on distance.
107 sub_id_array_list = []
108 sub_dist_array_list = []
116 sub_id_array = np.zeros((self.
_n_reference_n_reference - 1 - ref_id, 2),
118 sub_id_array[:, 0] = ref_id
119 sub_id_array[:, 1] = np.arange(ref_id + 1, self.
_n_reference_n_reference,
125 - ref_obj).astype(np.float32)
126 sub_dist_array = np.sqrt(sub_delta_array[:, 0] ** 2
127 + sub_delta_array[:, 1] ** 2
128 + sub_delta_array[:, 2] ** 2)
132 sub_id_array_list.append(sub_id_array)
133 sub_dist_array_list.append(sub_dist_array)
136 self.
_pair_id_array_pair_id_array[ref_id, ref_id:] = sub_id_array[:, 1]
148 sorted_pair_dist_args = self.
_pair_dist_array_pair_dist_array[ref_id, :].argsort()
150 ref_id, sorted_pair_dist_args]
152 ref_id, sorted_pair_dist_args]
155 unsorted_id_array = np.concatenate(sub_id_array_list)
156 unsorted_dist_array = np.concatenate(sub_dist_array_list)
160 sorted_dist_args = unsorted_dist_array.argsort()
161 self.
_dist_array_dist_array = unsorted_dist_array[sorted_dist_args]
162 self.
_id_array_id_array = unsorted_id_array[sorted_dist_args]
166 def match(self, source_array, n_check, n_match, n_agree,
167 max_n_patterns, max_shift, max_rotation, max_dist,
168 min_matches, pattern_skip_array=None):
169 """Match a given source catalog into the loaded reference catalog.
171 Given array of points on the unit sphere and tolerances, we
172 attempt to match a pinwheel like pattern between these input sources
173 and the reference objects this class was created with. This pattern
174 informs of the shift and rotation needed to align the input source
175 objects into the frame of the references.
179 source_array : `numpy.ndarray`, (N, 3)
180 An array of spherical x,y,z coordinates and a magnitude in units
181 of objects having a lower value for sorting. The array should be
184 Number of sources to create a pattern from. Not all objects may be
185 checked if n_match criteria is before looping through all n_check
188 Number of objects to use in constructing a pattern to match.
190 Number of found patterns that must agree on their shift and
191 rotation before exiting. Set this value to 1 to recover the
192 expected behavior of Optimistic Pattern Matcher B.
193 max_n_patters : `int`
194 Number of patterns to create from the input source objects to
195 attempt to match into the reference objects.
197 Maximum allowed shift to match patterns in arcseconds.
198 max_rotation : `float`
199 Maximum allowed rotation between patterns in degrees.
201 Maximum distance in arcseconds allowed between candidate spokes in
202 the source and reference objects. Also sets that maximum distance
203 in the intermediate verify, pattern shift/rotation agreement, and
205 pattern_skip_array : `int`
206 Patterns we would like to skip. This could be due to the pattern
207 being matched on a previous iteration that we now consider invalid.
208 This assumes the ordering of the source objects is the same
209 between different runs of the matcher which, assuming no object
210 has been inserted or the magnitudes have changed, it should be.
214 output_struct : `lsst.pipe.base.Struct`
215 Result struct with components
217 - ``matches`` : (N, 2) array of matched ids for pairs. Empty list if no
218 match found (`numpy.ndarray`, (N, 2) or `list`)
219 - ``distances_rad`` : Radian distances between the matched objects.
220 Empty list if no match found (`numpy.ndarray`, (N,))
221 - ``pattern_idx``: Index of matched pattern. None if no match found
223 - ``shift`` : Magnitude for the shift between the source and reference
224 objects in arcseconds. None if no match found (`float`).
228 sorted_source_array = source_array[source_array[:, -1].argsort(), :3]
229 n_source = len(sorted_source_array)
232 output_match_struct = pipeBase.Struct(
240 self.
loglog.
warn(
"Source object array is empty. Unable to match. "
254 max_cos_shift = np.cos(np.radians(max_shift / 3600.))
255 max_cos_rot_sq = np.cos(np.radians(max_rotation)) ** 2
256 max_dist_rad = np.radians(max_dist / 3600.)
260 for pattern_idx
in range(np.min((max_n_patterns,
261 n_source - n_match))):
265 if pattern_skip_array
is not None and \
266 np.any(pattern_skip_array == pattern_idx):
268 "Skipping previously matched bad pattern %i..." %
272 pattern = sorted_source_array[
273 pattern_idx: np.min((pattern_idx + n_check, n_source)), :3]
278 construct_return_struct = \
280 pattern, n_match, max_cos_shift, max_cos_rot_sq,
284 if construct_return_struct.ref_candidates
is None or \
285 construct_return_struct.shift_rot_matrix
is None or \
286 construct_return_struct.cos_shift
is None or \
287 construct_return_struct.sin_rot
is None:
291 ref_candidates = construct_return_struct.ref_candidates
292 shift_rot_matrix = construct_return_struct.shift_rot_matrix
293 cos_shift = construct_return_struct.cos_shift
294 sin_rot = construct_return_struct.sin_rot
298 if len(ref_candidates) < n_match:
304 tmp_rot_vect_list = []
305 for test_vect
in test_vectors:
306 tmp_rot_vect_list.append(np.dot(shift_rot_matrix, test_vect))
313 tmp_rot_vect_list.append(pattern_idx)
314 rot_vect_list.append(tmp_rot_vect_list)
323 match_struct = self.
_final_verify_final_verify(source_array[:, :3],
327 if match_struct.match_ids
is None or \
328 match_struct.distances_rad
is None or \
329 match_struct.max_dist_rad
is None:
333 shift = np.degrees(np.arccos(cos_shift)) * 3600.
335 self.
loglog.
debug(
"Succeeded after %i patterns." % pattern_idx)
336 self.
loglog.
debug(
"\tShift %.4f arcsec" % shift)
337 self.
loglog.
debug(
"\tRotation: %.4f deg" %
338 np.degrees(np.arcsin(sin_rot)))
341 output_match_struct.match_ids = \
342 match_struct.match_ids
343 output_match_struct.distances_rad = \
344 match_struct.distances_rad
345 output_match_struct.pattern_idx = pattern_idx
346 output_match_struct.shift = shift
347 output_match_struct.max_dist_rad = match_struct.max_dist_rad
348 return output_match_struct
350 self.
loglog.
debug(
"Failed after %i patterns." % (pattern_idx + 1))
351 return output_match_struct
353 def _compute_test_vectors(self, source_array):
354 """Compute spherical 3 vectors at the edges of the x, y, z extent
355 of the input source catalog.
359 source_array : `numpy.ndarray`, (N, 3)
360 array of 3 vectors representing positions on the unit
365 test_vectors : `numpy.ndarray`, (N, 3)
366 Array of vectors representing the maximum extents in x, y, z
367 of the input source array. These are used with the rotations
368 the code finds to test for agreement from different patterns
369 when the code is running in pessimistic mode.
373 if np.any(np.logical_not(np.isfinite(source_array))):
374 self.
loglog.
warn(
"Input source objects contain non-finite values. "
375 "This could end badly.")
376 center_vect = np.nanmean(source_array, axis=0)
380 xbtm_vect = np.array([np.min(source_array[:, 0]), center_vect[1],
381 center_vect[2]], dtype=np.float64)
382 xtop_vect = np.array([np.max(source_array[:, 0]), center_vect[1],
383 center_vect[2]], dtype=np.float64)
384 xbtm_vect /= np.sqrt(np.dot(xbtm_vect, xbtm_vect))
385 xtop_vect /= np.sqrt(np.dot(xtop_vect, xtop_vect))
387 ybtm_vect = np.array([center_vect[0], np.min(source_array[:, 1]),
388 center_vect[2]], dtype=np.float64)
389 ytop_vect = np.array([center_vect[0], np.max(source_array[:, 1]),
390 center_vect[2]], dtype=np.float64)
391 ybtm_vect /= np.sqrt(np.dot(ybtm_vect, ybtm_vect))
392 ytop_vect /= np.sqrt(np.dot(ytop_vect, ytop_vect))
394 zbtm_vect = np.array([center_vect[0], center_vect[1],
395 np.min(source_array[:, 2])], dtype=np.float64)
396 ztop_vect = np.array([center_vect[0], center_vect[1],
397 np.max(source_array[:, 2])], dtype=np.float64)
398 zbtm_vect /= np.sqrt(np.dot(zbtm_vect, zbtm_vect))
399 ztop_vect /= np.sqrt(np.dot(ztop_vect, ztop_vect))
402 return np.array([xbtm_vect, xtop_vect, ybtm_vect, ytop_vect,
403 zbtm_vect, ztop_vect])
405 def _construct_pattern_and_shift_rot_matrix(self, src_pattern_array,
406 n_match, max_cos_theta_shift,
407 max_cos_rot_sq, max_dist_rad):
408 """Test an input source pattern against the reference catalog.
410 Returns the candidate matched patterns and their
411 implied rotation matrices or None.
415 src_pattern_array : `numpy.ndarray`, (N, 3)
416 Sub selection of source 3 vectors to create a pattern from
418 Number of points to attempt to create a pattern from. Must be
419 >= len(src_pattern_array)
420 max_cos_theta_shift : `float`
421 Maximum shift allowed between two patterns' centers.
422 max_cos_rot_sq : `float`
423 Maximum rotation between two patterns that have been shifted
424 to have their centers on top of each other.
425 max_dist_rad : `float`
426 Maximum delta distance allowed between the source and reference
427 pair distances to consider the reference pair a candidate for
428 the source pair. Also sets the tolerance between the opening
429 angles of the spokes when compared to the reference.
433 output_matched_pattern : `lsst.pipe.base.Struct`
434 Result struct with components:
436 - ``ref_candidates`` : ids of the matched pattern in the internal
437 reference_array object (`list` of `int`).
438 - ``src_candidates`` : Pattern ids of the sources matched
440 - ``shift_rot_matrix_src_to_ref`` : 3x3 matrix specifying the full
441 shift and rotation between the reference and source objects.
442 Rotates source into reference frame. `None` if match is not
443 found. (`numpy.ndarray`, (3, 3))
444 - ``shift_rot_matrix_ref_to_src`` : 3x3 matrix specifying the full
445 shift and rotation of the reference and source objects. Rotates
446 reference into source frame. None if match is not found
447 (`numpy.ndarray`, (3, 3)).
448 - ``cos_shift`` : Magnitude of the shift found between the two
449 patten centers. `None` if match is not found (`float`).
450 - ``sin_rot`` : float value of the rotation to align the already
451 shifted source pattern to the reference pattern. `None` if no match
459 output_matched_pattern = pipeBase.Struct(
462 shift_rot_matrix=
None,
468 src_delta_array = np.empty((len(src_pattern_array) - 1, 3))
469 src_delta_array[:, 0] = (src_pattern_array[1:, 0]
470 - src_pattern_array[0, 0])
471 src_delta_array[:, 1] = (src_pattern_array[1:, 1]
472 - src_pattern_array[0, 1])
473 src_delta_array[:, 2] = (src_pattern_array[1:, 2]
474 - src_pattern_array[0, 2])
475 src_dist_array = np.sqrt(src_delta_array[:, 0]**2
476 + src_delta_array[:, 1]**2
477 + src_delta_array[:, 2]**2)
483 src_dist_array[0], self.
_dist_array_dist_array, max_dist_rad)
486 for ref_dist_idx
in ref_dist_index_array:
490 tmp_ref_pair_list = self.
_id_array_id_array[ref_dist_idx]
491 for pair_idx, ref_id
in enumerate(tmp_ref_pair_list):
492 src_candidates = [0, 1]
494 shift_rot_matrix =
None
501 cos_shift = np.dot(src_pattern_array[0], ref_center)
502 if cos_shift < max_cos_theta_shift:
506 ref_candidates.append(ref_id)
509 ref_candidates.append(
510 tmp_ref_pair_list[1])
514 ref_candidates.append(
515 tmp_ref_pair_list[0])
525 src_pattern_array[0], ref_center, src_delta_array[0],
526 ref_delta, cos_shift, max_cos_rot_sq)
527 if test_rot_struct.cos_rot_sq
is None or \
528 test_rot_struct.proj_ref_ctr_delta
is None or \
529 test_rot_struct.shift_matrix
is None:
533 cos_rot_sq = test_rot_struct.cos_rot_sq
534 proj_ref_ctr_delta = test_rot_struct.proj_ref_ctr_delta
535 shift_matrix = test_rot_struct.shift_matrix
546 src_pattern_array[0], src_delta_array, src_dist_array,
548 tmp_ref_dist_array, tmp_ref_id_array, max_dist_rad,
553 if len(pattern_spoke_struct.ref_spoke_list) < n_match - 2
or \
554 len(pattern_spoke_struct.src_spoke_list) < n_match - 2:
558 ref_candidates.extend(pattern_spoke_struct.ref_spoke_list)
559 src_candidates.extend(pattern_spoke_struct.src_spoke_list)
565 cos_rot_sq, shift_matrix, src_delta_array[0],
569 if shift_rot_struct.sin_rot
is None or \
570 shift_rot_struct.shift_rot_matrix
is None:
574 sin_rot = shift_rot_struct.sin_rot
575 shift_rot_matrix = shift_rot_struct.shift_rot_matrix
584 src_pattern_array[src_candidates],
586 shift_rot_matrix, max_dist_rad)
588 if fit_shift_rot_matrix
is not None:
590 output_matched_pattern.ref_candidates = ref_candidates
591 output_matched_pattern.src_candidates = src_candidates
592 output_matched_pattern.shift_rot_matrix = \
594 output_matched_pattern.cos_shift = cos_shift
595 output_matched_pattern.sin_rot = sin_rot
596 return output_matched_pattern
598 return output_matched_pattern
600 def _find_candidate_reference_pairs(self, src_dist, ref_dist_array,
602 """Wrap numpy.searchsorted to find the range of reference spokes
603 within a spoke distance tolerance of our source spoke.
605 Returns an array sorted from the smallest absolute delta distance
606 between source and reference spoke length. This sorting increases the
607 speed for the pattern search greatly.
612 float value of the distance we would like to search for in
613 the reference array in radians.
614 ref_dist_array : `numpy.ndarray`, (N,)
615 sorted array of distances in radians.
616 max_dist_rad : `float`
617 maximum plus/minus search to find in the reference array in
622 tmp_diff_array : `numpy.ndarray`, (N,)
623 indices lookup into the input ref_dist_array sorted by the
624 difference in value to the src_dist from absolute value
629 start_idx = np.searchsorted(ref_dist_array, src_dist - max_dist_rad)
630 end_idx = np.searchsorted(ref_dist_array, src_dist + max_dist_rad,
634 if start_idx == end_idx:
640 if end_idx > ref_dist_array.shape[0]:
641 end_idx = ref_dist_array.shape[0]
646 tmp_diff_array = np.fabs(ref_dist_array[start_idx:end_idx] - src_dist)
647 return tmp_diff_array.argsort() + start_idx
649 def _test_rotation(self, src_center, ref_center, src_delta, ref_delta,
650 cos_shift, max_cos_rot_sq):
651 """ Test if the rotation implied between the source
652 pattern and reference pattern is within tolerance. To test this
653 we need to create the first part of our spherical rotation matrix
654 which we also return for use later.
658 src_center : `numpy.ndarray`, (N, 3)
660 ref_center : `numpy.ndarray`, (N, 3)
661 3 vector defining the center of the candidate reference pinwheel
663 src_delta : `numpy.ndarray`, (N, 3)
664 3 vector delta between the source pattern center and the end of
666 ref_delta : `numpy.ndarray`, (N, 3)
667 3 vector delta of the candidate matched reference pair
669 Cosine of the angle between the source and reference candidate
671 max_cos_rot_sq : `float`
672 candidate reference pair after shifting the centers on top of each
673 other. The function will return None if the rotation implied is
674 greater than max_cos_rot_sq.
678 result : `lsst.pipe.base.Struct`
679 Result struct with components:
681 - ``cos_rot_sq`` : magnitude of the rotation needed to align the
682 two patterns after their centers are shifted on top of each
683 other. `None` if rotation test fails (`float`).
684 - ``shift_matrix`` : 3x3 rotation matrix describing the shift needed to
685 align the source and candidate reference center. `None` if rotation
686 test fails (`numpy.ndarray`, (N, 3)).
692 elif cos_shift < -1.0:
694 sin_shift = np.sqrt(1 - cos_shift ** 2)
700 rot_axis = np.cross(src_center, ref_center)
701 rot_axis /= sin_shift
703 rot_axis, cos_shift, sin_shift)
705 shift_matrix = np.identity(3)
709 rot_src_delta = np.dot(shift_matrix, src_delta)
710 proj_src_delta = (rot_src_delta
711 - np.dot(rot_src_delta, ref_center) * ref_center)
712 proj_ref_delta = (ref_delta
713 - np.dot(ref_delta, ref_center) * ref_center)
714 cos_rot_sq = (np.dot(proj_src_delta, proj_ref_delta) ** 2
715 / (np.dot(proj_src_delta, proj_src_delta)
716 * np.dot(proj_ref_delta, proj_ref_delta)))
718 if cos_rot_sq < max_cos_rot_sq:
719 return pipeBase.Struct(
721 proj_ref_ctr_delta=
None,
725 return pipeBase.Struct(
726 cos_rot_sq=cos_rot_sq,
727 proj_ref_ctr_delta=proj_ref_delta,
728 shift_matrix=shift_matrix,)
730 def _create_spherical_rotation_matrix(self, rot_axis, cos_rotation,
732 """Construct a generalized 3D rotation matrix about a given
737 rot_axis : `numpy.ndarray`, (3,)
738 3 vector defining the axis to rotate about.
739 cos_rotation : `float`
740 cosine of the rotation angle.
741 sin_rotation : `float`
742 sine of the rotation angle.
746 shift_matrix : `numpy.ndarray`, (3, 3)
747 3x3 spherical, rotation matrix.
750 rot_cross_matrix = np.array(
751 [[0., -rot_axis[2], rot_axis[1]],
752 [rot_axis[2], 0., -rot_axis[0]],
753 [-rot_axis[1], rot_axis[0], 0.]], dtype=np.float64)
754 shift_matrix = (cos_rotation*np.identity(3)
755 + sin_rotion*rot_cross_matrix
756 + (1. - cos_rotation)*np.outer(rot_axis, rot_axis))
760 def _create_pattern_spokes(self, src_ctr, src_delta_array, src_dist_array,
761 ref_ctr, ref_ctr_id, proj_ref_ctr_delta,
762 ref_dist_array, ref_id_array, max_dist_rad,
764 """ Create the individual spokes that make up the pattern now that the
765 shift and rotation are within tolerance.
767 If we can't create a valid pattern we exit early.
771 src_ctr : `numpy.ndarray`, (3,)
772 3 vector of the source pinwheel center
773 src_delta_array : `numpy.ndarray`, (N, 3)
774 Array of 3 vector deltas between the source center and the pairs
775 that make up the remaining spokes of the pinwheel
776 src_dist_array : `numpy.ndarray`, (N, 3)
777 Array of the distances of each src_delta in the pinwheel
778 ref_ctr : `numpy.ndarray`, (3,)
779 3 vector of the candidate reference center
781 id of the ref_ctr in the master reference array
782 proj_ref_ctr_delta : `numpy.ndarray`, (3,)
783 Plane projected 3 vector formed from the center point of the
784 candidate pin-wheel and the second point in the pattern to create
785 the first spoke pair. This is the candidate pair that was matched
786 in the main _construct_pattern_and_shift_rot_matrix loop
787 ref_dist_array : `numpy.ndarray`, (N,)
788 Array of vector distances for each of the reference pairs
789 ref_id_array : `numpy.ndarray`, (N,)
790 Array of id lookups into the master reference array that our
791 center id object is paired with.
792 max_dist_rad : `float`
793 Maximum search distance
795 Number of source deltas that must be matched into the reference
796 deltas in order to consider this a successful pattern match.
800 output_spokes : `lsst.pipe.base.Struct`
801 Result struct with components:
803 - ``ref_spoke_list`` : list of ints specifying ids into the master
804 reference array (`list` of `int`).
805 - ``src_spoke_list`` : list of ints specifying indices into the
806 current source pattern that is being tested (`list` of `int`).
809 output_spokes = pipeBase.Struct(
821 proj_src_ctr_delta = (src_delta_array[0]
822 - np.dot(src_delta_array[0], src_ctr) * src_ctr)
823 proj_src_ctr_dist_sq = np.dot(proj_src_ctr_delta, proj_src_ctr_delta)
826 proj_ref_ctr_dist_sq = np.dot(proj_ref_ctr_delta, proj_ref_ctr_delta)
829 for src_idx
in range(1, len(src_dist_array)):
830 if n_fail > len(src_dist_array) - (n_match - 1):
835 src_sin_tol = (max_dist_rad
836 / (src_dist_array[src_idx] + max_dist_rad))
843 if src_sin_tol > max_sin_tol:
844 src_sin_tol = max_sin_tol
849 src_delta_array[src_idx]
850 - np.dot(src_delta_array[src_idx], src_ctr) * src_ctr)
851 geom_dist_src = np.sqrt(
852 np.dot(proj_src_delta, proj_src_delta)
853 * proj_src_ctr_dist_sq)
856 cos_theta_src = (np.dot(proj_src_delta, proj_src_ctr_delta)
858 cross_src = (np.cross(proj_src_delta, proj_src_ctr_delta)
860 sin_theta_src = np.dot(cross_src, src_ctr)
865 src_dist_array[src_idx], ref_dist_array, max_dist_rad)
875 proj_ref_ctr_dist_sq,
885 ref_spoke_list.append(ref_id)
886 src_spoke_list.append(src_idx + 1)
890 if len(ref_spoke_list) >= n_match - 2:
892 output_spokes.ref_spoke_list = ref_spoke_list
893 output_spokes.src_spoke_list = src_spoke_list
898 def _test_spoke(self, cos_theta_src, sin_theta_src, ref_ctr, ref_ctr_id,
899 proj_ref_ctr_delta, proj_ref_ctr_dist_sq,
900 ref_dist_idx_array, ref_id_array, src_sin_tol):
901 """Test the opening angle between the first spoke of our pattern
902 for the source object against the reference object.
904 This method makes heavy use of the small angle approximation to perform
909 cos_theta_src : `float`
910 Cosine of the angle between the current candidate source spoke and
912 sin_theta_src : `float`
913 Sine of the angle between the current candidate source spoke and
915 ref_ctr : `numpy.ndarray`, (3,)
916 3 vector of the candidate reference center
918 id lookup of the ref_ctr into the master reference array
919 proj_ref_ctr_delta : `float`
920 Plane projected first spoke in the reference pattern using the
921 pattern center as normal.
922 proj_ref_ctr_dist_sq : `float`
923 Squared length of the projected vector.
924 ref_dist_idx_array : `numpy.ndarray`, (N,)
925 Indices sorted by the delta distance between the source
926 spoke we are trying to test and the candidate reference
928 ref_id_array : `numpy.ndarray`, (N,)
929 Array of id lookups into the master reference array that our
930 center id object is paired with.
931 src_sin_tol : `float`
932 Sine of tolerance allowed between source and reference spoke
938 If we can not find a candidate spoke we return `None` else we
939 return an int id into the master reference array.
943 for ref_dist_idx
in ref_dist_idx_array:
945 ref_delta = (self.
_reference_array_reference_array[ref_id_array[ref_dist_idx]]
949 proj_ref_delta = ref_delta - np.dot(ref_delta, ref_ctr) * ref_ctr
950 geom_dist_ref = np.sqrt(proj_ref_ctr_dist_sq
951 * np.dot(proj_ref_delta, proj_ref_delta))
952 cos_theta_ref = (np.dot(proj_ref_delta, proj_ref_ctr_delta)
957 if cos_theta_ref ** 2 < (1 - src_sin_tol ** 2):
958 cos_sq_comparison = ((cos_theta_src - cos_theta_ref) ** 2
959 / (1 - cos_theta_ref ** 2))
961 cos_sq_comparison = ((cos_theta_src - cos_theta_ref) ** 2
966 if cos_sq_comparison > src_sin_tol ** 2:
972 cross_ref = (np.cross(proj_ref_delta, proj_ref_ctr_delta)
974 sin_theta_ref = np.dot(cross_ref, ref_ctr)
978 if abs(cos_theta_src) < src_sin_tol:
979 sin_comparison = (sin_theta_src - sin_theta_ref) / src_sin_tol
982 (sin_theta_src - sin_theta_ref) / cos_theta_ref
984 if abs(sin_comparison) > src_sin_tol:
988 return ref_id_array[ref_dist_idx]
992 def _create_shift_rot_matrix(self, cos_rot_sq, shift_matrix, src_delta,
994 """ Create the final part of our spherical rotation matrix.
999 cosine of the rotation needed to align our source and reference
1001 shift_matrix : `numpy.ndarray`, (3, 3)
1002 3x3 rotation matrix for shifting the source pattern center on top
1003 of the candidate reference pattern center.
1004 src_delta : `numpy.ndarray`, (3,)
1005 3 vector delta of representing the first spoke of the source
1007 ref_ctr : `numpy.ndarray`, (3,)
1008 3 vector on the unit-sphere representing the center of our
1010 ref_delta : `numpy.ndarray`, (3,)
1011 3 vector delta made by the first pair of the reference pattern.
1015 result : `lsst.pipe.base.Struct`
1016 Result struct with components:
1018 - ``sin_rot`` : float sine of the amount of rotation between the
1019 source and reference pattern. We use sine here as it is
1020 signed and tells us the chirality of the rotation (`float`).
1021 - ``shift_rot_matrix`` : float array representing the 3x3 rotation
1022 matrix that takes the source pattern and shifts and rotates
1023 it to align with the reference pattern (`numpy.ndarray`, (3,3)).
1025 cos_rot = np.sqrt(cos_rot_sq)
1026 rot_src_delta = np.dot(shift_matrix, src_delta)
1027 delta_dot_cross = np.dot(np.cross(rot_src_delta, ref_delta), ref_ctr)
1029 sin_rot = np.sign(delta_dot_cross) * np.sqrt(1 - cos_rot_sq)
1031 ref_ctr, cos_rot, sin_rot)
1033 shift_rot_matrix = np.dot(rot_matrix, shift_matrix)
1035 return pipeBase.Struct(
1037 shift_rot_matrix=shift_rot_matrix,)
1039 def _intermediate_verify(self, src_pattern, ref_pattern, shift_rot_matrix,
1041 """ Perform an intermediate verify step.
1043 Rotate the matches references into the source frame and test their
1044 distances against tolerance. Only return true if all points are within
1049 src_pattern : `numpy.ndarray`, (N,3)
1050 Array of 3 vectors representing the points that make up our source
1052 ref_pattern : `numpy.ndarray`, (N,3)
1053 Array of 3 vectors representing our candidate reference pinwheel
1055 shift_rot_matrix : `numpy.ndarray`, (3,3)
1056 3x3 rotation matrix that takes the source objects and rotates them
1057 onto the frame of the reference objects
1058 max_dist_rad : `float`
1059 Maximum distance allowed to consider two objects the same.
1063 fit_shift_rot_matrix : `numpy.ndarray`, (3,3)
1064 Return the fitted shift/rotation matrix if all of the points in our
1065 source pattern are within max_dist_rad of their matched reference
1066 objects. Returns None if this criteria is not satisfied.
1068 if len(src_pattern) != len(ref_pattern):
1070 "Source pattern length does not match ref pattern.\n"
1071 "\t source pattern len=%i, reference pattern len=%i" %
1072 (len(src_pattern), len(ref_pattern)))
1075 src_pattern, ref_pattern, shift_rot_matrix, max_dist_rad):
1083 fit_shift_rot_matrix = least_squares(
1084 _rotation_matrix_chi_sq,
1085 x0=shift_rot_matrix.flatten(),
1086 args=(src_pattern, ref_pattern, max_dist_rad)
1090 src_pattern, ref_pattern, fit_shift_rot_matrix,
1092 return fit_shift_rot_matrix
1096 def _intermediate_verify_comparison(self, pattern_a, pattern_b,
1097 shift_rot_matrix, max_dist_rad):
1098 """Test the input rotation matrix against one input pattern and
1101 If every point in the pattern after rotation is within a distance of
1102 max_dist_rad to its candidate point in the other pattern, we return
1107 pattern_a : `numpy.ndarray`, (N,3)
1108 Array of 3 vectors representing the points that make up our source
1110 pattern_b : `numpy.ndarray`, (N,3)
1111 Array of 3 vectors representing our candidate reference pinwheel
1113 shift_rot_matrix : `numpy.ndarray`, (3,3)
1114 3x3 rotation matrix that takes the source objects and rotates them
1115 onto the frame of the reference objects
1116 max_dist_rad : `float`
1117 Maximum distance allowed to consider two objects the same.
1123 True if all rotated source points are within max_dist_rad of
1124 the candidate references matches.
1126 shifted_pattern_a = np.dot(shift_rot_matrix,
1127 pattern_a.transpose()).transpose()
1128 tmp_delta_array = shifted_pattern_a - pattern_b
1129 tmp_dist_array = (tmp_delta_array[:, 0] ** 2
1130 + tmp_delta_array[:, 1] ** 2
1131 + tmp_delta_array[:, 2] ** 2)
1132 return np.all(tmp_dist_array < max_dist_rad ** 2)
1134 def _test_pattern_lengths(self, test_pattern, max_dist_rad):
1135 """ Test that the all vectors in a pattern are unit length within
1138 This is useful for assuring the non unitary transforms do not contain
1139 too much distortion.
1143 test_pattern : `numpy.ndarray`, (N, 3)
1144 Test vectors at the maximum and minimum x, y, z extents.
1145 max_dist_rad : `float`
1146 maximum distance in radians to consider two points "agreeing" on
1154 dists = (test_pattern[:, 0] ** 2
1155 + test_pattern[:, 1] ** 2
1156 + test_pattern[:, 2] ** 2)
1158 np.logical_and((1 - max_dist_rad) ** 2 < dists,
1159 dists < (1 + max_dist_rad) ** 2))
1161 def _test_rotation_agreement(self, rot_vects, max_dist_rad):
1162 """ Test this rotation against the previous N found and return
1163 the number that a agree within tolerance to where our test
1168 rot_vects : `numpy.ndarray`, (N, 3)
1169 Arrays of rotated 3 vectors representing the maximum x, y,
1170 z extent on the unit sphere of the input source objects rotated by
1171 the candidate rotations into the reference frame.
1172 max_dist_rad : `float`
1173 maximum distance in radians to consider two points "agreeing" on
1179 Number of candidate rotations that agree for all of the rotated
1183 self.
loglog.
debug(
"Comparing pattern %i to previous %i rotations..." %
1184 (rot_vects[-1][-1], len(rot_vects) - 1))
1187 for rot_idx
in range(
max((len(rot_vects) - 1), 0)):
1189 for vect_idx
in range(len(rot_vects[rot_idx]) - 1):
1190 tmp_delta_vect = (rot_vects[rot_idx][vect_idx]
1191 - rot_vects[-1][vect_idx])
1192 tmp_dist_list.append(
1193 np.dot(tmp_delta_vect, tmp_delta_vect))
1194 if np.all(np.array(tmp_dist_list) < max_dist_rad ** 2):
1198 def _final_verify(self,
1203 """Match the all sources into the reference catalog using the shift/rot
1206 After the initial shift/rot matrix is found, we refit the shift/rot
1207 matrix using the matches the initial matrix produces to find a more
1212 source_array : `numpy.ndarray` (N, 3)
1213 3-vector positions on the unit-sphere representing the sources to
1215 shift_rot_matrix : `numpy.ndarray` (3, 3)
1216 Rotation matrix representing inferred shift/rotation of the
1217 sources onto the reference catalog. Matrix need not be unitary.
1218 max_dist_rad : `float`
1219 Maximum distance allowed for a match.
1221 Minimum number of matched objects required to consider the
1226 output_struct : `lsst.pipe.base.Struct`
1227 Result struct with components:
1229 - ``match_ids`` : Pairs of indexes into the source and reference
1230 data respectively defining a match (`numpy.ndarray`, (N, 2)).
1231 - ``distances_rad`` : distances to between the matched objects in
1232 the shift/rotated frame. (`numpy.ndarray`, (N,)).
1233 - ``max_dist_rad`` : Value of the max matched distance. Either
1234 returning the input value of the 2 sigma clipped value of the
1235 shift/rotated distances. (`float`).
1237 output_struct = pipeBase.Struct(
1244 match_sources_struct = self.
_match_sources_match_sources(source_array,
1246 cut_ids = match_sources_struct.match_ids[
1247 match_sources_struct.distances_rad < max_dist_rad]
1249 n_matched = len(cut_ids)
1251 match_sources_struct.distances_rad)
1252 n_matched_clipped = clipped_struct.n_matched_clipped
1254 if n_matched < min_matches
or n_matched_clipped < min_matches:
1255 return output_struct
1262 fit_shift_rot_matrix = least_squares(
1263 _rotation_matrix_chi_sq,
1264 x0=shift_rot_matrix.flatten(),
1265 args=(source_array[cut_ids[:, 0], :3],
1272 source_array, fit_shift_rot_matrix)
1277 match_sources_struct.distances_rad < max_dist_rad)
1279 match_sources_struct.distances_rad)
1280 n_matched_clipped = clipped_struct.n_matched_clipped
1281 clipped_max_dist = clipped_struct.clipped_max_dist
1283 if n_matched < min_matches
or n_matched_clipped < min_matches:
1284 return output_struct
1288 output_struct.match_ids = match_sources_struct.match_ids
1289 output_struct.distances_rad = match_sources_struct.distances_rad
1290 if clipped_max_dist < max_dist_rad:
1291 output_struct.max_dist_rad = clipped_max_dist
1293 output_struct.max_dist_rad = max_dist_rad
1295 return output_struct
1297 def _clip_distances(self, distances_rad):
1298 """Compute a clipped max distance and calculate the number of pairs
1299 that pass the clipped dist.
1303 distances_rad : `numpy.ndarray`, (N,)
1304 Distances between pairs.
1308 output_struct : `lsst.pipe.base.Struct`
1309 Result struct with components:
1311 - ``n_matched_clipped`` : Number of pairs that survive the
1312 clipping on distance. (`float`)
1313 - ``clipped_max_dist`` : Maximum distance after clipping.
1316 clipped_dists, _, clipped_max_dist = sigmaclip(
1322 if clipped_max_dist < 1e-16:
1323 clipped_max_dist = 1e-16
1324 n_matched_clipped = np.sum(distances_rad < clipped_max_dist)
1326 n_matched_clipped = len(clipped_dists)
1328 return pipeBase.Struct(n_matched_clipped=n_matched_clipped,
1329 clipped_max_dist=clipped_max_dist)
1331 def _match_sources(self,
1334 """ Shift both the reference and source catalog to the the respective
1335 frames and find their nearest neighbor using a kdTree.
1337 Removes all matches who do not agree when either the reference or
1338 source catalog is rotated. Cuts on a maximum distance are left to an
1343 source_array : `numpy.ndarray`, (N, 3)
1344 array of 3 vectors representing the source objects we are trying
1345 to match into the source catalog.
1346 shift_rot_matrix : `numpy.ndarray`, (3, 3)
1347 3x3 rotation matrix that performs the spherical rotation from the
1348 source frame into the reference frame.
1352 results : `lsst.pipe.base.Struct`
1353 Result struct with components:
1355 - ``matches`` : array of integer ids into the source and
1356 reference arrays. Matches are only returned for those that
1357 satisfy the distance and handshake criteria
1358 (`numpy.ndarray`, (N, 2)).
1359 - ``distances`` : Distances between each match in radians after
1360 the shift and rotation is applied (`numpy.ndarray`, (N)).
1362 shifted_references = np.dot(
1363 np.linalg.inv(shift_rot_matrix),
1365 shifted_sources = np.dot(
1367 source_array.transpose()).transpose()
1369 ref_matches = np.empty((len(shifted_references), 2),
1371 src_matches = np.empty((len(shifted_sources), 2),
1374 ref_matches[:, 1] = np.arange(len(shifted_references),
1376 src_matches[:, 0] = np.arange(len(shifted_sources),
1380 src_kdtree = cKDTree(source_array)
1382 ref_to_src_dist, tmp_ref_to_src_idx = \
1383 src_kdtree.query(shifted_references)
1384 src_to_ref_dist, tmp_src_to_ref_idx = \
1385 ref_kdtree.query(shifted_sources)
1387 ref_matches[:, 0] = tmp_ref_to_src_idx
1388 src_matches[:, 1] = tmp_src_to_ref_idx
1390 handshake_mask = self.
_handshake_match_handshake_match(src_matches, ref_matches)
1391 return pipeBase.Struct(
1392 match_ids=src_matches[handshake_mask],
1393 distances_rad=src_to_ref_dist[handshake_mask],)
1395 def _handshake_match(self, matches_src, matches_ref):
1396 """Return only those matches where both the source
1397 and reference objects agree they they are each others'
1402 matches_src : `numpy.ndarray`, (N, 2)
1403 int array of nearest neighbor matches between shifted and
1404 rotated reference objects matched into the sources.
1405 matches_ref : `numpy.ndarray`, (N, 2)
1406 int array of nearest neighbor matches between shifted and
1407 rotated source objects matched into the references.
1410 handshake_mask_array : `numpy.ndarray`, (N,)
1411 Return the array positions where the two match catalogs agree.
1413 handshake_mask_array = np.zeros(len(matches_src), dtype=bool)
1415 for src_match_idx, match
in enumerate(matches_src):
1416 ref_match_idx = np.searchsorted(matches_ref[:, 1], match[1])
1417 if match[0] == matches_ref[ref_match_idx, 0]:
1418 handshake_mask_array[src_match_idx] =
True
1419 return handshake_mask_array
def _intermediate_verify(self, src_pattern, ref_pattern, shift_rot_matrix, max_dist_rad)
def _build_distances_and_angles(self)
def _clip_distances(self, distances_rad)
def _test_rotation(self, src_center, ref_center, src_delta, ref_delta, cos_shift, max_cos_rot_sq)
def _handshake_match(self, matches_src, matches_ref)
def _test_pattern_lengths(self, test_pattern, max_dist_rad)
def _create_pattern_spokes(self, src_ctr, src_delta_array, src_dist_array, ref_ctr, ref_ctr_id, proj_ref_ctr_delta, ref_dist_array, ref_id_array, max_dist_rad, n_match)
def _test_rotation_agreement(self, rot_vects, max_dist_rad)
def _create_spherical_rotation_matrix(self, rot_axis, cos_rotation, sin_rotion)
def _find_candidate_reference_pairs(self, src_dist, ref_dist_array, max_dist_rad)
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 _match_sources(self, source_array, shift_rot_matrix)
def __init__(self, reference_array, log)
def _test_spoke(self, cos_theta_src, sin_theta_src, ref_ctr, ref_ctr_id, proj_ref_ctr_delta, proj_ref_ctr_dist_sq, ref_dist_idx_array, ref_id_array, src_sin_tol)
def _create_shift_rot_matrix(self, cos_rot_sq, shift_matrix, src_delta, ref_ctr, ref_delta)
def _compute_test_vectors(self, source_array)
def _intermediate_verify_comparison(self, pattern_a, pattern_b, shift_rot_matrix, max_dist_rad)
Angle abs(Angle const &a)