Quad matching in DIPOL polarimetry images

%autoawait off
%matplotlib inline
%run 01_notebook_configuration.py

# To reduce disk usage, only a central cut of of the full field of the DIPOL 
# camera is actually saved.

# polarimetry (cut)
from IPython.display import display
from iop4lib.db import RawFit, ReducedFit
from iop4lib.utils.plotting import imshow_w_sources
rf_pol = RawFit.by_fileloc(f"OSN-T090/2023-11-06/BLLAC_R_IAR-0760.fts")
display(rf_pol)
imshow_w_sources(rf_pol.mdata)

RawFit(id=79):
- telescope: OSN-T090
- night: 2023-11-06
- filename: BLLAC_R_IAR-0760.fts
- instrument: DIPOL
- imgtype: LIGHT
- imgbinning: 2x2
- size: 1104x896
- obsmode: POLARIMETRY
- band: R
- exptime: 20.0
- flags: DOWNLOADED,CLASSIFIED,BUILT_REDUCED

The position of the cut is defined in the header of the raw file

rf_pol.header['XORGSUBF'], rf_pol.header['YORGSUBF']

(1496, 1000)

We can compare it with the full field image of the photometry field

# photometry filed (full field)
rf_phot = RawFit.by_fileloc(f"OSN-T090/2023-11-06/BLLac_IAR-0001R.fit")
display(rf_phot)
imshow_w_sources(rf_phot.mdata)

RawFit(id=77):
- telescope: OSN-T090
- night: 2023-11-06
- filename: BLLac_IAR-0001R.fit
- instrument: DIPOL
- imgtype: LIGHT
- imgbinning: 2x2
- size: 4144x2822
- obsmode: PHOTOMETRY
- band: R
- exptime: 300.0
- flags: DOWNLOADED,CLASSIFIED,BUILT_REDUCED

Let’s build the photometry field first

redf_phot = ReducedFit.create(rf_phot)
redf_phot.build_file()

