Fix quad matching (#149)

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:

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
+    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)