Incorrect co-registration

[kanishk-aidash]

• arosics version: 1.3.0
• Python version: 3.9.2
• Operating System: MacOD Big Sur Version 11.4

Description

Co-registered two 4-band rasters (RGB-NIR), but the detected absolute shift vector length is not correct. Although, the flow is working fine but the final output isn't much corrected / co-registered w.r.t. the reference raster. QGIS tool is able to co-register the rasters, but requires a lot of tie-points to be selected

What I Did

im_reference = 'reference_4band.tiff'
im_target = 'target_4band.tiff'

{
            'grid_res': 64,
            'window_size': (32, 32),
            'min_reliability': 10,
            'tieP_filter_level': 3,
            'q': False,
            'v': True
            'max_shift': 128,
            'max_iter': 20,
            'r_b4match': 4,
            's_b4match': 4,
            'path_out': '/tmp/coregistration/output.tif',
            'fmt_out': 'GTIFF'
        }

CRL = COREG_LOCAL(im_reference, im_target, **kwargs)
CRL.correct_shifts()
CRL.view_CoRegPoints(figsize=(15, 15), backgroundIm='ref', savefigPath='/tmp/coregistration/tiepoints.png')

[danschef]

Hi, please first make sure you are using the latest version of AROSICS. Version 1.3.0 is quite outdated. However, it is possible that a package incompatibility prevented your package manager from installing the latest one. But this should be resolved in the latest version of AROSICS now (1.4.6).

I suggest to re-create the environment from scratch as this is safest method to get all the latest package versions. Note that using mamba (run conda install -c conda-forge mamba in your conda base environment) instead of conda will speed up the installation process a lot. Then just run mamba install -c conda-forge arosics or mamba create -n arosics -c conda-forge arosics.

[kanishk-aidash]

Hey @danschef ,

I re-ran the above code with the new AROSICS (1.4.6), but results are looking all the same!

I have experimented / fiddled with the kwargs parameters quite a bit, and this configuration gives best number of tie-points for my input rasters. In different configurations, only 4-5 tie-points remain and the target raster is left as it is, without any co-registration

[danschef]

Ok, could you please post the log output of AROSICS here? Also a figure which comes out of CRL.view_CoRegPoints(figsize=(15, 15), backgroundIm='tgt', hide_filtered=True) would be helpful.

[kanishk-aidash]

AROSICS run output Logs:

Coreg_points: Reference

Coreg_points: Target

[danschef]

Ok, what directly comes into my mind when looking at your input parameters is the very small size of the matching window. It seems like your input images have a quite high spatial resolution and you expect shifts up to 128 (which is really a lot!!). However, when using a matching window of 32x32 pixels, a shift in that dimension could never be detected as the corresponding positions in the reference and the target window would not be covered both.

So I suggest to use a matching window size of 256x256 pixels at least (which is the default). 512x512 could even work better. Also don´t set the reliability threshold to 10%, this would allow very unreliable matches to still be accepted. Try it with the default value and maybe decrease it a bit in case there are too many points filtered out by the reliability filtering. You can also see this by looking at the output of CRL.view_CoRegPoints(figsize=(15, 15), hide_filtered=False). Setting the maximum number of iterations to 20 should not improve the result much but just increase the computational effort.

[kanishk-aidash]

@danschef

