Fixed ACF parameters and added interpolation

surfalize/autocorrelation.py

@@ -4,7 +4,7 @@
 from mpl_toolkits.axes_grid1 import make_axes_locatable
 
 from .cache import CachedInstance, cache
-from .mathutils import interpolate_line_on_2d_array
+from .mathutils import interpolate_line_on_2d_array, argmin_all, argmax_all, argclosest
 
 
 class AutocorrelationFunction(CachedInstance):
@@ -24,7 +24,7 @@ def __init__(self, surface):
         self._surface = surface
         self._current_threshold = None
         self.data = self.calculate_autocorrelation()
-        self._center = np.array(self.data.shape) // 2
+        self.center = np.array(self.data.shape) // 2
 
     def calculate_autocorrelation(self):
         data = self._surface.center().data
@@ -33,9 +33,10 @@ def calculate_autocorrelation(self):
         acf_data = np.fft.fftshift(np.fft.ifft2(data_fft * np.conj(data_fft)).real / data.size)
         return acf_data
 
-    def old_calculate_distances(self, s=0.2):
+    @cache
+    def _calculate_decay_lengths(self, s):
         """
-        Calculates the distances of the 2d autocorrelation function of the surface height data and
+        Calculates the decay lengths of the 2d autocorrelation function of the surface height data and
         extracts the indices of the points of minimum and maximum decay.
 
         Parameters
@@ -50,63 +51,37 @@ def old_calculate_distances(self, s=0.2):
         -------
         None
         """
-        self.clear_cache()
-        self._current_threshold = s
-
-        threshold = s * self.acf_data.max()
-        mask = (self.acf_data < threshold)
-
-        # Find the center point of the array
-        self.center = np.array(self.acf_data.shape) // 2
-
-        # Invert mask because the function considers all 0 values to be background
-        labels, _ = ndimage.label(~mask)
-        feature_center_id = labels[self.center[0], self.center[1]]
-        mask = (labels != feature_center_id)
-
-        # Find the indices of the True values in the mask
-        indices = np.argwhere(mask)
-        # Calculate the Euclidean distance from each index to the center
-        distances = np.linalg.norm(indices - self.center, axis=1)
-        # Find the index with the smallest distance
-        self._idx_min = indices[np.argmin(distances)]
-
-        # Find the indices of the True values in the mask
-        indices = np.argwhere(~mask)
-        # Calculate the Euclidean distance from each index to the center
-        distances = np.linalg.norm(indices - self.center, axis=1)
-        self._idx_max = indices[np.argmax(distances)]
-
-    def _calculate_distances(self, s):
-        self.clear_cache()
-        self._current_threshold = s
-        self._absolute_threshold = s * self.data.max()
-        mask = (self.data < self._absolute_threshold)
-
-        # Invert mask because the function considers all 0 values to be background
-        # Use image segementation to identify contiguous regions in the data
-        labels = label(~mask)
-        feature_center_id = labels[self.center[0], self.center[1]]
-        # Remove all but the central region
-        mask = (labels == feature_center_id)
-        # Currently, the outer pixels of the center region are just above the threshold
-        # So we expand the mask by one pixel to obtain the pixels that are just below the threshold
-        expanded_mask = expand_labels(mask, distance=1)
-        # Using xor we obtain the edge of the region with values just below the threshold
-        edge = expanded_mask ^ mask
-
-        # Calculate the Euclidean distance from each edge pixel to the center
-        indices = np.argwhere(edge)
-
-        # We need to calculate the distances in units for this to work for surfaces with unequal spacing
-        # in x and y (looking at you OPD files, why???)
-        distances_xy_px = indices - self.center
-        distances_xy_units = distances_xy_px * np.array([self._surface.step_y, self._surface.step_x])
+        threshold = s * self.data.max()
+
+        mask = self.data > threshold
+        labels, _ = ndimage.label(mask)
+        region = labels == labels[*self.center]
+        edge = region ^ ndimage.binary_dilation(region, iterations=1)
+
+        idx_edge = np.argwhere(edge)
+        distances_xy_px = idx_edge - self.center
+        step_array = np.array([self._surface.step_y, self._surface.step_x])
+        distances_xy_units = distances_xy_px * step_array
         distances = np.linalg.norm(distances_xy_units, axis=1)
+        all_idx_min = idx_edge[argmin_all(distances)]
+        all_idx_max = idx_edge[argmax_all(distances)]
+
+        idx_min = all_idx_min[np.argmin(self.data[all_idx_min[:, 0], all_idx_min[:, 1]])]
+        idx_max = all_idx_max[np.argmin(self.data[all_idx_max[:, 0], all_idx_max[:, 1]])]
+
+        length_min = np.hypot(*((idx_min - self.center) * step_array))
+        length_max = np.hypot(*((idx_max - self.center) * step_array))
+
+        n_points = 1000
+        interpolated_line_x = np.linspace(0, length_min, n_points)
+        interpolated_line_y = interpolate_line_on_2d_array(self.data, self.center, idx_min, num_points=n_points)
+        shortest_decay_length = interpolated_line_x[argclosest(threshold, interpolated_line_y)]
+
+        interpolated_line_x = np.linspace(0, length_max, n_points)
+        interpolated_line_y = interpolate_line_on_2d_array(self.data, self.center, idx_max, num_points=n_points)
+        longest_decay_length = interpolated_line_x[argclosest(threshold, interpolated_line_y)]
 
-        # Find the index with the smallest and largest distance
-        self.idx_min = indices[np.argmin(distances)]
-        self.idx_max = indices[np.argmax(distances)]
+        return shortest_decay_length, longest_decay_length
 
 
     @cache
@@ -129,10 +104,7 @@ def Sal(self, s=0.2):
         Sal : float
             autocorrelation length.
         """
-        if self._current_threshold != s:
-            self._calculate_distances(s)
-        dy, dx = np.abs(self._idx_min[0] - self.center[0]), np.abs(self._idx_min[1] - self.center[1])
-        Sal = np.hypot(dx * self._surface.step_x, dy * self._surface.step_y) - self._surface.step_x/2
+        Sal, _ = self._calculate_decay_lengths(s)
         return Sal
 
     @cache
@@ -156,10 +128,8 @@ def Str(self, s=0.2):
         Str : float
             texture aspect ratio.
         """
-        if self._current_threshold != s:
-            self._calculate_distances(s)
-        dy, dx = np.abs(self._idx_max[0] - self.center[0]), np.abs(self._idx_max[1] - self.center[1])
-        Str = self.Sal() / (np.hypot(dx * self._surface.step_x, dy * self._surface.step_y) - self._surface.step_x/2)
+        shortest_decay_length, longest_decay_length = self._calculate_decay_lengths(s)
+        Str = shortest_decay_length / longest_decay_length
         return Str
 
     def plot_autocorrelation(self, ax=None, cmap='jet', show_cbar=True):

tests/test_parameters.py

@@ -41,11 +41,11 @@ def test_Sdq(surface):
 
 
 def test_Sal(surface):
-    assert surface.Sal() == pytest.approx(4.25, abs=EPSILON)
+    assert surface.Sal() == pytest.approx(4.274174, abs=EPSILON)
 
 
 def test_Str(surface):
-    assert surface.Str() == pytest.approx(0.120735, abs=EPSILON)
+    assert surface.Str() == pytest.approx(0.121940, abs=EPSILON)
 
 
 def test_Sk(surface):