2026-06-23 15:23 - pid 3911 - [reducedfit.py:123] - DEBUG - DB entry of ReducedFit for OSN-T090/2023-11-06/BLLac_IAR-0001R.fit already exists, it will be used instead.
2026-06-23 15:23 - pid 3911 - [instrument.py:433] - WARNING - <ReducedFit 77 | OSN-T090/2023-11-06/BLLac_IAR-0001R.fit>: masterflat in this epoch could not be found, attemptying adjacent epochs.
2026-06-23 15:23 - pid 3911 - [instrument.py:411] - WARNING - Master Flat from epoch <Epoch 4 | OSN-T090/2023-12-13> is more than 1 week away from epoch <Epoch 6 | OSN-T090/2023-11-06>.
2026-06-23 15:23 - pid 3911 - [instrument.py:540] - DEBUG - <ReducedFit 77 | OSN-T090/2023-11-06/BLLac_IAR-0001R.fit>: building file
2026-06-23 15:23 - pid 3911 - [dipol.py:248] - DEBUG - <ReducedFit 77 | OSN-T090/2023-11-06/BLLac_IAR-0001R.fit>: applying masters
2026-06-23 15:23 - pid 3911 - [instrument.py:546] - DEBUG - <ReducedFit 77 | OSN-T090/2023-11-06/BLLac_IAR-0001R.fit>: performing astrometric calibration
2026-06-23 15:23 - pid 3911 - [astrometry.py:214] - DEBUG - <ReducedFit 77 | OSN-T090/2023-11-06/BLLac_IAR-0001R.fit>: 9 different combinations of parameters to try.
2026-06-23 15:23 - pid 3911 - [astrometry.py:217] - DEBUG - <ReducedFit 77 | OSN-T090/2023-11-06/BLLac_IAR-0001R.fit>: attempt 1 / 9, ({'keep_n_seg': 300, 'border_margin_px': 20, 'output_logodds_threshold': 14, 'n_rms_seg': 1.2, 'bkg_filter_size': 11, 'bkg_box_size': 32, 'seg_fwhm': 1.0, 'npixels': 32, 'seg_kernel_size': None, 'position_hint': PositionHint(ra_deg=330.6804125, dec_deg=42.27766361111111, radius_deg=0.23137333333333335), 'size_hint': SizeHint(lower_arcsec_per_pixel=0.1273, upper_arcsec_per_pixel=0.14070000000000002), 'disp_sign_mean': array([-209.,   13.]), 'disp_sign_err': array([12., 16.])}) ...
2026-06-23 15:23 - pid 3911 - [astrometry.py:299] - DEBUG - <ReducedFit 77 | OSN-T090/2023-11-06/BLLac_IAR-0001R.fit>: len(pos_seg)=220
2026-06-23 15:23 - pid 3911 - [astrometry.py:306] - DEBUG - <ReducedFit 77 | OSN-T090/2023-11-06/BLLac_IAR-0001R.fit>: Removing segments within 20 px from border.
2026-06-23 15:23 - pid 3911 - [astrometry.py:316] - DEBUG - <ReducedFit 77 | OSN-T090/2023-11-06/BLLac_IAR-0001R.fit>: seg pairs xy -> 171, disp_sign_xy=[-209.   13.]
2026-06-23 15:23 - pid 3911 - [astrometry.py:318] - DEBUG - <ReducedFit 77 | OSN-T090/2023-11-06/BLLac_IAR-0001R.fit>: seg pairs xy best -> 89 (42.4%), seg_disp_sign_xy_best=[-209.32879376   14.21187073]
2026-06-23 15:23 - pid 3911 - [astrometry.py:349] - DEBUG - <ReducedFit 77 | OSN-T090/2023-11-06/BLLac_IAR-0001R.fit>: trying to solve astrometry with Seg Best XY Pairs (n=89) (output_logodds_threshold=14).
2026-06-23 15:23 - pid 3943 - [__init__.py:215] - INFO - loaded 240 index files
2026-06-23 15:23 - pid 3943 - [__init__.py:297] - INFO - solve 1: start
2026-06-23 15:23 - pid 3943 - [__init__.py:297] - INFO - solve 1: slice=[0, 25) (1 / 3), index="5200/index-5200-14.fits" (1 / 5)
2026-06-23 15:23 - pid 3943 - [__init__.py:297] - INFO - solve 1: slice=[0, 25) (1 / 3), index="5200/index-5201-14.fits" (2 / 5)
2026-06-23 15:23 - pid 3943 - [__init__.py:297] - INFO - solve 1: logodds=56.9665, matches=29, conflicts=0, distractors=43, ra=330.684, dec=42.2756, scale=0.133194, index="5200/index-5201-14.fits"
2026-06-23 15:23 - pid 3943 - [__init__.py:297] - INFO - solve 1: slice=[0, 25) (1 / 3), index="5200/index-5202-14.fits" (3 / 5)
2026-06-23 15:23 - pid 3943 - [__init__.py:297] - INFO - solve 1: slice=[0, 25) (1 / 3), index="5200/index-5203-14.fits" (4 / 5)
2026-06-23 15:23 - pid 3943 - [__init__.py:297] - INFO - solve 1: slice=[0, 25) (1 / 3), index="5200/index-5204-14.fits" (5 / 5)
2026-06-23 15:23 - pid 3943 - [__init__.py:297] - INFO - solve 1: slice=[25, 50) (2 / 3), index="5200/index-5200-14.fits" (1 / 5)
2026-06-23 15:23 - pid 3943 - [__init__.py:297] - INFO - solve 1: logodds=229.138, matches=57, conflicts=18, distractors=9, ra=330.684, dec=42.2756, scale=0.133993, index="5200/index-5200-14.fits"
2026-06-23 15:23 - pid 3911 - [astrometry.py:364] - DEBUG - <ReducedFit 77 | OSN-T090/2023-11-06/BLLac_IAR-0001R.fit>: Seg Best XY Pairs (n=89) worked.
2026-06-23 15:23 - pid 3911 - [astrometry.py:365] - DEBUG - <ReducedFit 77 | OSN-T090/2023-11-06/BLLac_IAR-0001R.fit>: bm.index_path=PosixPath('/mnt/astrometry_cache/5200/index-5200-14.fits')
2026-06-23 15:23 - pid 3911 - [astrometry.py:366] - DEBUG - <ReducedFit 77 | OSN-T090/2023-11-06/BLLac_IAR-0001R.fit>: bm.center_ra_deg=330.6843284755797
2026-06-23 15:23 - pid 3911 - [astrometry.py:367] - DEBUG - <ReducedFit 77 | OSN-T090/2023-11-06/BLLac_IAR-0001R.fit>: bm.center_dec_deg=42.275589764486774
2026-06-23 15:23 - pid 3911 - [astrometry.py:368] - DEBUG - <ReducedFit 77 | OSN-T090/2023-11-06/BLLac_IAR-0001R.fit>: bm.scale_arcsec_per_pixel=0.13399307634286572
2026-06-23 15:23 - pid 3911 - [astrometry.py:369] - DEBUG - <ReducedFit 77 | OSN-T090/2023-11-06/BLLac_IAR-0001R.fit>: bm.logodds=229.13844299316406
2026-06-23 15:23 - pid 3911 - [astrometry.py:222] - DEBUG - <ReducedFit 77 | OSN-T090/2023-11-06/BLLac_IAR-0001R.fit>: WCS built with attempt 1 / 9 ({'keep_n_seg': 300, 'border_margin_px': 20, 'output_logodds_threshold': 14, 'n_rms_seg': 1.2, 'bkg_filter_size': 11, 'bkg_box_size': 32, 'seg_fwhm': 1.0, 'npixels': 32, 'seg_kernel_size': None, 'position_hint': PositionHint(ra_deg=330.6804125, dec_deg=42.27766361111111, radius_deg=0.23137333333333335), 'size_hint': SizeHint(lower_arcsec_per_pixel=0.1273, upper_arcsec_per_pixel=0.14070000000000002), 'disp_sign_mean': array([-209.,   13.]), 'disp_sign_err': array([12., 16.])}).
2026-06-23 15:23 - pid 3911 - [astrometry.py:234] - DEBUG - <ReducedFit 77 | OSN-T090/2023-11-06/BLLac_IAR-0001R.fit>: building summary images.
2026-06-23 15:23 - pid 3911 - [plotting.py:475] - DEBUG - <ReducedFit 77 | OSN-T090/2023-11-06/BLLac_IAR-0001R.fit>: plotting astrometry summary image of background substraction results
2026-06-23 15:24 - pid 3911 - [plotting.py:487] - DEBUG - <ReducedFit 77 | OSN-T090/2023-11-06/BLLac_IAR-0001R.fit>: plotting astrometry summary image of segmentation results
2026-06-23 15:24 - pid 3911 - [plotting.py:587] - DEBUG - <ReducedFit 77 | OSN-T090/2023-11-06/BLLac_IAR-0001R.fit>: plotting astrometry summary image of astrometry results
2026-06-23 15:24 - pid 3911 - [plotting.py:213] - DEBUG - <ReducedFit 77 | OSN-T090/2023-11-06/BLLac_IAR-0001R.fit>: found 17 catalog sources in field: <AstroSourceQuerySet [AstroSource.objects.get(name='2200+420'), AstroSource.objects.get(name='PanSTARRS 330.70127 42.31286'), AstroSource.objects.get(name='PanSTARRS 330.65424 42.22716'), AstroSource.objects.get(name='PanSTARRS 330.69304 42.28230'), AstroSource.objects.get(name='PanSTARRS 330.60371 42.27227'), AstroSource.objects.get(name='PanSTARRS 330.78031 42.29269'), AstroSource.objects.get(name='PanSTARRS 330.74253 42.30545'), AstroSource.objects.get(name='PanSTARRS 330.68965 42.26675'), AstroSource.objects.get(name='PanSTARRS 330.61553 42.24694'), AstroSource.objects.get(name='PanSTARRS 330.71458 42.25859'), AstroSource.objects.get(name='PanSTARRS 330.65030 42.28163'), AstroSource.objects.get(name='PanSTARRS 330.67245 42.31685'), AstroSource.objects.get(name='PanSTARRS 330.63916 42.24036'), AstroSource.objects.get(name='PanSTARRS 330.58896 42.22607'), AstroSource.objects.get(name='PanSTARRS 330.67147 42.25858'), AstroSource.objects.get(name='PanSTARRS 330.77732 42.24798'), AstroSource.objects.get(name='PanSTARRS 330.73073 42.23001')]>
2026-06-23 15:24 - pid 3911 - [connectionpool.py:1062] - DEBUG - Starting new HTTPS connection (1): simbad.cds.unistra.fr:443
2026-06-23 15:24 - pid 3911 - [connectionpool.py:544] - DEBUG - https://simbad.cds.unistra.fr:443 "GET /simbad/sim-tap/capabilities HTTP/1.1" 200 None
2026-06-23 15:24 - pid 3911 - [connectionpool.py:544] - DEBUG - https://simbad.cds.unistra.fr:443 "POST /simbad/sim-tap/sync HTTP/1.1" 200 None
2026-06-23 15:24 - pid 3911 - [connectionpool.py:544] - DEBUG - https://simbad.cds.unistra.fr:443 "POST /simbad/sim-tap/sync HTTP/1.1" 200 None
2026-06-23 15:24 - pid 3911 - [connectionpool.py:544] - DEBUG - https://simbad.cds.unistra.fr:443 "POST /simbad/sim-tap/sync HTTP/1.1" 200 None
2026-06-23 15:24 - pid 3911 - [logger.py:235] - WARNING - AstropyDeprecationWarning: "epoch" was deprecated in version 0.4.8 and will be removed in a future version. 
2026-06-23 15:24 - pid 3911 - [connectionpool.py:544] - DEBUG - https://simbad.cds.unistra.fr:443 "POST /simbad/sim-tap/sync HTTP/1.1" 200 None
2026-06-23 15:24 - pid 3911 - [instrument.py:485] - DEBUG - <ReducedFit 77 | OSN-T090/2023-11-06/BLLac_IAR-0001R.fit>: saving WCSs to FITS header.
2026-06-23 15:24 - pid 3911 - [instrument.py:557] - DEBUG - <ReducedFit 77 | OSN-T090/2023-11-06/BLLac_IAR-0001R.fit>: astrometric calibration was successful.
2026-06-23 15:24 - pid 3911 - [instrument.py:560] - DEBUG - <ReducedFit 77 | OSN-T090/2023-11-06/BLLac_IAR-0001R.fit>: searching for sources in field...
2026-06-23 15:24 - pid 3911 - [instrument.py:563] - DEBUG - <ReducedFit 77 | OSN-T090/2023-11-06/BLLac_IAR-0001R.fit>: found 17 sources in field.

