added filtering in the cost function. It can be (de)activted with the...

added filtering in the cost function. It can be (de)activted with the filtered_cost argument which defaults to True

nabu/preproc/alignment.py

@@ -611,7 +611,12 @@ class CenterOfRotationAdaptiveSearch(CenterOfRotation):
     the apodisation function, which is a gaussian, is moved to a new guess position.
     The lenght of the step, by which the gaussian is moved, and its sigma are
     obtained by multiplying the shortest distance from the left or right border with
-    a self.step_fraction and  self.sigma_fraction factors which ensure global overlapping
+    a self.step_fraction and  self.sigma_fraction factors which ensure global overlapping.
+    for each step a region around the CoR  of each image is selected, and the regions of the two images
+    are compared to  calculate a cost function. The value of the cost function, at its minimum
+    is used to select the best step at which the CoR is taken as final result.
+    The option filtered_cost= True (default) triggers the filtering (accordin to low_pass and high_pass)
+    of the two images which are used for he cost function.
     """
 
     sigma_fraction = 1.0 / 4.0
@@ -628,6 +633,7 @@ class CenterOfRotationAdaptiveSearch(CenterOfRotation):
         high_pass=None,
         low_pass=None,
         margins=None,
+        filtered_cost=True,
     ):
         """Find the Center of Rotation (CoR), given two images.
 
@@ -670,13 +676,15 @@ class CenterOfRotationAdaptiveSearch(CenterOfRotation):
         peak_fit_radius: int, optional
             Radius size around the max correlation pixel, for sub-pixel fitting.
             Minimum and default value is 1.
-        low_pass: float or sequence of two floats
+        low_pass: float or sequence of two floats.
             Low-pass filter properties, as described in `nabu.misc.fourier_filters`
         high_pass: float or sequence of two floats
-            High-pass filter properties, as described in `nabu.misc.fourier_filters`
+            High-pass filter properties, as described in `nabu.misc.fourier_filters`.
         margins:  None or a couple of floats or ints
             if margins is None or in the form of  (margin1,margin2) the search is done between margin1 and  dim_x-1-margin2.
-            If left to None then by default (margin1,margin2)  = ( 10, 10 )
+            If left to None then by default (margin1,margin2)  = ( 10, 10 ).
+        filtered_cost: boolean.
+            True by default. It triggers the use of filtered images in the calculation of the cost function.
 
         Raises
         ------
@@ -716,9 +724,27 @@ class CenterOfRotationAdaptiveSearch(CenterOfRotation):
 
         roi_yxhw = self._determine_roi(img_1.shape, roi_yxhw)
 
+        if filtered_cost and (low_pass is not None or high_pass is not None):
+            img_filter = fourier_filters.get_bandpass_filter(
+                img_1.shape[-2:],
+                cutoff_lowpass=low_pass,
+                cutoff_highpass=high_pass,
+                use_rfft=__have_scipy__,
+                data_type=self.data_type,
+            )
+            # fft2 and iff2 use axes=(-2, -1) by default
+            img_filtered_1 = local_ifftn(local_fftn(img_1, axes=(-2, -1)) * img_filter, axes=(-2, -1)).real
+            img_filtered_2 = local_ifftn(local_fftn(img_2, axes=(-2, -1)) * img_filter, axes=(-2, -1)).real
+        else:
+            img_filtered_1 = img_1
+            img_filtered_2 = img_2
+
         img_1 = self._prepare_image(img_1, roi_yxhw=roi_yxhw, median_filt_shape=median_filt_shape)
         img_2 = self._prepare_image(img_2, roi_yxhw=roi_yxhw, median_filt_shape=median_filt_shape)
 
+        img_filtered_1 = self._prepare_image(img_filtered_1, roi_yxhw=roi_yxhw, median_filt_shape=median_filt_shape)
+        img_filtered_2 = self._prepare_image(img_filtered_2, roi_yxhw=roi_yxhw, median_filt_shape=median_filt_shape)
+
         dim_radio = img_1.shape[1]
 
         if margins is None:
@@ -734,7 +760,7 @@ class CenterOfRotationAdaptiveSearch(CenterOfRotation):
 
         if lim_2 <= lim_1:
             message = (
-                "Image shape or cropped selection too small for global search. After removal of the margins the search limit collide. The cropped size is %d\n"
+                "Image shape or cropped selection too small for global search. After removal of the margins the search limits collide. The cropped size is %d\n"
                 % (dim_radio)
             )
             raise ValueError(message)
@@ -789,8 +815,8 @@ class CenterOfRotationAdaptiveSearch(CenterOfRotation):
             M1 = int(round(cor_position + img_1.shape[1] // 2)) - int(round(tmp_sigma))
             M2 = int(round(cor_position + img_1.shape[1] // 2)) + int(round(tmp_sigma))
 
-            piece_1 = img_1[:, M1:M2]
-            piece_2 = img_2[:, img_1.shape[1] - M2 : img_1.shape[1] - M1]
+            piece_1 = img_filtered_1[:, M1:M2]
+            piece_2 = img_filtered_2[:, img_1.shape[1] - M2 : img_1.shape[1] - M1]
             piece_1 = piece_1 - piece_1.mean()
             piece_2 = piece_2 - piece_2.mean()
             energy = np.array(piece_1 * piece_1 + piece_2 * piece_2, "d").sum()

nabu/preproc/tests/test_alignment.py

@@ -369,7 +369,7 @@ class TestCor(object):
 
         CoR_calc = alignment.CenterOfRotationAdaptiveSearch()
 
-        cor_position = CoR_calc.find_shift(radio1, radio2, low_pass=1, high_pass=20)
+        cor_position = CoR_calc.find_shift(radio1, radio2, low_pass=1, high_pass=20, filtered_cost=True)
 
         print("Found cor_position", cor_position)
 
@@ -391,7 +391,7 @@ class TestCor(object):
 
         CoR_calc = alignment.CenterOfRotationAdaptiveSearch()
 
-        cor_position = CoR_calc.find_shift(radio1, radio2, low_pass=1, high_pass=20, margins=(100, 10))
+        cor_position = CoR_calc.find_shift(radio1, radio2, low_pass=1, high_pass=20, margins=(100, 10), filtered_cost=False)
 
         print("Found cor_position", cor_position)