From 36b734bb7ce1ceaad409c367d75507e3350fda50 Mon Sep 17 00:00:00 2001 From: Juan Escudero Date: Tue, 29 Oct 2024 11:21:24 +0100 Subject: [PATCH] Fix quad matching (#149) --- iop4lib/instruments/dipol.py | 43 ++++++++++++++++++++++++++--------- iop4lib/utils/quadmatching.py | 38 +++++++++++++++++++++++++++++-- 2 files changed, 68 insertions(+), 13 deletions(-) diff --git a/iop4lib/instruments/dipol.py b/iop4lib/instruments/dipol.py index 268d0b09..bbb46db0 100644 --- a/iop4lib/instruments/dipol.py +++ b/iop4lib/instruments/dipol.py @@ -753,7 +753,10 @@ def PolyArea(x,y): logger.debug(f"Selected the quads [{best_i},{best_j}]") - logger.debug(f"t = {t_L[np.argmin(delta_t)]}") + t = t_L[np.argmin(delta_t)] + R = R_L[np.argmin(delta_t)] + logger.debug(f"t = {t}, R = {R}") + logger.debug(f"det R = {np.linalg.det(R)}") # build_summary_images, replace indices_selected by all_indices if no R,t filtering was done. if summary_kwargs['build_summary_images']: @@ -792,23 +795,42 @@ def PolyArea(x,y): quads_1 = [qorder_ish(quad) for quad in quads_1] quads_2 = [qorder_ish(quad) for quad in quads_2] - # get the pre wcs with the target in the center of the image + # save the flipped status of both images + + is_redf_pol_flipped = 'FLIPSTAT' in redf_pol.rawfit.header and redf_pol.rawfit.header['FLIPSTAT'] == "Flip" + is_redf_phot_flipped = 'FLIPSTAT' in redf_phot.rawfit.header and redf_phot.rawfit.header['FLIPSTAT'] == "Flip" + + # get the pre wcs with the target in the center of the image (if the image is flipped, the angle is negative) angle_mean, angle_std = get_angle_from_history(redf_pol, target_src) - if 'FLIPSTAT' in redf_pol.rawfit.header: - # TODO: check that indeed the mere presence of this keyword means that the image is flipped, without the need of checking the value. - # FLIPSTAT is a MaximDL thing only, but it seems that the iamge is flipped whenever the keyword is present, regardless of the value. + if is_redf_pol_flipped: angle = - angle_mean else: angle = angle_mean logger.debug(f"Using {angle=} for pre wcs.") - pre_wcs = build_wcs_centered_on((redf_pol.width//2,redf_pol.height//2), redf=redf_phot, angle=angle) + # Get the appropiate transformation depending on whether both images are flipped or not + + from iop4lib.utils.quadmatching import find_best_transformation, distance_to_y_flip, distance_to_identity - # get the pixel arrays in the polarimetry field and in the FULL photometry field to relate them + if is_redf_pol_flipped != is_redf_phot_flipped: + R, t, perm = find_best_transformation(quads_1[best_i], quads_2[best_j], distance_to_y_flip) + else: + R, t, perm = find_best_transformation(quads_1[best_i], quads_2[best_j], distance_to_identity) + + # fit a wcs centered on the target source + + pre_wcs = build_wcs_centered_on((redf_pol.width//2,redf_pol.height//2), redf=redf_phot, angle=angle) + + # Get the pixel position of the quad points in the (small) polarimetry field pix_array_1 = np.array(list(zip(*[(x,y) for x,y in quads_1[best_i]]))) - pix_array_2 = np.array(list(zip(*[(x+x_start,y+y_start) for x,y in quads_2[best_j]]))) + + # Get pixel positions of the quad points in the (full) photometry field + #pix_array_2 = np.array(list(zip(*[(x+x_start,y+y_start) for x,y in quads_2[best_j]]))) + # instead of quads_2[best_j], transform quads_1[best_i] with the linear transformation, to avoid + # incorrect permutations + pix_array_2 = np.array(list(zip(*[(x+x_start,y+y_start) for x,y in (np.dot(R, np.array(quads_1[best_i]).T).T + t)]))) # fit the WCS so the pixel arrays in 1 correspond to the ra/dec of the pixel array in 2 wcs1 = fit_wcs_from_points(pix_array_1, redf_phot.wcs1.pixel_to_world(*pix_array_2), projection=pre_wcs) @@ -944,8 +966,7 @@ def _build_wcs_for_polarimetry_from_target_O_and_E(cls, redf: 'ReducedFit', summ logger.debug(f"Using angle {angle_mean:.2f} +- {angle_std:.2f} deg") # DIPOL polarimery images seem to be flipped vertically, which results in negative angle - # TODO: watch this FLIP thing, check that indeed this is the behaviour - if 'FLIPSTAT' in redf.rawfit.header: + if 'FLIPSTAT' in redf.rawfit.header and redf.rawfit.header['FLIPSTAT'] == "Flip": angle = - angle_mean else: angle = angle_mean @@ -978,7 +999,7 @@ def estimate_common_apertures(cls, reducedfits, reductionmethod=None, fit_boxsiz sigma = fit_res_dict['sigma'] fwhm = fit_res_dict["mean_fwhm"] - return 1.1*fwhm, 6*fwhm, 10*fwhm, fit_res_dict + return 1.1*fwhm, 6*fwhm, 10*fwhm, fit_res_dict @classmethod def get_instrumental_polarization(cls, reducedfit) -> dict: diff --git a/iop4lib/utils/quadmatching.py b/iop4lib/utils/quadmatching.py index 7f62a8b8..3e425f94 100644 --- a/iop4lib/utils/quadmatching.py +++ b/iop4lib/utils/quadmatching.py @@ -1,4 +1,4 @@ -from itertools import combinations +import itertools import numpy as np def distance(h1,h2): @@ -186,4 +186,38 @@ def find_linear_transformation(P1, P2): # Translation Vector t = P2_mean - R @ P1_mean - return R, t \ No newline at end of file + return R, t + +def find_best_transformation(P1, P2, distance_function): + """Find the linear transformation of P1 to P2 that minimizes the distance function, searching over all permutations of P2.""" + + P1, P2 = np.array(P1), np.array(P2) + + # Store all permutations and their corresponding transformations (R, t) + permutations = list(itertools.permutations(P2)) + transformations = [find_linear_transformation(P1, perm) for perm in permutations] + + # Compute distances to the target transformation (using distance_function) + distances = [distance_function(R) for R, _ in transformations] + + # Find the permutation with the minimum distance + best_index = np.argmin(distances) + + best_R, best_t = transformations[best_index] + best_perm = permutations[best_index] + + return best_R, best_t, best_perm + +def distance_to_y_flip(R): + # Target matrix for a single Y-axis flip + R_target = np.array([[1, 0], [0, -1]]) + + # Compute Frobenius norm of the difference + return np.linalg.norm(R - R_target) + +def distance_to_identity(R): + # Target matrix is the identity matrix + R_target = np.eye(2) # Identity matrix of size 2x2 + + # Compute Frobenius norm of the difference + return np.linalg.norm(R - R_target) \ No newline at end of file