GFZ / arosics

AROSICS - Automated and Robust Open-Source Image Co-Registration Software
https://git.gfz-potsdam.de/danschef/arosics
Apache License 2.0
146 stars 26 forks source link

Tie Point Grid Calculation takes ages #36

Closed eikeschuett closed 7 months ago

eikeschuett commented 8 months ago

Description

I am currently trying to perform a local co-registration of an hyperspectral image (410x410 pixels) to a PleiadesNeo scene (subsetted to ~1000x1000 pixels). scene. However, the calculation of the tie point grid takes much longer than I would expect. The last printed message says it is using all 24 CPU cores, but my total CPU usage is ~1%. It is also strange that I still can't see the progres bar of the tie point grid calculation, although the code is running for more than 50 minutes.

I already tried different grid_res, but the result was always the same...

What I Did

from geoarray import GeoArray
import rioxarray
import numpy as np
from arosics import COREG_LOCAL

with rioxarray.open_rasterio("session_000_001k_185_sprad_clipped_ORTHO.tif") as xds:
    xds = xds.fillna(-999)

    data = xds.sel(band = 50).values

    im_target = GeoArray(
        data,
        geotransform = xds.rio.transform().to_gdal(),
        projection = xds.rio.crs.to_wkt(),
        nodata = -999)

with rioxarray.open_rasterio("PNEO3_2023_09_04_10_33_30_L2W_rhow_563_ortho.tif") as im_reference_xds:
    im_reference_xds = im_reference_xds.sel(
                                    x = slice(xds.x.min()-200,xds.x.min()+200),
                                    y = slice(xds.y.max()+200,xds.y.min()-200))

    im_reference_xds = im_reference_xds.fillna(-999)

    im_reference = GeoArray(
        im_reference_xds.data[0,:,:],
        geotransform = im_reference_xds.rio.transform().to_gdal(),
        projection = xds.rio.crs.to_wkt(),
        nodata = -999)

    kwargs = {
        "grid_res": 50,
        "window_size": (256, 256),
        "path_out": None,
        "q": False,
    }

    CRL = COREG_LOCAL(im_reference, im_target, **kwargs)

    CRL.correct_shifts()

And in the console:

Polygonize progress     |==================================================| 100.0% Complete  => 0:00:00
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:
    (587518.874, 6032418.557, 587918.4739999999, 6033077.357)
Calculating footprint polygon and actual data corner coordinates for image to be shifted...
Bounding box of calculated footprint for image to be shifted:
    (587718.1291309658, 6032405.083383153, 588192.1291309658, 6032878.083383153)
Matching window position (X,Y): 587839.3055236468/6032649.451760621
Target window size (256, 256) not possible due to too small overlap area or window position too close to an image edge. New matching window size: (108, 109).
Initializing tie points grid...
Calculating tie point grid (28 points) using 24 CPU cores...
danschef commented 8 months ago

Computing 28 tie points on 24 CPU cores should be a matter of seconds. This looks rather like the code hangs somewhere in multiprocessing for some reason. To figure out what happened exactly, I would need the input data and the exact way how you installed arosics. As a workaround for now, just disable multiprocessing (COREG_LOCAL(..., CPUs=1)) and see what happens.

Also the rest of the above code seems a bit complicated for what it actually does. I would rather go for something like this:

from arosics import COREG_LOCAL

COREG_LOCAL(
    "session_000_001k_185_sprad_clipped_ORTHO.tif",
    "PNEO3_2023_09_04_10_33_30_L2W_rhow_563_ortho.tif",
    grid_res=50,
    window_size=(256, 256),
    r_b4match=50, 
    path_out=None,
    CPUs=1,
    q=False
).correct_shifts()

No-data values are auto-detected, otherwise you can also set them with the nodata parameter. Same applies to the spatial overlap of reference and target image.

danschef commented 8 months ago

Probably, the code hangs because of this: https://git.gfz-potsdam.de/danschef/arosics/-/issues/93. But that is just a guess and I need to have a closer look. However, to me it looks like it is a Windows-specific issue.

eikeschuett commented 8 months ago

Setting the CPUs=1 does the trick! Without multiprocessing everything works fine.

I installed arosics from conda-forge into a fresh anaconda environment and later installed rioxarray, also from conda-forge.

And yes I know, the code is a bit complicated at the moment,. I wanted to have the option to enhance the image quality of my hyperspectral images on the fly because this data contains few objects with little contrast. I may need to come back to you with this problem at some point in the future.

danschef commented 8 months ago

You could let arosics ignore these areas by setting mask_baddata_ref or mask_baddata_tgt, see here.

danschef commented 7 months ago

AROSICS v1.10.2 replaces the former multiprocessing implementation with joblib-based parallelization which should (I think) fix your issue. Would be nice if you could try it out and provide me some feedback.

eikeschuett commented 7 months ago

Works fine for me. All cores are used at 100 % and the computation time reduces significantly, in my particular case ~ factor 5. Thank you very much!

danschef commented 7 months ago

Cool, thanks for trying!

leikas227 commented 6 months ago

I have no problem with co-registration using the Global approach. However, the error message below occurs when trying the local registration, even when I use only 1 CPU. Different grid_res same error. Any suggestions? arosics version: 1.10.2 Python version: 3.11.8 Operating System: Windows 11

>>> COREG_LOCAL(
...     "path_to/RefImage.tif",
...     "path_to/SourceImage.tif",
...     grid_res=50,
...     window_size=(256, 256),
...     r_b4match=1,
...     path_out=None,
...     CPUs=1,
...     q=False
... ).correct_shifts()
Calculating footprint polygon and actual data corner coordinates for reference image...
Bounding box of calculated footprint for reference image:
        (474683.0, 6686905.0, 476442.0, 6691721.0)
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:
        (474682.96922200004, 6686904.749473, 476441.569222, 6691721.099473)
Matching window position (X,Y): 475562.28606702085/6689313.000074867
Initializing tie points grid...
Equalizing pixel grids and projections of reference and target image...
Warping progress     |==================================================| 100.0% Complete  => 0:00:18
Warping progress     |==================================================| 100.0% Complete  => 0:00:01
Calculating tie point grid (16692 points) using 1 CPU cores...
Traceback (most recent call last):
  File "<stdin>", line 10, in <module>
  File "d:\pygis\Lib\site-packages\arosics\CoReg_local.py", line 821, in correct_shifts
    self.calculate_spatial_shifts()
  File "d:\pygis\Lib\site-packages\arosics\CoReg_local.py", line 491, in calculate_spatial_shifts
    self._tiepoint_grid.get_CoRegPoints_table()
  File "d:\pygis\Lib\site-packages\arosics\Tie_Point_Grid.py", line 345, in get_CoRegPoints_table
    Parallel(n_jobs=self.CPUs, **kw_parallel)(
    ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
TypeError: Parallel.__init__() got an unexpected keyword argument 'return_as'
danschef commented 6 months ago

Thanks for reporting, looks like an incompatibility with the latest version of joblib. Could you check which version of joblib is installed (conda list joblib).

leikas227 commented 6 months ago

That was version: 1.3.2

leikas227 commented 6 months ago

I tried it also on my home computer, and it seems to work there, at least for a small area. The joblib version I read out earlier was from that computer. The laptop, where it doesn’t work, has the joblib version 1.2.0