# B-Spline interpolation and its mask handling GridR implements Cardinal B-Spline interpolation following {cite}`briand2018theory`. This page focuses on what makes B-Spline masking different from the other interpolators: the recursive prefiltering step propagates invalid samples in a predictable, exponentially decaying pattern that GridR handles automatically via mask dilation. **What you'll learn** - Why B-Spline prefiltering needs special mask handling - The role of the ``epsilon`` parameter - The role of the ``mask_influence_threshold`` parameter - How to set the threshold for outlier values - Reading the diamond-shaped invalidity region in your output ## Masking with a Cardinal B-Spline Interpolator GridR implements Cardinal B-spline interpolation based on {cite}`briand2018theory`. While we won't detail their algorithms here, we focus on GridR's unique mask handling features not covered in the original paper. The interpolation process requires both causal and anti-causal recursive prefiltering applied to the input image. While these filters theoretically require infinite sums, GridR approximates them using finite sums as proposed in the reference article. This approximation relies on either the image's immediate neighborhood or extrapolated data through boundary conditions (detailed further below). The recursive filters used in B-spline interpolation exhibit an **exponential decay** property. When a pixel is invalid, its influence propagates to neighboring pixels through the recursive filtering process, but this influence decreases exponentially with distance. For a pole $z_k$ (where $-1 < z_k < 0$), the exponential filter has the form: $$h^{(z_k)}_j = \frac{z_k}{z_k^2 - 1} z_k^{|j|} (Eq. 21)$$ This shows that the 1D influence at distance $j$ pixels decays as $|z_k|^{|j|}$. In practice, the dominant pole (largest absolute value, typically $z_1$) determines the decay rate. To determine at what distance $d$ the 1D influence becomes negligible, we solve: $$|z_k|^{|d|} \le s$$ where $d$ is the distance in pixels and $s$ is the acceptable residual influence threshold. This gives: $$d \ge \frac{ln(s)}{ln(|z_i|)}$$ Since the prefiltering is performed separably on rows and column, the 2D influence can be expressed as : $$Influence(i,j) \propto |z_k|^{|i|} \times |z_k|^{|j|} = |z_k|^{|i| + |j]}$$ In the 2D case we then have to consider the Manhattan distance (not Euclidian). In practice, this yields a diamond-shaped (45° oriented square) influence zone centered around the invalid pixel's position. Instead of filtering the validity mask (which can produce overshoots artifacts), a more robust approach is to **dilate the invalid area within the validity mask** by the computed influence radius using a Manhattan distance structuring element. Additionally, mask elements corresponding to image pixels used to approximate the infinite sums - refered as domain extension in the article - are automatically invalidated during the mask prefiltering process. This invalidation occurs after the previously mentioned propagation of invalid mask elements. **It is important to note that the residual influence threshold $s$ is relative to the invalid pixel value itself. For extreme outlier values (e.g. contaminated pixel with 10000 in a 0-255 image), even a small relative threshold (e.g., $s = 10^{-3} = 0.1%$) can correspond to a significant absolute contamination ($10000 \times 10^{-3} = 10$), requiring a much larger dilatation radius than for typical data values. When dealing with aberrant values, the threshold $s$ should be chosen accodingly to ensure the absolute contamination remains acceptable for your application.** Let's examine how GridR handles masking for B-spline prefiltering. Additionnaly we will notice we need to pass two additional arguments for B-spline interpolation: `epsilon` and `mask_influence_threshold`. The `epsilon` parameter defines the acceptable error when approximating the finite sums, directly affecting the required margins for prefiltering on each side (referred to as *truncation index* in the reference article). For a fixed B-spline order, the smaller this value, the larger the required margin becomes. Note that the total required margin is the sum of both the prefiltering margin and the margin needed for pure interpolation, corresponding to the interpolation kernel radius. The `mask_influence_threshold` parameter is particularly important here and corresponds to the the residual influence threshold $s$ previously mentioned. We now demonstrate how an aberrant point propagates through both the image and validity mask, comparing results with and without GridR's specialized B-spline prefiltering mask handling. ```python # Create a grid centered on the masked position # Here we use a larger input domain as previously in order to demonstrate the # behaviour of prefiltering. zoom = 4 y = masked_pos[0] + np.linspace(-20, 20, 40*zoom+1) x = masked_pos[1] + np.linspace(-20, 20, 40*zoom+1) xx, yy = np.meshgrid(x, y) # We will use the same input as before but we need to copy it as the B-spline # prefiltering will operate in place. array_in_truth = np.copy(array_in) # First compute resampling with all valid data without altering # the element going to be invalid array_out_womask_truth, array_mask_out_womask_truth = array_grid_resampling( interp="bspline3", interp_kwargs={"epsilon": 1e-6, "mask_influence_threshold": 1}, array_in=array_in_truth, grid_row=yy, grid_col=xx, grid_resolution=(1,1), array_out=None, array_out_mask=True, nodata_out=0, ) # Perform a second computation after making the target element aberrant, # without activating array masking, to observe how the unmasked aberrant # value affects the result. anomalous_value = 1e9 array_in_womask = np.copy(array_in) array_in_womask[*masked_pos] = anomalous_value array_out_womask, array_mask_out_womask = array_grid_resampling( interp="bspline3", interp_kwargs={"epsilon": 1e-6, "mask_influence_threshold": 1}, array_in=array_in_womask, grid_row=yy, grid_col=xx, grid_resolution=(1,1), array_out=None, array_out_mask=True, nodata_out=0, ) # For plot purpose we will set to nodata values that differs more than 1 # in absolute value from the reference output ref_masked = np.abs(array_out_womask - array_out_womask_truth) > 1 array_out_womask_custom = np.copy(array_out_womask) array_out_womask_custom[ref_masked] = 0 ``` /home/docs/checkouts/readthedocs.org/user_builds/gridr/checkouts/latest/notebooks/../python/gridr/core/grid/grid_resampling.py:1784: UserWarning: Standalone mode is enabled : trust_padding will be ignored with `boundary_condition` set as `None` ) = standalone_preprocessing( /tmp/ipykernel_4428/495858054.py:3: RuntimeWarning: invalid value encountered in log 'wo_in_mask_log': np.log(array_out_womask), ![005_bspline_masking_01.png](grid_resampling_005_bspline_masking_files/005_bspline_masking_01.png) Once more we have several plots to comment. From left to right, - The first plot corresponds to a reference resampling with no abberation at all. - The second plot shows the resampling result after modifying the input data to create an aberrant center element. While we can observe some propagation of this aberrant point, the extreme value ($10^9$ - exaggerated for demonstration purposes) makes it impossible to display the result's full dynamic range using linear scaling in an 8-bit grayscale format. Thus motivating the next plot. - The third plot shows the same data than the previous one but using a logarithm scaling for display. It clearly shows the propagation of the aberration while making it possible to guess the region were the image is valid. You can also notice oscillation within the propagation area: they are caracterial of prefiltering around a disruptive point - negative values are here showed as NaN explaining why the same gray color is used for all of them. - The fourth plot comes as complement of the third one. Here it shows the same data but using a logarithm scaling on its absolute value for display. This serves as observing the exponential decay of the influence of the anomalous value. - The last plot shows the same result data but this time the pixels whose value differs from the reference of 0.1 are masked out for display. While the value of 0.1 is arbitrary, we can experimentally check the obtained distance for the corresponding threshold $s = 0.1 * 1e^{-9} = 1e^{-10}$. ```python # Experienced distance - compute it from vertical diameter and considering the zoom_factor (4x) zoom = 1. / (y[1]-y[0]) left = np.min(np.where(ref_masked == 1)[0]) right = np.max(np.where(ref_masked == 1)[0]) d = (right - left) / (2 * zoom) print(f"Measured influence distance : {d}") # Theoretical distance bspline3_z1 = -2.679491924311227e-1 d_th = np.log(1./anomalous_value) / np.log(np.abs(bspline3_z1)) print(f"Theretical influence distance : {d_th}") ``` Measured influence distance : 15.5 Theretical influence distance : 15.735708700586859 An interesting point observed in the last plot is that, since the grid performs an exact 4x zoom, every fourth point targets a whole source coordinate. For these points, the interpolated values are not contaminated by aberrant values. Let's now examine how GridR handles masks by activating the usage of an input mask. We set the `mask_influence_threshold` parameter such that the absolute residual does not exceed a value of 1, aligning with the previously obtained experimental masking. Since this threshold is relative to the anomalous value, we set it to $1. / 10^{-9} = 10^{-10}$, resulting in a computed distance of 15.73 as shown earlier. ```python # Perform a computation after making the target element aberrant, # with activating array masking array_in_wmask_mask = np.ones(array_in.shape, dtype=np.uint8) array_in_wmask_mask[*masked_pos] = 0 # Copy mask as it will be updated in place array_in_wmask = np.copy(array_in) array_in_wmask[*masked_pos] = anomalous_value # Here we set mask_influence_threshold to match the targeted accepted # residual of the anomalous value. # We want it to correspond to 0.1 in absolute value array_out_wmask, array_mask_out_wmask = array_grid_resampling( interp="bspline3", interp_kwargs={"epsilon": 1e-6, "mask_influence_threshold": 1./anomalous_value}, array_in=array_in_wmask, array_in_mask=array_in_wmask_mask, grid_row=yy, grid_col=xx, grid_resolution=(1,1), array_out=None, array_out_mask=True, nodata_out=0, ) array_out_wmask_masked = np.copy(array_out_wmask) array_out_wmask_masked[array_mask_out_wmask == 0] = 0 ``` ![005_bspline_masking_02.png](grid_resampling_005_bspline_masking_files/005_bspline_masking_02.png) In this case, the diamond-shaped mask is internally computed during the prefiltering stage and subsequently used as an updated mask for the interpolation process. Consequently, any invalid data within the mask are automatically set to the nodata value specified through the `nodata_out` parameter.