WARNING: AstropyDeprecationWarning: "epoch" was deprecated in version 0.4.8 and will be removed in a future version.  [iop4lib.utils]

As we can see, the astrometric calibration was done with the blind astrometry.net solver as before.

Next, if we go to the polarimetry field, we see that there are not enough stars in the image to do the blind astrometric calibration.

The calibration of this file will use the previously calibrated photometry field and compare quads of stars in both images to find the transformation between them, as we can see in the log:

redf_pol = ReducedFit.create(rf_pol)
redf_pol.build_file()

2026-06-23 15:24 - pid 3911 - [reducedfit.py:123] - DEBUG - DB entry of ReducedFit for OSN-T090/2023-11-06/BLLAC_R_IAR-0760.fts already exists, it will be used instead.
2026-06-23 15:24 - pid 3911 - [dipol.py:194] - WARNING - I could not find a master with same size, searching for any size
2026-06-23 15:24 - pid 3911 - [dipol.py:205] - WARNING - Checking for flats of any imgsize
2026-06-23 15:24 - pid 3911 - [dipol.py:194] - WARNING - I could not find a master with same size, searching for any size
2026-06-23 15:24 - pid 3911 - [dipol.py:205] - WARNING - Checking for flats of any imgsize
2026-06-23 15:24 - pid 3911 - [dipol.py:194] - WARNING - I could not find a master with same size, searching for any size
2026-06-23 15:24 - pid 3911 - [dipol.py:205] - WARNING - Checking for flats of any imgsize
2026-06-23 15:24 - pid 3911 - [instrument.py:433] - WARNING - <ReducedFit 79 | OSN-T090/2023-11-06/BLLAC_R_IAR-0760.fts>: masterflat in this epoch could not be found, attemptying adjacent epochs.
2026-06-23 15:24 - pid 3911 - [instrument.py:404] - DEBUG - No Master Flat for {'instrument': 'DIPOL', 'obsmode': 'POLARIMETRY', 'imgsize': '1104x896', 'imgbinning': '2x2', 'band': 'R', 'rotangle': 157.5} in DB, None will be returned.
2026-06-23 15:24 - pid 3911 - [dipol.py:194] - WARNING - I could not find a master with same size, searching for any size
2026-06-23 15:24 - pid 3911 - [dipol.py:205] - WARNING - Checking for flats of any imgsize
2026-06-23 15:24 - pid 3911 - [instrument.py:540] - DEBUG - <ReducedFit 79 | OSN-T090/2023-11-06/BLLAC_R_IAR-0760.fts>: building file
2026-06-23 15:24 - pid 3911 - [dipol.py:248] - DEBUG - <ReducedFit 79 | OSN-T090/2023-11-06/BLLAC_R_IAR-0760.fts>: applying masters
2026-06-23 15:24 - pid 3911 - [instrument.py:546] - DEBUG - <ReducedFit 79 | OSN-T090/2023-11-06/BLLAC_R_IAR-0760.fts>: performing astrometric calibration
2026-06-23 15:24 - pid 3911 - [logger.py:235] - WARNING - AstropyDeprecationWarning: "epoch" was deprecated in version 0.4.8 and will be removed in a future version. 
2026-06-23 15:24 - pid 3911 - [connectionpool.py:544] - DEBUG - https://simbad.cds.unistra.fr:443 "POST /simbad/sim-tap/sync HTTP/1.1" 200 None
2026-06-23 15:24 - pid 3911 - [dipol.py:491] - DEBUG - target_src.srctype='blazar'
2026-06-23 15:24 - pid 3911 - [dipol.py:492] - DEBUG - n_estimate=11
2026-06-23 15:24 - pid 3911 - [dipol.py:493] - DEBUG - n_estimate_centered=5
2026-06-23 15:24 - pid 3911 - [dipol.py:494] - DEBUG - redf_phot=ReducedFit.objects.get(id=77)
2026-06-23 15:24 - pid 3911 - [dipol.py:495] - DEBUG - n_expected_simbad_sources=0
2026-06-23 15:24 - pid 3911 - [dipol.py:496] - DEBUG - n_expected_calibrators=30
2026-06-23 15:24 - pid 3911 - [dipol.py:550] - DEBUG - Trying attempt: {'method': '_build_wcs_for_polarimetry_images_photo_quads', 'conds': {'redf_phot__ne': None}, 'args': {'n_seg_threshold': [1.1, 1.0, 0.9], 'npixels': [64, 32], 'min_quad_distance': [4.0, 8.0]}} for <ReducedFit 79 | OSN-T090/2023-11-06/BLLAC_R_IAR-0760.fts>.
2026-06-23 15:24 - pid 3911 - [dipol.py:562] - INFO - Attempt: {'method': '_build_wcs_for_polarimetry_images_photo_quads', 'conds': {'redf_phot__ne': None}, 'args': {'n_seg_threshold': [1.1, 1.0, 0.9], 'npixels': [64, 32], 'min_quad_distance': [4.0, 8.0]}}: _build_wcs_for_polarimetry_images_photo_quads with {'n_seg_threshold': 1.1, 'npixels': 64, 'min_quad_distance': 4.0} for <ReducedFit 79 | OSN-T090/2023-11-06/BLLAC_R_IAR-0760.fts>.
2026-06-23 15:24 - pid 3911 - [dipol.py:703] - DEBUG - Invoked with n_seg_threshold=1.1, npixels=64
2026-06-23 15:24 - pid 3911 - [dipol.py:727] - DEBUG - Using 10 sources in polarimetry field and 10 in photometry field.
2026-06-23 15:24 - pid 3911 - [dipol.py:730] - DEBUG - Building summary image for astrometry detected sources.
2026-06-23 15:24 - pid 3911 - [dipol.py:819] - DEBUG - Selected 5 quads with distance < 4.0. I will get the one with less deviation from the median linear transform.
2026-06-23 15:24 - pid 3911 - [dipol.py:827] - DEBUG - Filtering out big translations (<1400 px)
2026-06-23 15:24 - pid 3911 - [dipol.py:831] - DEBUG - Filtering large transformations
2026-06-23 15:24 - pid 3911 - [dipol.py:839] - DEBUG - Filtered to 5 quads with distance < 4.0 and translation < 1400 px.
2026-06-23 15:24 - pid 3911 - [dipol.py:880] - DEBUG - median_t=array([100.24747473,  27.39696465])
2026-06-23 15:24 - pid 3911 - [dipol.py:886] - DEBUG - Selected the quads [7,1]
2026-06-23 15:24 - pid 3911 - [dipol.py:890] - DEBUG - t = [100.18744966  27.39696465], R = [[ 9.99999897e-01 -4.53157715e-04]
 [ 4.53157715e-04  9.99999897e-01]]