Yes, our data is 50cm client data, and have ~5m shifts. I set these parameters after quite a bit of trial and error. Selecting a bigger matching window size significantly reduces the detected tie-points (from > 400 to < 50, and only 3-4 remain after 3 level of filtering.

Also, seems like I am still confused about the 'gridres' parameter. It says 'tie point grid resolution in pixels of the target image (x-direction)'_, not entirely sure what to make of it.

[danschef]

This looks like is was not possible to find a match at most of the points. You can check if there are any special reasons for this by taking a look at the "LAST_ERR" column of CRL.CoRegPoints_table which contains the error messages for all points of the tie point grid. Also looking at an exemplary matching window may help to get an idea what is the problem:

from arosics import COREG
CR = COREG(im_reference, im_target)
CR.calculate_spatial_shifts()
CR.show_matchWin(interactive=True)

I could imagine that the problem is related to the acquisition and illumination geometries of your input images. To me, the images look like they are acquired at quite different day times which means that you have a quite different length of shadows in your image. In the frequency domain (where AROSICS detects the coregistration point), this will reduce the similarity of frequency features which can then cause the tie point detection to fail. This especially applies to urban surfaces where you have strong effects of shadow lengths and acquisition directions on your images. AROSICS is mainly developed for satellite data with spatial resolutions of, lets say, 5 meters or more. There, the effect of shadows is much lower compared to high resolution image data. So, to make it work in your case, you could try to replace the reference or the target image by one with a more similar illumination or acquisition geometry. Note, that one of the requirements of AROSICS is to use orthorectified images (as stated in the documentation).

If that does not help, I could offer to take a look at it by myself, if you send me the input data.

[kanishk-aidash]

Sorry, I can't share the data with you, it's proprietary.

I am not entirely sure if that's what's happening in my case. Because if I change the parameters, I am getting a good number of matches spread across the raster.

When it gets these matches the absolute vector shift value is not coming out right, although the tie point is correct.

Is there a way, I can add manually detected Tie-Points (using feature matching) to the algorithm in a quality controlled manner?

[kanishk-aidash]

@danschef Any update?

[danschef]

When it gets these matches the absolute vector shift value is not coming out right, although the tie point is correct.

Could you explain that a bit more? Does this mean you checked an exemplary tie point with the code snippet from https://github.com/GFZ/arosics/issues/5#issuecomment-868380803 and the misregistration was correctly removed but the computed length of the shift vector (as displayed in the log output) was wrong?

Is there a way, I can add manually detected Tie-Points (using feature matching) to the algorithm in a quality controlled manner?

So far, there is not directly a function or method to provide custom tie-points on your own. However, you could still just add them to CRL.CoRegPoints_table before calling CRL.correct_shifts().

[kanishk-aidash]

@danschef Yes! .. The misregistration is getting corrected to some extent not fully! I am doing Local coregistration and the code snippet doesn't work with COREG_LOCAL. I have to cross-verify the correction via QGIS.

Adding tie-points this way doesn't compute the shift vector correctly. I have tried adding tie-points detected via keypoints matching to the table, but the shift computed isn't correct/aligned with the one computed via the arosics algorithm

[danschef]

the code snippet doesn't work with COREG_LOCAL

Right, the intention of the snippet in https://github.com/GFZ/arosics/issues/5#issuecomment-868380803 was to find out what goes wrong in your case. Therefore, it makes sense to reproduce the issue with the most simple co-registration case, i.e., for a single tie point / global co-registration. As soon as that works properly, you can go on with local co-registation. I think it would be helpful to take a look at the output of CR.show_matchWin(interactive=True). You have to run this is in a Jupyter notebook to get an interactive figure. Otherwise, set interactive to False. If possible, post a screenshot of the figure here. Otherwise, it is a bit hard for me to find out what goes wrong.

Also make sure to use the latest version of AROSICS. There were some changes since 1.4.6 that might also affect this issue.

[kanishk-aidash]

Ran the snippet in #5

The interactive window doesn't show much improvement:

Have tried running local co-registration on a small raster where there is a good amount of shift, but the the module doesn't seem to find any good tie-points or correct the shift in local co-registration.

Attaching a sample data:

data.zip

Yes, I am on arosics 1.5.0 but don't see much differenct

[danschef]

Thanks for the sample data, that makes it a lot easier. When I pass your data to the global coregistration, the first thing I get is a RuntimeError that tells you that the projections of your input data are unequal. This is currently not directly supported which is documented here. The reference image is in WGS84/Pseudo Mercator projection and the target image has Longitude/Latitude coordinates (EPSG 4626):

p_ref = '/path/reference.tif'
p_tgt = '/path/target.tif'

CR = COREG(ref, tgt, max_shift=20)
CR.calculate_spatial_shifts()

------------
~/scheffler/python/geoarray/geoarray/baseclasses.py:764: UserWarning: 
Automatic nodata value detection returned the value 0.0 for GeoArray 'target' but this seems to be unreliable (occurs in only 1 image corner). To avoid automatic detection, just pass the correct nodata value.
  warnings.warn("\nAutomatic nodata value detection returned the value %s for GeoArray '%s' but this "
Polygonize progress     |==================================================| 100.0% Complete  => 0:00:00
Calculating footprint polygon and actual data corner coordinates for reference image...
Bounding box of calculated footprint for reference image:
    (-11721421.8588, 4355165.4727, -11720851.8618, 4355461.9149)
Calculating footprint polygon and actual data corner coordinates for image to be shifted...
Automatically detected nodata value for GeoArray_CoReg 'target': 0.0
Bounding box of calculated footprint for image to be shifted:
    (-105.295324074, 36.395402778, -105.290203704, 36.397546296)
---------------------------------------------------------------------------
RuntimeError                              Traceback (most recent call last)
/tmp/ipykernel_1118122/2090732032.py in <module>
----> 1 CR = COREG(ref, tgt, max_shift=20, path_out='/path/coreg_global.tif')
      2 CR.calculate_spatial_shifts()

~/scheffler/python/arosics/arosics/CoReg.py in __init__(self, im_ref, im_tgt, path_out, fmt_out, out_crea_options, r_b4match, s_b4match, wp, ws, max_iter, max_shift, align_grids, match_gsd, out_gsd, target_xyGrid, resamp_alg_deshift, resamp_alg_calc, footprint_poly_ref, footprint_poly_tgt, data_corners_ref, data_corners_tgt, nodata, calc_corners, binary_ws, mask_baddata_ref, mask_baddata_tgt, CPUs, force_quadratic_win, progress, v, path_verbose_out, q, ignore_errors)
    409         gdal.AllRegister()
    410         self._check_and_handle_metaRotation()
--> 411         self._get_image_params()
    412         self._set_outpathes(im_ref, im_tgt)
    413         self.grid2use = 'ref' if self.shift.xgsd <= self.ref.xgsd else 'shift'

~/scheffler/python/arosics/arosics/CoReg.py in _get_image_params(self)
    555                 (crs_ref.name, crs_shift.name) if not crs_ref.name == crs_shift.name else (crs_ref.srs, crs_shift.srs)
    556 
--> 557             raise RuntimeError(
    558                 'Input projections are not equal. Different projections are currently not supported. '
    559                 'Got %s vs. %s.' % (name_ref, name_shift))

RuntimeError: Input projections are not equal. Different projections are currently not supported. Got WGS 84 / Pseudo-Mercator vs. WGS 84.

After reprojecting your reference image to EPSG 4326, everything works fine:

from geoarray  import GeoArray

ref = GeoArray(p_ref)
tgt = GeoArray(p_tgt)

ref.reproject_to_new_grid(prototype=tgt)

CR = COREG(ref, tgt, max_shift=20, path_out='/path/coreg_global.tif')
CR.calculate_spatial_shifts()

-------
~/scheffler/python/geoarray/geoarray/baseclasses.py:764: UserWarning: 
Automatic nodata value detection returned the value 0.0 for GeoArray 'target' but this seems to be unreliable (occurs in only 1 image corner). To avoid automatic detection, just pass the correct nodata value.
  warnings.warn("\nAutomatic nodata value detection returned the value %s for GeoArray '%s' but this "
Polygonize progress     |==================================================| 100.0% Complete  => 0:00:00
Calculating footprint polygon and actual data corner coordinates for reference image...
Bounding box of calculated footprint for reference image:
    (-105.29531944437072, 36.395407407630074, -105.29019907437826, 36.397546296)
Calculating footprint polygon and actual data corner coordinates for image to be shifted...
Automatically detected nodata value for GeoArray_CoReg 'target': 0.0
Bounding box of calculated footprint for image to be shifted:
    (-105.295324074, 36.395402778, -105.290203704, 36.397546296)
Matching window position (X,Y): -105.29327269469178/36.39671560139622
Detected integer shifts (X/Y):                            2/9
Detected subpixel shifts (X/Y):                           -0.2697507613216568/0.4160501972531371
Calculated total shifts in fft pixel units (X/Y):         1.7302492386783432/9.416050197253137
Calculated total shifts in reference pixel units (X/Y):   1.7302492386783432/9.416050197253137
Calculated total shifts in target pixel units (X/Y):      1.7302492386783432/9.416050197253137
Calculated map shifts (X,Y):                  8.010412557268864e-06/-4.3592814442661165e-05
Calculated absolute shift vector length in map units:     4.432268245909693e-05
Calculated angle of shift vector in degrees from North:   349.58775256870763
Original map info: ['Geographic Lat/Lon', 1.0, 1.0, -105.295324074, 36.397546296, 4.6296292947560644e-06, 4.629628509727879e-06, 'WGS-84']
Updated map info:  ['Geographic Lat/Lon', 1.0, 1.0, '-105.29531606358745', '36.39750270318556', 4.6296292947560644e-06, 4.629628509727879e-06, 'WGS-84']
Image similarity within the matching window (SSIM before/after correction): 0.1720 => 0.3572
Estimated reliability of the calculated shifts:  72.4 %
'success'

When looking at the co-registered matching window, I can´t see an issue (you have to move the slider on right): CR.show_matchWin(interactive=True)

(top: before, bottom: after co-registration, note that the co-registration point is in the center of the image)

Also the local co-registration detects a lot of tie points:

CRL = COREG_LOCAL(ref, tgt, 25, max_shift=20, nodata=(0, 0))
CRL.calculate_spatial_shifts()
CRL.view_CoRegPoints(hide_filtered=False)

-------------------
Polygonize progress     |==================================================| 100.0% Complete  => 0:00:00
Polygonize progress     |==========================------------------------| 51.0% Complete  => 0:00:00
Calculating footprint polygon and actual data corner coordinates for reference image...
Bounding box of calculated footprint for reference image:
    (-105.29531944437072, 36.395407407630074, -105.29020370400755, 36.397546296)
Calculating footprint polygon and actual data corner coordinates for image to be shifted...
Polygonize progress     |==================================================| 100.0% Complete  => 0:00:00
Bounding box of calculated footprint for image to be shifted:
    (-105.295324074, 36.395402778, -105.290203704, 36.397546296)
Matching window position (X,Y): -105.29327269469354/36.396715601395975
Initializing tie points grid...
Calculating tie point grid (582 points) using 32 CPU cores...
    progress: |==================================================| 100.0% Complete  => 0:00:17
Found 236 matches.
Performing validity checks...
118 tie points flagged by level 1 filtering (reliability).
2 tie points flagged by level 2 filtering (SSIM).
0 tie points flagged by level 3 filtering (RANSAC)
118 valid tie points remain after filtering.

Here are the detected tie points with grid_res=25, including false-positives flagged by the three outlier detection algorithms:

Note that the area without tie points at the bottom is likely related to the relatively large shift and the automatic reduction of the matching window size due to the nodata area. You may reduce this tie-point-less area by first running a global co-registration with align_grids=False (to avoid resampling) and then running a local co-registration to fine-tune the result.

Attached are the co-registration results as produced by the code above. data.zip

[kanishk-aidash]

Hey @danschef Thanks for the detailed response. Sorry for the wrong data, I am re-projecting both to EPSG:4326 at my end, but forgot to attach the re-projected rasters.

The configuration you just shared with me, I have already tried this set of params. I used to get this error earlier as well, which I am still getting after trying out your code snippet:

I get the same error for COREG_LOCAL as well.

This error, seemingly tells that widow size is too large for our image. But, if I try smaller window sizes .. Number of tie-points I get goes down to l.t. 5

Could it have something to do with the OS? I am running the code in a jupyter notebook on macOS 11.4, Python 3.9

I am not sure at this point, what I am missing. I have tried tons of different configs for both COREG and COREG_LOCAL, and have gone through most of the the code of this package! I am sure it's something really stupid that I am unable to pin-point

[danschef]

Well, this just an issue with this single tie point where you try to find a co-registration via the global co-registration approach. In general, it is not possible to find a valid match at each point within the input image overlap. This may have several reasons, e.g., that the center position of the matching window is too close to an area of no-data values. In this case, AROSICS will automatically reduce the size of the matching window to avoid disturbing features in the frequency domain caused by the image edge. At some point, a further reduction makes no sense anymore and the tie point detection fails for this position. Just take a look at CR.show_matchWin() for this specific point. I guess you will see that the corresponding matching window is directly at the image edge. Use another window position to find a valid match then.

Note that the accurracy and reliability of the found tie points depends on the matching window size. You will get a good reliability when using window sizes of 256x256 or 512x512 pixels but the reliability decreases when using smaller windows. This is also the reason why all of the tie points are filtered out by the reliability filtering in your screenshot above. With a window size of 32x32, you get quite bad reliabilities which are all below the reliability threshold used in the filtering. So my advise is to keep the window size at 256x256 in your case which may result in some missing tie points at the bottom of the image overlap but produces reasonable tie points for the rest of the image.

[kanishk-aidash]

@danschef
I guess I missed out on something.. I am trying to run the exact same code snippet you shared in your previous detailed comment on the same data (grid_res=25, window_size=(256, 256) (default), max_shift=20, nodata=(0, 0). In your run you seem to get very good tie-points but for me .. this configuration errors out in Global as well as Local Co-rgistration approach.

Yes, I understand that reducing the window size will reduce the reliability. It becomes apparent in the co-registered raster as it becomes distorted if I don't remove these unreliable points.

My concerns here are:

1. As mentioned in the description. The value of absolute shift vector isn't correct. Some places it just co-registers partially.
2. Why are runs in your setup v/s mine giving such different results. The error screen shot is for running on the same data with same parameters as you have used in your comment. Your Global Run Succeeds, whereas mine fails.

I have also noticed, that passing the file path(matched CRS) v/s passing the geoarray directly also changes the behaviour of the runs

[danschef]

The value of absolute shift vector isn't correct. Some places it just co-registers partially.

In this case, please provide me a location (window position) where a shift is successfully computed but is not correct. If I can reproduce that here, I may have a closer look.

Why are runs in your setup v/s mine giving such different results.

If you compare the matching window position, you will see that the global co-registration runs at a different location than in my case. Mine is running at Matching window position (X,Y): -105.29327269469178/36.39671560139622, yours at -105.2933.../36.3965.... By default, the matching window position is placed at the center of the image overlap polygon. I guess, in your case, the computation of the target image footprint polygon fails and therefore the full target image extent is used. The center point of it is then directly at the image edge which causes the tie point correction to fail. For computing the image footprint, the packages geoarray and py-tools-ds are involved. Please make sure that they are up-to-date. I have geoarray 0.14.2 and py-tools-ds 0.17.7 installed.

I have also noticed, that passing the file path(matched CRS) v/s passing the geoarray directly also changes the behaviour of the runs.

There must be a difference in your matched CRS version and the one you create with ref.reproject_to_new_grid(prototype=tgt).

[danschef]

If you have old versions of the dependencies of arosics installed, I recommend to just re-create your environment with mamba create -c conda-forge --name arosics_env 'arosics==1.5.0'. You might have to install mamba with conda install -n base -c conda-forge mamba beforehand.

[danschef]

Closing this due to inactivity.