160 def match(self, source_array, n_check, n_match, n_agree,
161 max_n_patterns, max_shift, max_rotation, max_dist,
162 min_matches, pattern_skip_array=None):
163 """Match a given source catalog into the loaded reference catalog.
165 Given array of points on the unit sphere and tolerances, we
166 attempt to match a pinwheel like pattern between these input sources
167 and the reference objects this class was created with. This pattern
168 informs of the shift and rotation needed to align the input source
169 objects into the frame of the references.
173 source_array : `numpy.ndarray`, (N, 3)
174 An array of spherical x,y,z coordinates and a magnitude in units
175 of objects having a lower value for sorting. The array should be
178 Number of sources to create a pattern from. Not all objects may be
179 checked if n_match criteria is before looping through all n_check
182 Number of objects to use in constructing a pattern to match.
184 Number of found patterns that must agree on their shift and
185 rotation before exiting. Set this value to 1 to recover the
186 expected behavior of Optimistic Pattern Matcher B.
187 max_n_patters : `int`
188 Number of patterns to create from the input source objects to
189 attempt to match into the reference objects.
191 Maximum allowed shift to match patterns in arcseconds.
192 max_rotation : `float`
193 Maximum allowed rotation between patterns in degrees.
195 Maximum distance in arcseconds allowed between candidate spokes in
196 the source and reference objects. Also sets that maximum distance
197 in the intermediate verify, pattern shift/rotation agreement, and
199 pattern_skip_array : `int`
200 Patterns we would like to skip. This could be due to the pattern
201 being matched on a previous iteration that we now consider invalid.
202 This assumes the ordering of the source objects is the same
203 between different runs of the matcher which, assuming no object
204 has been inserted or the magnitudes have changed, it should be.
208 output_struct : `lsst.pipe.base.Struct`
209 Result struct with components
211 - ``matches`` : (N, 2) array of matched ids for pairs. Empty list if no
212 match found (`numpy.ndarray`, (N, 2) or `list`)
213 - ``distances_rad`` : Radian distances between the matched objects.
214 Empty list if no match found (`numpy.ndarray`, (N,))
215 - ``pattern_idx``: Index of matched pattern. None if no match found
217 - ``shift`` : Magnitude for the shift between the source and reference
218 objects in arcseconds. None if no match found (`float`).
221 sorted_source_array = source_array[source_array[:, -1].argsort(), :3]
222 n_source = len(sorted_source_array)
225 output_match_struct = pipeBase.Struct(
233 self.
log.warning(
"Source object array is empty. Unable to match. Exiting matcher.")
246 max_cos_shift = np.cos(np.radians(max_shift / 3600.))
247 max_cos_rot_sq = np.cos(np.radians(max_rotation)) ** 2
248 max_dist_rad = np.radians(max_dist / 3600.)
252 for pattern_idx
in range(np.min((max_n_patterns,
253 n_source - n_match))):
257 if pattern_skip_array
is not None and \
258 np.any(pattern_skip_array == pattern_idx):
260 "Skipping previously matched bad pattern %i...",
264 pattern = sorted_source_array[
265 pattern_idx: np.min((pattern_idx + n_check, n_source)), :3]
270 construct_return_struct = \
272 pattern, n_match, max_cos_shift, max_cos_rot_sq,
276 if construct_return_struct.ref_candidates
is None or \
277 construct_return_struct.shift_rot_matrix
is None or \
278 construct_return_struct.cos_shift
is None or \
279 construct_return_struct.sin_rot
is None:
283 ref_candidates = construct_return_struct.ref_candidates
284 shift_rot_matrix = construct_return_struct.shift_rot_matrix
285 cos_shift = construct_return_struct.cos_shift
286 sin_rot = construct_return_struct.sin_rot
290 if len(ref_candidates) < n_match:
296 tmp_rot_vect_list = []
297 for test_vect
in test_vectors:
298 tmp_rot_vect_list.append(np.dot(shift_rot_matrix, test_vect))
305 tmp_rot_vect_list.append(pattern_idx)
306 rot_vect_list.append(tmp_rot_vect_list)
319 if match_struct.match_ids
is None or \
320 match_struct.distances_rad
is None or \
321 match_struct.max_dist_rad
is None:
325 shift = np.degrees(np.arccos(cos_shift)) * 3600.
327 self.
log.debug(
"Succeeded after %i patterns.", pattern_idx)
328 self.
log.debug(
"\tShift %.4f arcsec", shift)
329 self.
log.debug(
"\tRotation: %.4f deg",
330 np.degrees(np.arcsin(sin_rot)))
333 output_match_struct.match_ids = \
334 match_struct.match_ids
335 output_match_struct.distances_rad = \
336 match_struct.distances_rad
337 output_match_struct.pattern_idx = pattern_idx
338 output_match_struct.shift = shift
339 output_match_struct.max_dist_rad = match_struct.max_dist_rad
340 return output_match_struct
342 self.
log.debug(
"Failed after %i patterns.", pattern_idx + 1)
343 return output_match_struct
346 """Compute spherical 3 vectors at the edges of the x, y, z extent
347 of the input source catalog.
351 source_array : `numpy.ndarray`, (N, 3)
352 array of 3 vectors representing positions on the unit
357 test_vectors : `numpy.ndarray`, (N, 3)
358 Array of vectors representing the maximum extents in x, y, z
359 of the input source array. These are used with the rotations
360 the code finds to test for agreement from different patterns
361 when the code is running in pessimistic mode.
364 if np.any(np.logical_not(np.isfinite(source_array))):
365 self.
log.warning(
"Input source objects contain non-finite values. "
366 "This could end badly.")
367 center_vect = np.nanmean(source_array, axis=0)
371 xbtm_vect = np.array([np.min(source_array[:, 0]), center_vect[1],
372 center_vect[2]], dtype=np.float64)
373 xtop_vect = np.array([np.max(source_array[:, 0]), center_vect[1],
374 center_vect[2]], dtype=np.float64)
375 xbtm_vect /= np.sqrt(np.dot(xbtm_vect, xbtm_vect))
376 xtop_vect /= np.sqrt(np.dot(xtop_vect, xtop_vect))
378 ybtm_vect = np.array([center_vect[0], np.min(source_array[:, 1]),
379 center_vect[2]], dtype=np.float64)
380 ytop_vect = np.array([center_vect[0], np.max(source_array[:, 1]),
381 center_vect[2]], dtype=np.float64)
382 ybtm_vect /= np.sqrt(np.dot(ybtm_vect, ybtm_vect))
383 ytop_vect /= np.sqrt(np.dot(ytop_vect, ytop_vect))
385 zbtm_vect = np.array([center_vect[0], center_vect[1],
386 np.min(source_array[:, 2])], dtype=np.float64)
387 ztop_vect = np.array([center_vect[0], center_vect[1],
388 np.max(source_array[:, 2])], dtype=np.float64)
389 zbtm_vect /= np.sqrt(np.dot(zbtm_vect, zbtm_vect))
390 ztop_vect /= np.sqrt(np.dot(ztop_vect, ztop_vect))
393 return np.array([xbtm_vect, xtop_vect, ybtm_vect, ytop_vect,
394 zbtm_vect, ztop_vect])
397 n_match, max_cos_theta_shift,
398 max_cos_rot_sq, max_dist_rad):
399 """Test an input source pattern against the reference catalog.
400 Returns the candidate matched patterns and their
401 implied rotation matrices or None.
405 src_pattern_array : `numpy.ndarray`, (N, 3)
406 Sub selection of source 3 vectors to create a pattern from
408 Number of points to attempt to create a pattern from. Must be
409 >= len(src_pattern_array)
410 max_cos_theta_shift : `float`
411 Maximum shift allowed between two patterns' centers.
412 max_cos_rot_sq : `float`
413 Maximum rotation between two patterns that have been shifted
414 to have their centers on top of each other.
415 max_dist_rad : `float`
416 Maximum delta distance allowed between the source and reference
417 pair distances to consider the reference pair a candidate for
418 the source pair. Also sets the tolerance between the opening
419 angles of the spokes when compared to the reference.
423 output_matched_pattern : `lsst.pipe.base.Struct`
424 Result struct with components:
426 - ``ref_candidates`` : ids of the matched pattern in the internal
427 reference_array object (`list` of `int`).
428 - ``src_candidates`` : Pattern ids of the sources matched
430 - ``shift_rot_matrix_src_to_ref`` : 3x3 matrix specifying the full
431 shift and rotation between the reference and source objects.
432 Rotates source into reference frame. `None` if match is not
433 found. (`numpy.ndarray`, (3, 3))
434 - ``shift_rot_matrix_ref_to_src`` : 3x3 matrix specifying the full
435 shift and rotation of the reference and source objects. Rotates
436 reference into source frame. None if match is not found
437 (`numpy.ndarray`, (3, 3)).
438 - ``cos_shift`` : Magnitude of the shift found between the two
439 patten centers. `None` if match is not found (`float`).
440 - ``sin_rot`` : float value of the rotation to align the already
441 shifted source pattern to the reference pattern. `None` if no match
446 src_delta_array = np.empty((len(src_pattern_array) - 1, 3))
447 src_delta_array[:, 0] = (src_pattern_array[1:, 0]
448 - src_pattern_array[0, 0])
449 src_delta_array[:, 1] = (src_pattern_array[1:, 1]
450 - src_pattern_array[0, 1])
451 src_delta_array[:, 2] = (src_pattern_array[1:, 2]
452 - src_pattern_array[0, 2])
453 src_dist_array = np.sqrt(src_delta_array[:, 0]**2
454 + src_delta_array[:, 1]**2
455 + src_delta_array[:, 2]**2)
457 pattern_result = construct_pattern_and_shift_rot_matrix(
458 src_pattern_array, src_delta_array, src_dist_array,
460 max_cos_theta_shift, max_cos_rot_sq, max_dist_rad)
462 if pattern_result.success:
463 candidate_array = np.array(pattern_result.candidate_pairs)
465 src_pattern_array[candidate_array[:, 1]],
467 pattern_result.shift_rot_matrix, max_dist_rad)
468 return pipeBase.Struct(ref_candidates=candidate_array[:, 0].tolist(),
469 src_candidates=candidate_array[:, 1].tolist(),
470 shift_rot_matrix=fit_shift_rot_matrix,
471 cos_shift=pattern_result.cos_shift,
472 sin_rot=pattern_result.sin_rot)
473 return pipeBase.Struct(ref_candidates=[],
475 shift_rot_matrix=
None,
481 """ Perform an intermediate verify step.
482 Rotate the matches references into the source frame and test their
483 distances against tolerance. Only return true if all points are within
488 src_pattern : `numpy.ndarray`, (N,3)
489 Array of 3 vectors representing the points that make up our source
491 ref_pattern : `numpy.ndarray`, (N,3)
492 Array of 3 vectors representing our candidate reference pinwheel
494 shift_rot_matrix : `numpy.ndarray`, (3,3)
495 3x3 rotation matrix that takes the source objects and rotates them
496 onto the frame of the reference objects
497 max_dist_rad : `float`
498 Maximum distance allowed to consider two objects the same.
502 fit_shift_rot_matrix : `numpy.ndarray`, (3,3)
503 Fitted shift/rotation matrix if all of the points in our
504 source pattern are within max_dist_rad of their matched reference
505 objects. Returns None if this criteria is not satisfied.
507 if len(src_pattern) != len(ref_pattern):
509 "Source pattern length does not match ref pattern.\n"
510 "\t source pattern len=%i, reference pattern len=%i" %
511 (len(src_pattern), len(ref_pattern)))
514 src_pattern, ref_pattern, shift_rot_matrix, max_dist_rad):
522 fit_shift_rot_matrix = least_squares(
523 _rotation_matrix_chi_sq,
524 x0=shift_rot_matrix.flatten(),
525 args=(src_pattern, ref_pattern, max_dist_rad)
529 src_pattern, ref_pattern, fit_shift_rot_matrix,
531 return fit_shift_rot_matrix
537 """Construct a generalized 3D rotation matrix about a given
542 rot_axis : `numpy.ndarray`, (3,)
543 3 vector defining the axis to rotate about.
544 cos_rotation : `float`
545 cosine of the rotation angle.
546 sin_rotation : `float`
547 sine of the rotation angle.
551 shift_matrix : `numpy.ndarray`, (3, 3)
552 3x3 spherical, rotation matrix.
554 rot_cross_matrix = np.array(
555 [[0., -rot_axis[2], rot_axis[1]],
556 [rot_axis[2], 0., -rot_axis[0]],
557 [-rot_axis[1], rot_axis[0], 0.]], dtype=np.float64)
558 shift_matrix = (cos_rotation*np.identity(3)
559 + sin_rotion*rot_cross_matrix
560 + (1. - cos_rotation)*np.outer(rot_axis, rot_axis))
565 shift_rot_matrix, max_dist_rad):
566 """Test the input rotation matrix against one input pattern and
569 If every point in the pattern after rotation is within a distance of
570 max_dist_rad to its candidate point in the other pattern, we return
575 pattern_a : `numpy.ndarray`, (N,3)
576 Array of 3 vectors representing the points that make up our source
578 pattern_b : `numpy.ndarray`, (N,3)
579 Array of 3 vectors representing our candidate reference pinwheel
581 shift_rot_matrix : `numpy.ndarray`, (3,3)
582 3x3 rotation matrix that takes the source objects and rotates them
583 onto the frame of the reference objects
584 max_dist_rad : `float`
585 Maximum distance allowed to consider two objects the same.
591 True if all rotated source points are within max_dist_rad of
592 the candidate references matches.
594 shifted_pattern_a = np.dot(shift_rot_matrix,
595 pattern_a.transpose()).transpose()
596 tmp_delta_array = shifted_pattern_a - pattern_b
597 tmp_dist_array = (tmp_delta_array[:, 0] ** 2
598 + tmp_delta_array[:, 1] ** 2
599 + tmp_delta_array[:, 2] ** 2)
600 return np.all(tmp_dist_array < max_dist_rad ** 2)
603 """ Test that the all vectors in a pattern are unit length within
606 This is useful for assuring the non unitary transforms do not contain
611 test_pattern : `numpy.ndarray`, (N, 3)
612 Test vectors at the maximum and minimum x, y, z extents.
613 max_dist_rad : `float`
614 maximum distance in radians to consider two points "agreeing" on
622 dists = (test_pattern[:, 0] ** 2
623 + test_pattern[:, 1] ** 2
624 + test_pattern[:, 2] ** 2)
626 np.logical_and((1 - max_dist_rad) ** 2 < dists,
627 dists < (1 + max_dist_rad) ** 2))
630 """ Test this rotation against the previous N found and return
631 the number that a agree within tolerance to where our test
636 rot_vects : `numpy.ndarray`, (N, 3)
637 Arrays of rotated 3 vectors representing the maximum x, y,
638 z extent on the unit sphere of the input source objects rotated by
639 the candidate rotations into the reference frame.
640 max_dist_rad : `float`
641 maximum distance in radians to consider two points "agreeing" on
647 Number of candidate rotations that agree for all of the rotated
650 self.
log.debug(
"Comparing pattern %i to previous %i rotations...",
651 rot_vects[-1][-1], len(rot_vects) - 1)
654 for rot_idx
in range(
max((len(rot_vects) - 1), 0)):
656 for vect_idx
in range(len(rot_vects[rot_idx]) - 1):
657 tmp_delta_vect = (rot_vects[rot_idx][vect_idx]
658 - rot_vects[-1][vect_idx])
659 tmp_dist_list.append(
660 np.dot(tmp_delta_vect, tmp_delta_vect))
661 if np.all(np.array(tmp_dist_list) < max_dist_rad ** 2):
670 """Match the all sources into the reference catalog using the shift/rot
673 After the initial shift/rot matrix is found, we refit the shift/rot
674 matrix using the matches the initial matrix produces to find a more
679 source_array : `numpy.ndarray` (N, 3)
680 3-vector positions on the unit-sphere representing the sources to
682 shift_rot_matrix : `numpy.ndarray` (3, 3)
683 Rotation matrix representing inferred shift/rotation of the
684 sources onto the reference catalog. Matrix need not be unitary.
685 max_dist_rad : `float`
686 Maximum distance allowed for a match.
688 Minimum number of matched objects required to consider the
693 output_struct : `lsst.pipe.base.Struct`
694 Result struct with components:
696 - ``match_ids`` : Pairs of indexes into the source and reference
697 data respectively defining a match (`numpy.ndarray`, (N, 2)).
698 - ``distances_rad`` : distances to between the matched objects in
699 the shift/rotated frame. (`numpy.ndarray`, (N,)).
700 - ``max_dist_rad`` : Value of the max matched distance. Either
701 returning the input value of the 2 sigma clipped value of the
702 shift/rotated distances. (`float`).
704 output_struct = pipeBase.Struct(
713 cut_ids = match_sources_struct.match_ids[
714 match_sources_struct.distances_rad < max_dist_rad]
716 n_matched = len(cut_ids)
718 match_sources_struct.distances_rad)
719 n_matched_clipped = clipped_struct.n_matched_clipped
721 if n_matched < min_matches
or n_matched_clipped < min_matches:
729 fit_shift_rot_matrix = least_squares(
730 _rotation_matrix_chi_sq,
731 x0=shift_rot_matrix.flatten(),
732 args=(source_array[cut_ids[:, 0], :3],
739 source_array, fit_shift_rot_matrix)
744 match_sources_struct.distances_rad < max_dist_rad)
746 match_sources_struct.distances_rad)
747 n_matched_clipped = clipped_struct.n_matched_clipped
748 clipped_max_dist = clipped_struct.clipped_max_dist
750 if n_matched < min_matches
or n_matched_clipped < min_matches:
755 output_struct.match_ids = match_sources_struct.match_ids
756 output_struct.distances_rad = match_sources_struct.distances_rad
757 if clipped_max_dist < max_dist_rad:
758 output_struct.max_dist_rad = clipped_max_dist
760 output_struct.max_dist_rad = max_dist_rad
765 """Compute a clipped max distance and calculate the number of pairs
766 that pass the clipped dist.
770 distances_rad : `numpy.ndarray`, (N,)
771 Distances between pairs.
775 output_struct : `lsst.pipe.base.Struct`
776 Result struct with components:
778 - ``n_matched_clipped`` : Number of pairs that survive the
779 clipping on distance. (`float`)
780 - ``clipped_max_dist`` : Maximum distance after clipping.
783 clipped_dists, _, clipped_max_dist = sigmaclip(
789 if clipped_max_dist < 1e-16:
790 clipped_max_dist = 1e-16
791 n_matched_clipped = np.sum(distances_rad < clipped_max_dist)
793 n_matched_clipped = len(clipped_dists)
795 return pipeBase.Struct(n_matched_clipped=n_matched_clipped,
796 clipped_max_dist=clipped_max_dist)
801 """ Shift both the reference and source catalog to the the respective
802 frames and find their nearest neighbor using a kdTree.
804 Removes all matches who do not agree when either the reference or
805 source catalog is rotated. Cuts on a maximum distance are left to an
810 source_array : `numpy.ndarray`, (N, 3)
811 array of 3 vectors representing the source objects we are trying
812 to match into the source catalog.
813 shift_rot_matrix : `numpy.ndarray`, (3, 3)
814 3x3 rotation matrix that performs the spherical rotation from the
815 source frame into the reference frame.
819 results : `lsst.pipe.base.Struct`
820 Result struct with components:
822 - ``matches`` : array of integer ids into the source and
823 reference arrays. Matches are only returned for those that
824 satisfy the distance and handshake criteria
825 (`numpy.ndarray`, (N, 2)).
826 - ``distances`` : Distances between each match in radians after
827 the shift and rotation is applied (`numpy.ndarray`, (N)).
829 shifted_references = np.dot(
830 np.linalg.inv(shift_rot_matrix),
832 shifted_sources = np.dot(
834 source_array.transpose()).transpose()
836 ref_matches = np.empty((len(shifted_references), 2),
838 src_matches = np.empty((len(shifted_sources), 2),
841 ref_matches[:, 1] = np.arange(len(shifted_references),
843 src_matches[:, 0] = np.arange(len(shifted_sources),
847 src_kdtree = cKDTree(source_array)
849 ref_to_src_dist, tmp_ref_to_src_idx = \
850 src_kdtree.query(shifted_references)
851 src_to_ref_dist, tmp_src_to_ref_idx = \
852 ref_kdtree.query(shifted_sources)
854 ref_matches[:, 0] = tmp_ref_to_src_idx
855 src_matches[:, 1] = tmp_src_to_ref_idx
858 return pipeBase.Struct(
859 match_ids=src_matches[handshake_mask],
860 distances_rad=src_to_ref_dist[handshake_mask],)
863 """Return only those matches where both the source
864 and reference objects agree they they are each others'
869 matches_src : `numpy.ndarray`, (N, 2)
870 int array of nearest neighbor matches between shifted and
871 rotated reference objects matched into the sources.
872 matches_ref : `numpy.ndarray`, (N, 2)
873 int array of nearest neighbor matches between shifted and
874 rotated source objects matched into the references.
878 handshake_mask_array : `numpy.ndarray`, (N,)
879 Array positions where the two match catalogs agree.
881 handshake_mask_array = np.zeros(len(matches_src), dtype=bool)
883 for src_match_idx, match
in enumerate(matches_src):
884 ref_match_idx = np.searchsorted(matches_ref[:, 1], match[1])
885 if match[0] == matches_ref[ref_match_idx, 0]:
886 handshake_mask_array[src_match_idx] =
True
887 return handshake_mask_array