2026-06-23 15:24 - pid 3911 - [dipol.py:891] - DEBUG - det R = 1.0000000000000002
2026-06-23 15:24 - pid 3911 - [dipol.py:895] - DEBUG - Building summary image for quad matching.
2026-06-23 15:24 - pid 3911 - [__init__.py:181] - WARNING - Less than 5 calibrated fits for DIPOL 2200+420, using all of them
2026-06-23 15:24 - pid 3911 - [dipol.py:938] - DEBUG - Using angle=181.3358155263559 for pre wcs.
2026-06-23 15:24 - pid 3911 - [dipol.py:960] - DEBUG - Building summary image for astrometry.
2026-06-23 15:24 - pid 3911 - [plotting.py:213] - DEBUG - <ReducedFit 79 | OSN-T090/2023-11-06/BLLAC_R_IAR-0760.fts>: found 4 catalog sources in field: <AstroSourceQuerySet [AstroSource.objects.get(name='2200+420'), AstroSource.objects.get(name='PanSTARRS 330.69304 42.28230'), AstroSource.objects.get(name='PanSTARRS 330.68965 42.26675'), AstroSource.objects.get(name='PanSTARRS 330.67147 42.25858')]>
2026-06-23 15:24 - pid 3911 - [logger.py:235] - WARNING - AstropyDeprecationWarning: "epoch" was deprecated in version 0.4.8 and will be removed in a future version. 
2026-06-23 15:24 - pid 3911 - [connectionpool.py:544] - DEBUG - https://simbad.cds.unistra.fr:443 "POST /simbad/sim-tap/sync HTTP/1.1" 200 None
2026-06-23 15:24 - pid 3911 - [instrument.py:485] - DEBUG - <ReducedFit 79 | OSN-T090/2023-11-06/BLLAC_R_IAR-0760.fts>: saving WCSs to FITS header.
2026-06-23 15:24 - pid 3911 - [instrument.py:557] - DEBUG - <ReducedFit 79 | OSN-T090/2023-11-06/BLLAC_R_IAR-0760.fts>: astrometric calibration was successful.
2026-06-23 15:24 - pid 3911 - [instrument.py:560] - DEBUG - <ReducedFit 79 | OSN-T090/2023-11-06/BLLAC_R_IAR-0760.fts>: searching for sources in field...
2026-06-23 15:24 - pid 3911 - [instrument.py:563] - DEBUG - <ReducedFit 79 | OSN-T090/2023-11-06/BLLAC_R_IAR-0760.fts>: found 4 sources in field.

