3 from scipy.optimize
import least_squares
4 from scipy.spatial
import cKDTree
5 from scipy.stats
import sigmaclip
7 from .pessimisticPatternMatcherUtils
import (find_candidate_reference_pair_range,
8 create_pattern_spokes,)
12 def _rotation_matrix_chi_sq(flattened_rot_matrix,
16 """Compute the squared differences for least squares fitting.
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.
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.
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.
46 rot_matrix = flattened_rot_matrix.reshape((3, 3))
48 rot_pattern_a = np.dot(rot_matrix, pattern_a.transpose()).transpose()
49 diff_pattern_a_to_b = rot_pattern_a - pattern_b
52 return diff_pattern_a_to_b.flatten() / max_dist_rad
56 """Class implementing a pessimistic version of Optimistic Pattern Matcher
57 B (OPMb) from Tabur 2007. See `DMTN-031 <http://ls.st/DMTN-031`_
61 reference_array : `numpy.ndarray`, (N, 3)
62 spherical points x, y, z of to use as reference objects for
64 log : `lsst.log.Log` or `logging.Logger`
65 Logger for outputting debug info.
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
87 raise ValueError(
"No reference objects supplied")
89 def _build_distances_and_angles(self):
90 """Create the data structures we will use to search for our pattern
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.
111 endIdx = startIdx + self.
_n_reference_n_reference - 1 - ref_id
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,
123 self.
_dist_array_dist_array[startIdx:endIdx] = np.sqrt(
125 - ref_obj) ** 2).sum(axis=1))
131 sorted_dist_args = self.
_dist_array_dist_array.argsort()
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.
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.
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
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
157 Number of objects to use in constructing a pattern to match.
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.
166 Maximum allowed shift to match patterns in arcseconds.
167 max_rotation : `float`
168 Maximum allowed rotation between patterns in degrees.
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
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.
183 output_struct : `lsst.pipe.base.Struct`
184 Result struct with components
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
192 - ``shift`` : Magnitude for the shift between the source and reference
193 objects in arcseconds. None if no match found (`float`).
197 sorted_source_array = source_array[source_array[:, -1].argsort(), :3]
198 n_source = len(sorted_source_array)
201 output_match_struct = pipeBase.Struct(
209 self.
loglog.
warning(
"Source object array is empty. Unable to match. Exiting matcher.")
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.)
228 for pattern_idx
in range(np.min((max_n_patterns,
229 n_source - n_match))):
233 if pattern_skip_array
is not None and \
234 np.any(pattern_skip_array == pattern_idx):
236 "Skipping previously matched bad pattern %i...",
240 pattern = sorted_source_array[
241 pattern_idx: np.min((pattern_idx + n_check, n_source)), :3]
246 construct_return_struct = \
248 pattern, n_match, max_cos_shift, max_cos_rot_sq,
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:
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
266 if len(ref_candidates) < n_match:
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))
281 tmp_rot_vect_list.append(pattern_idx)
282 rot_vect_list.append(tmp_rot_vect_list)
291 match_struct = self.
_final_verify_final_verify(source_array[:, :3],
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:
301 shift = np.degrees(np.arccos(cos_shift)) * 3600.
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)))
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
318 self.
loglog.
debug(
"Failed after %i patterns.", pattern_idx + 1)
319 return output_match_struct
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.
327 source_array : `numpy.ndarray`, (N, 3)
328 array of 3 vectors representing positions on the unit
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.
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)
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))
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))
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))
370 return np.array([xbtm_vect, xtop_vect, ybtm_vect, ytop_vect,
371 zbtm_vect, ztop_vect])
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.
378 Returns the candidate matched patterns and their
379 implied rotation matrices or None.
383 src_pattern_array : `numpy.ndarray`, (N, 3)
384 Sub selection of source 3 vectors to create a pattern from
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.
401 output_matched_pattern : `lsst.pipe.base.Struct`
402 Result struct with components:
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
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
427 output_matched_pattern = pipeBase.Struct(
430 shift_rot_matrix=
None,
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)
454 src_dist_array[0], self.
_dist_array_dist_array, max_dist_rad)
456 def generate_ref_dist_indexes(low, high):
457 """Generator to loop outward from the midpoint between two values.
462 Minimum of index range.
464 Maximum of index range.
471 mid = (high + low) // 2
472 for idx
in range(high - low):
481 for ref_dist_idx
in generate_ref_dist_indexes(candidate_range[0],
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]
490 shift_rot_matrix =
None
497 cos_shift = np.dot(src_pattern_array[0], ref_center)
498 if cos_shift < max_cos_theta_shift:
502 ref_candidates.append(ref_id)
505 ref_candidates.append(
506 tmp_ref_pair_list[1])
510 ref_candidates.append(
511 tmp_ref_pair_list[0])
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:
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
538 tmp_ref_dist_array = np.sqrt(
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]
549 src_pattern_array[0], src_delta_array, src_dist_array,
551 tmp_ref_dist_array, tmp_ref_id_array, self.
_reference_array_reference_array,
552 max_dist_rad, n_match)
556 if len(pattern_spokes) < n_match - 2:
560 ref_candidates.extend([cand[0]
for cand
in pattern_spokes])
561 src_candidates.extend([cand[1]
for cand
in pattern_spokes])
567 cos_rot_sq, shift_matrix, src_delta_array[0],
571 if shift_rot_struct.sin_rot
is None or \
572 shift_rot_struct.shift_rot_matrix
is None:
576 sin_rot = shift_rot_struct.sin_rot
577 shift_rot_matrix = shift_rot_struct.shift_rot_matrix
586 src_pattern_array[src_candidates],
588 shift_rot_matrix, max_dist_rad)
590 if fit_shift_rot_matrix
is not None:
592 output_matched_pattern.ref_candidates = ref_candidates
593 output_matched_pattern.src_candidates = src_candidates
594 output_matched_pattern.shift_rot_matrix = \
596 output_matched_pattern.cos_shift = cos_shift
597 output_matched_pattern.sin_rot = sin_rot
598 return output_matched_pattern
600 return output_matched_pattern
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.
611 src_center : `numpy.ndarray`, (N, 3)
613 ref_center : `numpy.ndarray`, (N, 3)
614 3 vector defining the center of the candidate reference pinwheel
616 src_delta : `numpy.ndarray`, (N, 3)
617 3 vector delta between the source pattern center and the end of
619 ref_delta : `numpy.ndarray`, (N, 3)
620 3 vector delta of the candidate matched reference pair
622 Cosine of the angle between the source and reference candidate
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.
631 result : `lsst.pipe.base.Struct`
632 Result struct with components:
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)).
645 elif cos_shift < -1.0:
647 sin_shift = np.sqrt(1 - cos_shift ** 2)
653 rot_axis = np.cross(src_center, ref_center)
654 rot_axis /= sin_shift
656 rot_axis, cos_shift, sin_shift)
658 shift_matrix = np.identity(3)
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)))
671 if cos_rot_sq < max_cos_rot_sq:
672 return pipeBase.Struct(
674 proj_ref_ctr_delta=
None,
678 return pipeBase.Struct(
679 cos_rot_sq=cos_rot_sq,
680 proj_ref_ctr_delta=proj_ref_delta,
681 shift_matrix=shift_matrix,)
683 def _create_spherical_rotation_matrix(self, rot_axis, cos_rotation,
685 """Construct a generalized 3D rotation matrix about a given
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.
699 shift_matrix : `numpy.ndarray`, (3, 3)
700 3x3 spherical, rotation matrix.
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))
713 def _create_shift_rot_matrix(self, cos_rot_sq, shift_matrix, src_delta,
715 """ Create the final part of our spherical rotation matrix.
720 cosine of the rotation needed to align our source and reference
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
728 ref_ctr : `numpy.ndarray`, (3,)
729 3 vector on the unit-sphere representing the center of our
731 ref_delta : `numpy.ndarray`, (3,)
732 3 vector delta made by the first pair of the reference pattern.
736 result : `lsst.pipe.base.Struct`
737 Result struct with components:
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)).
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)
750 sin_rot = np.sign(delta_dot_cross) * np.sqrt(1 - cos_rot_sq)
752 ref_ctr, cos_rot, sin_rot)
754 shift_rot_matrix = np.dot(rot_matrix, shift_matrix)
756 return pipeBase.Struct(
758 shift_rot_matrix=shift_rot_matrix,)
760 def _intermediate_verify(self, src_pattern, ref_pattern, shift_rot_matrix,
762 """ Perform an intermediate verify step.
764 Rotate the matches references into the source frame and test their
765 distances against tolerance. Only return true if all points are within
770 src_pattern : `numpy.ndarray`, (N,3)
771 Array of 3 vectors representing the points that make up our source
773 ref_pattern : `numpy.ndarray`, (N,3)
774 Array of 3 vectors representing our candidate reference pinwheel
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.
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.
789 if len(src_pattern) != len(ref_pattern):
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)))
796 src_pattern, ref_pattern, shift_rot_matrix, max_dist_rad):
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)
811 src_pattern, ref_pattern, fit_shift_rot_matrix,
813 return fit_shift_rot_matrix
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
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
828 pattern_a : `numpy.ndarray`, (N,3)
829 Array of 3 vectors representing the points that make up our source
831 pattern_b : `numpy.ndarray`, (N,3)
832 Array of 3 vectors representing our candidate reference pinwheel
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.
844 True if all rotated source points are within max_dist_rad of
845 the candidate references matches.
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)
855 def _test_pattern_lengths(self, test_pattern, max_dist_rad):
856 """ Test that the all vectors in a pattern are unit length within
859 This is useful for assuring the non unitary transforms do not contain
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
875 dists = (test_pattern[:, 0] ** 2
876 + test_pattern[:, 1] ** 2
877 + test_pattern[:, 2] ** 2)
879 np.logical_and((1 - max_dist_rad) ** 2 < dists,
880 dists < (1 + max_dist_rad) ** 2))
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
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
900 Number of candidate rotations that agree for all of the rotated
904 self.
loglog.
debug(
"Comparing pattern %i to previous %i rotations...",
905 rot_vects[-1][-1], len(rot_vects) - 1)
908 for rot_idx
in range(
max((len(rot_vects) - 1), 0)):
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):
919 def _final_verify(self,
924 """Match the all sources into the reference catalog using the shift/rot
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
933 source_array : `numpy.ndarray` (N, 3)
934 3-vector positions on the unit-sphere representing the sources to
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.
942 Minimum number of matched objects required to consider the
947 output_struct : `lsst.pipe.base.Struct`
948 Result struct with components:
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`).
958 output_struct = pipeBase.Struct(
965 match_sources_struct = self.
_match_sources_match_sources(source_array,
967 cut_ids = match_sources_struct.match_ids[
968 match_sources_struct.distances_rad < max_dist_rad]
970 n_matched = len(cut_ids)
972 match_sources_struct.distances_rad)
973 n_matched_clipped = clipped_struct.n_matched_clipped
975 if n_matched < min_matches
or n_matched_clipped < min_matches:
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],
993 source_array, fit_shift_rot_matrix)
998 match_sources_struct.distances_rad < max_dist_rad)
1000 match_sources_struct.distances_rad)
1001 n_matched_clipped = clipped_struct.n_matched_clipped
1002 clipped_max_dist = clipped_struct.clipped_max_dist
1004 if n_matched < min_matches
or n_matched_clipped < min_matches:
1005 return output_struct
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
1014 output_struct.max_dist_rad = max_dist_rad
1016 return output_struct
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.
1024 distances_rad : `numpy.ndarray`, (N,)
1025 Distances between pairs.
1029 output_struct : `lsst.pipe.base.Struct`
1030 Result struct with components:
1032 - ``n_matched_clipped`` : Number of pairs that survive the
1033 clipping on distance. (`float`)
1034 - ``clipped_max_dist`` : Maximum distance after clipping.
1037 clipped_dists, _, clipped_max_dist = sigmaclip(
1043 if clipped_max_dist < 1e-16:
1044 clipped_max_dist = 1e-16
1045 n_matched_clipped = np.sum(distances_rad < clipped_max_dist)
1047 n_matched_clipped = len(clipped_dists)
1049 return pipeBase.Struct(n_matched_clipped=n_matched_clipped,
1050 clipped_max_dist=clipped_max_dist)
1052 def _match_sources(self,
1055 """ Shift both the reference and source catalog to the the respective
1056 frames and find their nearest neighbor using a kdTree.
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
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.
1073 results : `lsst.pipe.base.Struct`
1074 Result struct with components:
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)).
1083 shifted_references = np.dot(
1084 np.linalg.inv(shift_rot_matrix),
1086 shifted_sources = np.dot(
1088 source_array.transpose()).transpose()
1090 ref_matches = np.empty((len(shifted_references), 2),
1092 src_matches = np.empty((len(shifted_sources), 2),
1095 ref_matches[:, 1] = np.arange(len(shifted_references),
1097 src_matches[:, 0] = np.arange(len(shifted_sources),
1101 src_kdtree = cKDTree(source_array)
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)
1108 ref_matches[:, 0] = tmp_ref_to_src_idx
1109 src_matches[:, 1] = tmp_src_to_ref_idx
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],)
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'
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.
1131 handshake_mask_array : `numpy.ndarray`, (N,)
1132 Return the array positions where the two match catalogs agree.
1134 handshake_mask_array = np.zeros(len(matches_src), dtype=bool)
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
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 _test_rotation_agreement(self, rot_vects, max_dist_rad)
def _create_spherical_rotation_matrix(self, rot_axis, cos_rotation, sin_rotion)
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 _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)
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...