WARNING: AstropyDeprecationWarning: "epoch" was deprecated in version 0.4.8 and will be removed in a future version.  [iop4lib.utils]
WARNING: AstropyDeprecationWarning: "epoch" was deprecated in version 0.4.8 and will be removed in a future version.  [iop4lib.utils]

!ls {redf_pol.filedpropdir}

from IPython.display import Image
Image(filename=f"{redf_pol.filedpropdir}/astrometry_matched_quads.png")

BLLAC_R_IAR-0760.fts		 astrometry_matched_quads.png
astrometry_detected_sources.png  astrometry_summary.png
astrometry_info

The image shows the quads of stars that were matched in both images, their distance in the hash space and the distance in the image space.

This quad matching is necessary because the number of stars in the polarimetry field is not enough for the astrometry.net solver to work, plus the scale of the images is too small for the default index files. This quad matching needs at least 4 (ideally a few more) stars in the image, and a previously calibrated photometry field. If any of these conditions are not met, other methods must be used to calibrate the image (if there are only two bright sources at the right distance, can be assumed that they are the target star, and if there are between 3 and 4-6 stars, we can try to match to known sources in the field).

The process is implemented in the DIPOL._build_wcs_for_polarimetry_images_photo_quads method, which is called by the ReducedFit.astrometric_calibration method. We could also manually follow the process:

import numpy as np
import itertools
from iop4lib.instruments import DIPOL
from iop4lib.utils.astrometry import BuildWCSResult
from photutils.aperture import CircularAperture
from pathlib import Path
import matplotlib as mplt
import matplotlib.pyplot as plt

# get the subframe of the photometry field that corresponds to this polarimetry field, (approx)
x_start = redf_pol.rawfit.header['XORGSUBF']
y_start = redf_pol.rawfit.header['YORGSUBF']

x_end = x_start + redf_pol.rawfit.header['NAXIS1']
y_end = y_start + redf_pol.rawfit.header['NAXIS2']

idx = np.s_[y_start:y_end, x_start:x_end]

photdata_subframe = redf_phot.mdata[idx]

# find 10 brightest sources in each field

sets_L = list()

for redf, data in zip([redf_pol, redf_phot], [redf_pol.mdata, photdata_subframe]):

    positions = DIPOL._estimate_positions_from_segments(redf=redf, data=data, n_seg_threshold=1.5, npixels=32, centering=None, fwhm=1.0)
    positions = positions[:10]

    sets_L.append(positions)

fig = plt.figure(figsize=(12,6), dpi=iop4conf.mplt_default_dpi)
axs = fig.subplots(nrows=1, ncols=2)

for ax, data, positions in zip(axs, [redf_pol.mdata, photdata_subframe], sets_L):
    imshow_w_sources(data, pos1=positions, ax=ax)
    candidates_aps = CircularAperture(positions[:2], r=10.0)
    candidates_aps.plot(ax, color="b")
    for i, (x,y) in enumerate(positions):
        ax.text(x, y, f"{i}", color="orange", fontdict={"size":14, "weight":"bold"})#, verticalalignment="center", horizontalalignment="center") 
    ax.plot([data.shape[1]//2], [data.shape[0]//2], '+', color='y', markersize=10)
    
axs[0].set_title("Polarimetry field")
axs[1].set_title("Photometry field")
fig.suptitle("astrometry detected_sources")
plt.show()

# Build the quads for each field
quads_1 = np.array(list(itertools.combinations(sets_L[0], 4)))
quads_2 = np.array(list(itertools.combinations(sets_L[1], 4)))

from iop4lib.utils.quadmatching import hash_ish, distance, order, qorder_ish, find_linear_transformation
hash_func, qorder = hash_ish, qorder_ish

hashes_1 = np.array([hash_func(quad) for quad in quads_1])
hashes_2 = np.array([hash_func(quad) for quad in quads_2])

all_indices = np.array(list(itertools.product(range(len(quads_1)),range(len(quads_2)))))
all_distances = np.array([distance(hashes_1[i], hashes_2[j]) for i,j in all_indices])

idx = np.argsort(all_distances)
all_indices = all_indices[idx]
all_distances = all_distances[idx]

indices_selected = all_indices[:6]
colors = ["r", "g", "b", "c", "m", "y"]

figs, axs = zip(*[plt.subplots(figsize=(6,4), dpi=iop4conf.mplt_default_dpi) for _ in range(2)])

for (i,j), color in list(zip(indices_selected, colors)): 
    
    tij = find_linear_transformation(qorder(quads_1[i]), qorder(quads_2[j]))[1]

    for ax, fig, data, quad, positions in zip(axs, figs, [redf_pol.mdata, photdata_subframe], [quads_1[i], quads_2[j]], sets_L):
        imshow_w_sources(data, pos1=positions, ax=ax)
        x, y = np.array(order(quad)).T
        ax.fill(x, y, edgecolor='k', fill=True, facecolor=mplt.colors.to_rgba(color, alpha=0.2))

plt.show()