# Using the GridR's core grid resampling API This guide shows how to use the core `array_grid_resampling` GridR's feature. It covers the following topics: - Performing basic grid-based resampling - Working with spatial windows - Controlling output grid resolution - Handling nodata values and generating output validity masks - Using masks to control the validity of input grids and input data, with a dedicated focus on mask handling for the cardinal B-Spline - Boundary condition (related to the standalone mode) - Advanced Usages: Pipeline integration ## Setting things up First we have to import the `array_grid_resampling` method. We also import the well known mandrill raster image. ```python from gridr.core.grid.grid_resampling import array_grid_resampling from gridr.misc.mandrill import mandrill ``` Please note the mandrill image has 3 color's channels. Here we limit the image to the first channel and convert its type to 64 bits floating points. ```python image = mandrill[0] grid_dtype = np.float64 data_in_dtype = np.float64 data_out_dtype = np.float64 array_in = image.astype(data_in_dtype) ``` ## Identity grid transform We start by defining the idendity transform grid. That's not really impressive but we will go further to illustrate some variants (windowing, zooming, ...) Let's first create the idendity transformation grid. ```python # create identity grid if image.ndim == 2: x = np.arange(0, image.shape[0], dtype=grid_dtype) y = np.arange(0, image.shape[1], dtype=grid_dtype) xx, yy = np.meshgrid(x, y) ``` ### Simple idendity transform The `array_grid_resampling` method provides the `array_out` argument in order to use a preallocated C-contiguous array as output. It is quite convenient when implementing a tiling chain that uses that method. For now, we choose to keep it simple : by setting `array_out` to `None` we let the method allocate the output buffer for us. We have defined the rows and columns grids to be at full resolution, ie. there is one coordinate in the grids for each pixel in the input image. To stay at that resolution we set `grid_resolution` to 1 in both directions. The `array_grid_resampling` function returns a tuple of two NumPy arrays. For now, we will focus on the first one — it contains the allocated output buffer, which holds the result of the interpolation. For the interpolation method, we will primarily use the `cubic` identifier, which refers to the Optimized Bicubic interpolation algorithm. Other available interpolators include `nearest`, `linear`, `bspline3`, `bspline5`, `bspline7`, `bspline9`, and `bspline11`. Note that the B-spline interpolators require additional parameters, which are detailed further in this guide within the masking section. ```python # Lets call the grid resampling array_out_id, _ = array_grid_resampling( interp="cubic", array_in=image.astype(data_in_dtype), grid_row=yy, grid_col=xx, grid_resolution=(1,1), array_out=None, win=None, array_in_mask=None, grid_mask=None, array_out_mask=None, nodata_out=None, array_in_origin=None, ) ``` ![001.png](grid_resampling_001_files/001.png) In the following test we compare input with output. It must succeed ! ```python np.testing.assert_almost_equal(array_in, array_out_id, decimal=10) ``` ### Apply a window The `array_grid_resampling` method provides the `win` parameter in order to set a window to limit the region to compute. The indices used in the `win` definition refers to indices considering the full resolution output. Below, we first define a window centered on the mandrill's left eye. Since we are considering an identity grid, the grid coordinates align with the input array's image coordinates ```python win_center = np.array((58, 175)) # left eye win_shape = np.array((100, 100)) win = np.array((win_center - win_shape//2, win_center - win_shape//2 + win_shape - 1)).T ``` ![002.png](grid_resampling_001_files/002.png) ```python array_out_id_win_1_1, _ = array_grid_resampling( interp="cubic", array_in=array_in, grid_row=yy, grid_col=xx, grid_resolution=(1,1), array_out=None, win=win, array_in_mask=None, grid_mask=None, array_out_mask=None, nodata_out=None, ) ``` ![003.png](grid_resampling_001_files/003.png) ### Playing with resolution to apply a zoom The `array_grid_resampling` method uses the `grid_resolution` parameter to insert additional nodes between existing grid nodes. This insertion is achieved through bilinear interpolation, with the oversampling factors for both rows and columns specified in the grid_resolution parameter. While the primary purpose of the `grid_resolution` is to control computation time and grid storage, in this case, we use it to perform a zoom operation, increasing the row resolution by a factor of 2 and the column resolution by a factor of 5. At the same time, we want to maintain the same input region as before. To achieve this, the window definition must be adjusted to accommodate the new full-resolution output. ```python # Set resolution to the target zoom factors resolution = (2, 3) # The full resolution grid has (Nrow-1)*resolution_row + 1 rows (apply same logic on column) full_output_shape = (np.asarray(array_in.shape) - 1)*resolution + 1 # setting output window to only consider a region of interest # here we construct a window centered on the left eye in the full resolution output grid win_center = np.array((58, 175)) * np.array(resolution) win_shape = np.array((100, 100)) * np.array(resolution) win = np.array((win_center - win_shape//2, win_center - win_shape//2 + win_shape - 1)).T # Lets call the grid resampling array_out_id_zoom_2_5, _ = array_grid_resampling( interp="cubic", array_in=array_in, grid_row=yy, grid_col=xx, grid_resolution=resolution, array_out=None, win=win, array_in_mask=None, grid_mask=None, array_out_mask=None, nodata_out=None, ) ``` ![004.png](grid_resampling_001_files/004.png) In order to illustrate the change of interpolation method, we will use here the `nearest` identifier in order to adress a nearest neighbor interpolation. ```python # Lets call the grid resampling array_out_id_zoom_2_5_nearest, _ = array_grid_resampling( interp="nearest", array_in=array_in, grid_row=yy, grid_col=xx, grid_resolution=resolution, array_out=None, win=win, array_in_mask=None, grid_mask=None, array_out_mask=None, nodata_out=None, ) # Let's show it plot_im(array_out_id_zoom_2_5_nearest, None, prefix='004_nearest') ``` ![004_nearest.png](grid_resampling_001_files/004_nearest.png) ## Geometry transformation examples In that part we provide some examples of simple non identity transformation. ### Translation transform Let's begin with a simple translation on the column coordinates by shifting the image to the left. ```python # First create identity grid if image.ndim == 2: x = np.arange(0, image.shape[0], dtype=grid_dtype) y = np.arange(0, image.shape[1], dtype=grid_dtype) xx, yy = np.meshgrid(x, y) xx += 50.5 # Lets call the grid resampling array_out_translation_left, _ = array_grid_resampling( interp="cubic", array_in=array_in, grid_row=yy, grid_col=xx, grid_resolution=(1,1), array_out=None, win=None, array_in_mask=None, grid_mask=None, array_out_mask=None, nodata_out=None, ) ``` ![005.png](grid_resampling_001_files/005.png) ```python array_out_translation_left[:, -53:] ``` array([[124.875 , nan, nan, ..., nan, nan, nan], [130. , nan, nan, ..., nan, nan, nan], [102.8125, nan, nan, ..., nan, nan, nan], ..., [102.5625, nan, nan, ..., nan, nan, nan], [ 86.5625, nan, nan, ..., nan, nan, nan], [ 3.5625, nan, nan, ..., nan, nan, nan]], shape=(512, 53)) As you can see, the right strip border is filled with `NaN` values. It is due to the fact that the grid target coordinates went outside of the image's definition domain. Please note that the filling with `NaN` is not the default behaviour and it occured because we set the `nodata_out` parameter to `None`. Let's replay it by setting `nodata_out` to its default value, ie 0. As you will see, the right strip border is now filled with zeros. We also take this opportunity to introduce the use of the `array_mask_out` parameter and examine the second array returned by the function. We will verify that the mask correctly matches the `nodata` region of the resampled image. Here we set the `array_out_mask` parameter to `True` ; that's a convenient way to let the `array_grid_resampling` method allocate the right sized `unsigned byte` (ie. uint8) array and return it. Please note that we can also directly pass a preallocated array similarly to the `array_out` parameter. ```python array_out_translation_left, array_out_translation_left_mask = array_grid_resampling( interp="cubic", array_in=array_in, grid_row=yy, grid_col=xx, grid_resolution=(1,1), array_out=None, win=None, array_in_mask=None, grid_mask=None, array_out_mask=True, nodata_out=0., ) array_out_translation_left[:, -53:] ``` array([[124.875 , 0. , 0. , ..., 0. , 0. , 0. ], [130. , 0. , 0. , ..., 0. , 0. , 0. ], [102.8125, 0. , 0. , ..., 0. , 0. , 0. ], ..., [102.5625, 0. , 0. , ..., 0. , 0. , 0. ], [ 86.5625, 0. , 0. , ..., 0. , 0. , 0. ], [ 3.5625, 0. , 0. , ..., 0. , 0. , 0. ]], shape=(512, 53)) Let's display the output validity mask ```python array_out_translation_left_mask[:, -53:] ``` array([[1, 0, 0, ..., 0, 0, 0], [1, 0, 0, ..., 0, 0, 0], [1, 0, 0, ..., 0, 0, 0], ..., [1, 0, 0, ..., 0, 0, 0], [1, 0, 0, ..., 0, 0, 0], [1, 0, 0, ..., 0, 0, 0]], shape=(512, 53), dtype=uint8) ![005-mask-out.png](grid_resampling_001_files/005-mask-out.png) ### Rotation transform Here we define a grid in order to apply a 45 degree rotation. In order to have the full image rotated please note that we have to extend the grid size. ```python center = np.asarray(array_in.shape) // 2 theta = np.pi/4. # we know the image is square - here we extend the output grid in order to cover # the full image ext = image.shape[0] * (2**0.5 - 1) // 2 # First create identity grid if image.ndim == 2: x = np.arange(0, image.shape[0] + 2 * ext, dtype=grid_dtype) - ext y = np.arange(0, image.shape[1] + 2 * ext, dtype=grid_dtype) - ext xx, yy = np.meshgrid(x , y) # Apply the rotation xrot = (xx - center[1]) * np.cos(theta) - (yy - center[0]) * np.sin(theta) + center[1] yrot = (xx - center[1]) * np.sin(theta) + (yy - center[0]) * np.cos(theta) + center[0] ``` As a demonstration purpose we here preallocate the array_out buffer. ```python # We have to init the out buffer with the right shape array_out = np.zeros(xx.shape, dtype=data_out_dtype) # Lets call the grid resampling _, mask_out = array_grid_resampling( interp="cubic", array_in=array_in, grid_row=yrot, grid_col=xrot, grid_resolution=(1,1), array_out=array_out, win=None, array_in_mask=None, grid_mask=None, array_out_mask=True, nodata_out=None, ) ``` ![006.png](grid_resampling_001_files/006.png) ### Specifying an output target window In the previous step we opt to pre-allocate the output buffer for the image. What if the requested production window is smaller that the `array_out` buffer ? ```python # Set a window to focus on the left eye (we get the coordinates from the previous resampled image) win = np.asarray([[250, 320], [120, 220]]) # We have to init the out buffer with an sufficient shape array_out = np.zeros((200, 200), dtype=data_out_dtype) # Lets call the grid resampling _, mask_out = array_grid_resampling( interp="cubic", array_in=array_in, grid_row=yrot, grid_col=xrot, grid_resolution=(1,1), array_out=array_out, win=win, array_in_mask=None, grid_mask=None, array_out_mask=True, nodata_out=None, ) ``` ![007.png](grid_resampling_001_files/007.png) It’s worth noting that even without passing a preallocated output validity mask, the internally generated one will still match the shape of the `array_out` buffer. As shown above, the resampled window is written starting at the origin of `array_out` (i.e., at [0, 0]). In some cases - for example, when working with tiles that share a common output buffer - it may be necessary to write the result to a specific region within `array_out`. This is precisely the purpose of the `array_out_win` parameter. ```python # Set a window to focus on the left eye (we get the coordinates from the previous resampled image) win = np.asarray([[250, 320], [120, 220]]) # We have to init the out buffer with an sufficient shape # Here we will also demonstrate that the resampling only modify the target window # For that purpose we do not initialize with zeros but with arange array_out = np.arange(200*200, dtype=data_out_dtype).reshape((200, 200)) % 255 # We also preallocate the validity mask for the same purpose mask_out = np.arange(array_out.size, dtype=np.uint8).reshape(array_out.shape) % 10 # Let's define where we want to place the resampled window array_out_win_origin = np.asarray([25, 40]) array_out_win = np.asarray([ [array_out_win_origin[0], array_out_win_origin[0] + win[0][1] - win[0][0]], [array_out_win_origin[1], array_out_win_origin[1] + win[1][1] - win[1][0]]]) # Lets call the grid resampling _, _ = array_grid_resampling( interp="cubic", array_in=array_in, grid_row=yrot, grid_col=xrot, grid_resolution=(1,1), array_out=array_out, array_out_win=array_out_win, # <= new parameter win=win, array_in_mask=None, grid_mask=None, array_out_mask=mask_out, nodata_out=None, ) ``` ![008.png](grid_resampling_001_files/008.png) Note that this example also demonstrates that operating on a specific window only affects the corresponding region in the output array, leaving all other areas unchanged. ## Working with masks In this part we will see how to use masking. The `array_grid_resampling` method provides several kind of masking for both the grid and the input array. ### General convention GridR resampling functionalities uses the following convention for mask : - the mask are considered unsigned 8 bits integers taking value in {0, 1}. - a value of 0 is considered as non valid - a value of 1 is considered as valid ### Grid masking #### Grid masking with a raster mask Here we use an 8 bit raster having the same shape of the grids in order to devalidate points in the grid. To keep it simple we will demonstrate using an identity transform. ```python # First create identity grid if image.ndim == 2: x = np.arange(0, image.shape[0], dtype=grid_dtype) y = np.arange(0, image.shape[1], dtype=grid_dtype) xx, yy = np.meshgrid(x, y) ``` We define the mask as following : * the first 10 lines and last 20 lines are masked * the first 30 columns and last 40 columns are masked * a 80x80 window positionned at the grid center are masked Please note the mask applies to the grid and not to the input image. By convention, the mask is considered binary. A value of 0 corresponds to a masked cell, a value of 1 corresponds to a valid cell. ```python # setting a mask to devalidate margins and center grid_mask_valid_value = 1 grid_mask_invalid_value = 0 # the final multiplication here is not necessary but we keep it to show the mechanic grid_mask = np.ones(xx.shape, dtype=np.uint8) * grid_mask_valid_value # masked pixels are positiv ones grid_mask[0:10] = grid_mask_invalid_value grid_mask[-20:] = grid_mask_invalid_value grid_mask[:,0:30] = grid_mask_invalid_value grid_mask[:,-40:] = grid_mask_invalid_value # devalidate center grid_mask[image.shape[0]//2-40:image.shape[0]//2+40, image.shape[1]//2-40:image.shape[1]//2+40] = grid_mask_invalid_value ``` The `nodata_out` is used to set the output value affected by the masking operation. ```python # Lets call the grid resampling array_out, mask_out = array_grid_resampling( interp="cubic", array_in=array_in, grid_row=yy, grid_col=xx, grid_resolution=(1,1), array_out=None, win=None, array_in_mask=None, grid_mask=grid_mask, grid_mask_valid_value=grid_mask_valid_value, array_out_mask=True, nodata_out=0, ) ``` ![010.png](grid_resampling_001_files/010.png) As in the previous example, we pass True to the function to obtain the validity mask output. We then check that it aligns correctly with the grid mask. #### Masking with an inplaced grid nodata value In this case, we avoid providing an additional mask, primarily to limit memory usage. The `array_grid_resampling` method allows specifying a `grid_nodata` parameter, which is used to invalidate grid cells. Cells in the grid whose values match this nodata value (within a fixed absolute tolerance of 1e-5) are treated as masked. Providing a non-`None` value for `grid_nodata` activates the masking mechanism. Note that it is not necessary to assign the nodata value to both grid coordinate arrays (row and column); setting it in either one is sufficient for the corresponding cell to be considered masked. ```python # Set a nodata value grid_nodata = -999. # Apply it to a window of the grid yy[50:200,100:350] = grid_nodata ``` ```python # Lets call the grid resampling # Spoiler alert: this will fail ! try: array_out, _ = array_grid_resampling( interp="cubic", array_in=array_in, grid_row=yy, grid_col=xx, grid_resolution=(1,1), array_out=None, win=None, array_in_mask=None, grid_nodata=grid_nodata, grid_mask=grid_mask, grid_mask_valid_value=grid_mask_valid_value, array_out_mask=None, nodata_out=0, ) except ValueError as e: display(e) ``` ValueError('Only one of `grid_mask` or `grid_nodata` may be provided, not both.') The call failed, and the ValueError exception provides the reason: both masking modes, grid_mask and grid_nodata, cannot be used simultaneously. ```python # Lets remove the grid_mask parameter's definition array_out, _ = array_grid_resampling( interp="cubic", array_in=array_in, grid_row=yy, grid_col=xx, grid_resolution=(1,1), array_out=None, win=None, array_in_mask=None, grid_nodata=grid_nodata, grid_mask=None, grid_mask_valid_value=None, array_out_mask=None, nodata_out=0, ) ``` ![011.png](grid_resampling_001_files/011.png) #### Effect of Grid Resolution on Grid Masking Both the raster mask and the invalid_value grid mask remain consistent with the sampling of the grid. When the resolution differs from 1 in a given direction, interpolation must occur within each grid cell — that is, within the mesh defined by four neighboring grid nodes. In such cases, if at least one non-zero weighted node of the mesh is masked, the interpolated result is marked as invalid. This results in nodata values in the resampled image and corresponding invalid entries (value 0) in the output mask. Let's illustrate this behavior by zooming in on a small region. ```python # create identity grid on the eye if image.ndim == 2: x = np.arange(130, 220, dtype=grid_dtype) y = np.arange(30, 100, dtype=grid_dtype) xx, yy = np.meshgrid(x, y) resolution = np.asarray((2,3)) # lets define the mask using grid_nodata grid_nodata = -99999 yy[30, 35] = grid_nodata array_out, mask_out = array_grid_resampling( interp="cubic", array_in=array_in, grid_row=yy, grid_col=xx, grid_resolution=tuple(resolution), array_out=None, win=None, array_in_mask=None, grid_nodata=grid_nodata, array_out_mask=True, nodata_out=0, ) ``` ![012.png](grid_resampling_001_files/012.png) To invalidate a specific grid position, we've configured the input (30, 35). This grid coordinate directly corresponds to the (60, 105) output position in the output coordinate system. ```python print(mask_out[60-2:60+3,105-3:105+4]) print(f"Number of masked data : {np.where(mask_out==0)[0].shape}") ``` [[1 1 1 1 1 1 1] [1 0 0 0 0 0 1] [1 0 0 0 0 0 1] [1 0 0 0 0 0 1] [1 1 1 1 1 1 1]] Number of masked data : (15,) The masked area, as illustrated, exceeds a single pixel due to a resolution factor of 2 for rows and 3 for columns. The displayed array is centered on the position of the invalid grid node. Notably, this invalidation propagates to its immediate neighborhood, but exclusively for interpolated values where its weight is non-zero. For clearer comprehension, in the schema below, an `x` denotes an invalid grid node, an `o` represents other grid nodes, and both `-` and `*` signify oversampled positions. It's crucial to note that only the `*` positions utilize the `x` (invalid grid node) in their bilinear interpolation calculations. o - - o - - o - * * * * * - o * * x * * o - * * * * * - o - - o - - o ### Input Array masking #### Array masking with a raster mask Here we use an 8 bit raster having the same shape of the input array in order to devalidate points in the input array. Basically if input data is devalidate, interpolation that involves those data must be considered as not valid resulting in output valued to `nodata_out` and the optional `array_out_mask` set to 0 instead of 1. To keep it simple we will demonstrate using a simple zoom transformation transform. ```python array_in_mask = np.ones(array_in.shape, dtype=np.uint8) masked_pos = (60, 160) array_in_mask[*masked_pos] = 0 # In order to check the good interpolation we will also set the array_in to hold an aberrant value save_array_in = array_in[*masked_pos] # Create a grid centered on the masked position y = masked_pos[0] + np.linspace(-3, 3, 100) x = masked_pos[1] + np.linspace(-3, 3, 100) xx, yy = np.meshgrid(x, y) # First compute a reference array_out_womask_truth, _ = array_grid_resampling( interp="cubic", array_in=array_in, grid_row=yy, grid_col=xx, grid_resolution=(1,1), array_out=None, nodata_out=0, ) # Second compute without array masking to watch the effect of the non masked aberrant value array_in[*masked_pos] = 1e7 array_out_womask, _ = array_grid_resampling( interp="cubic", array_in=array_in, grid_row=yy, grid_col=xx, grid_resolution=(1,1), array_out=None, array_out_mask=True, nodata_out=0, ) aberrant_radiation = np.zeros(array_out_womask.shape) aberrant_radiation[np.where(array_out_womask - array_out_womask_truth != 0)] = 1 # Recompute with the mask array_out, mask_out = array_grid_resampling( interp="cubic", array_in=array_in, grid_row=yy, grid_col=xx, grid_resolution=(1,1), array_in_mask=array_in_mask, array_out=None, array_out_mask=True, nodata_out=150, # Here we set nodata_out to 150 in order to have some texture for the display ) # restore array_in array_in[*masked_pos] = save_array_in ``` ![013-array-masked-data.png](grid_resampling_001_files/013-array-masked-data.png) We have several plots to analyze here. The first plot reveals what happens when the anomalous value is not masked. While the 8-bit display dynamic range prevents us from fully appreciating the valid interpolation, this plot effectively showcases the influence of an anomalous value (exaggerated for illustrative purposes) on the interpolation process. The second plot, a binary representation, indicates where the anomalous value "spreads" or "irradiates." Interestingly, some unaffected lines can be seen traversing the "irradiated" zone. The third plot demonstrates the outcome when the anomalous value is successfully masked. We can see a clear correspondence between the no-data pixels and the "irradiated" area. Lastly, the final plot presents the output validity mask, confirming its consistency with our earlier findings. ```python # min and max indices of masked data upper_left = np.min(np.where(mask_out==0), axis=1) lower_right = np.max(np.where(mask_out==0), axis=1) print(f'upper left corner of masked data : {upper_left}') print(f'lower right corner of masked data : {lower_right}') # corresponding grid coordinates print(f'masked data upper left corner grid coordinates : {yy[*upper_left], xx[*upper_left]}') print(f'masked data lower right corner grid coordinates : {yy[*lower_right], xx[*lower_right]}') ``` upper left corner of masked data : [17 17] lower right corner of masked data : [82 82] masked data upper left corner grid coordinates : (np.float64(58.03030303030303), np.float64(158.03030303030303)) masked data lower right corner grid coordinates : (np.float64(61.96969696969697), np.float64(161.96969696969697)) The optimized bicubic interpolator uses a spatial support of ±2 around the target coordinates. This directly explains the masked area's boundaries: row indices [58,62] and column indices [158,162]. Now, let's discuss the unmasked lines that cross the invalid area. To do this, we will plot the grid coordinates, differentiating between the valid and invalid ones. ![013-mask-out-grid-scatter.png](grid_resampling_001_files/013-mask-out-grid-scatter.png) We've used a distinct style in the plot above to highlight coordinates that fall on exact integer values. This specific scenario, where we have whole integer coordinates within our grid, serves to prove that the `GridR` implementation aligns correctly with the "irradiated" area. It also shows that it avoids inadvertently masking pixels that should remain unmasked. What makes these integer coordinate positions unique? For these particular points, the subsequent interpolation process can utilize a directly available value in at least one direction, meaning no actual interpolation is required along that axis. In these instances, the neighboring points that would normally contribute to the interpolation—especially any anomalous ones—are effectively bypassed. Their associated coefficients are set to zero, which then yields a perfectly valid interpolated value. #### 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 ``` /tmp/ipykernel_3956/3703061842.py:3: RuntimeWarning: invalid value encountered in log 'wo_in_mask_log': np.log(array_out_womask), ![bspline_masking_01.png](grid_resampling_001_files/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 ``` ![bspline_masking_02.png](grid_resampling_001_files/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. ## Boundary condition In this section, we will examine what occurs when interpolating data at the boundaries of the source array. To focus on a specific region of interest, we will crop the source array and store the relevant portion in a new buffer. ```python win_center = np.array((58, 175)) # left eye win_shape = np.array((80, 80)) win = np.array((win_center - win_shape//2, win_center - win_shape//2 + win_shape - 1)).T array_in_cropped = np.empty(win_shape, dtype=array_in.dtype, order='C') array_in_cropped[:] = array_in[win[0][0]:win[0][1]+1, win[1][0]:win[1][1]+1] ``` ![bc_1.png](grid_resampling_001_files/bc_1.png) We generate an identity grid that spans the entire domain of the source image. By utilizing the `resolution` parameter, we apply a zoom factor of 4x. ```python # First create identity grid zoom = 4 if image.ndim == 2: x = np.arange(0, array_in_cropped.shape[0], dtype=grid_dtype) y = np.arange(0, array_in_cropped.shape[1], dtype=grid_dtype) xx, yy = np.meshgrid(x, y) # Lets call the grid resampling array_out_cropped_zoom, array_out_cropped_zoom_mask = array_grid_resampling( interp="cubic", array_in=array_in_cropped, grid_row=yy, grid_col=xx, grid_resolution=(zoom,zoom), array_out=None, win=None, array_in_mask=None, grid_mask=None, array_out_mask=True, nodata_out=None, ) ``` ![bc_2.png](grid_resampling_001_files/bc_2.png) We have just performed a 4x zoom in both directions using a cubic interpolator. This interpolator has a radius of 2, which means that when interpolating near the edges of the image, it needs to access positions outside the source image domain. ```python # Left output edge np.round(array_out_cropped_zoom[4:9, 0:7], 2), array_out_cropped_zoom_mask[4:9, 0:7] ``` (array([[127. , nan, nan, nan, 118. , 102.01, 79.31], [118.77, nan, nan, nan, 116.18, 101.47, 80.01], [107.62, nan, nan, nan, 116.44, 103.69, 83.93], [ 96.67, nan, nan, nan, 115.73, 105.36, 88.05], [ 89. , nan, nan, nan, 111. , 103.18, 89.31]]), array([[1, 0, 0, 0, 1, 1, 1], [1, 0, 0, 0, 1, 1, 1], [1, 0, 0, 0, 1, 1, 1], [1, 0, 0, 0, 1, 1, 1], [1, 0, 0, 0, 1, 1, 1]], dtype=uint8)) As observed, the cubic interpolator returns `NaN` (the value set for nodata_out) for output column indexes 1, 2, and 3. These indexes correspond to source column positions at 0.25, 0.5, and 0.75, respectively. The interpolation requires accessing column values at the source's negative column position -1, which is not available. This invalidates the output mask for these positions, resulting in `NaN` values in the output image. To address these edge cases, the `boundary_condition` parameter (default = `None`) can be set. This parameter handles scenarios where the interpolation requires positions outside the source image domain. The accepted values are a subset of the `numpy.pad` accepted values for the mode parameter, including: - 'edge': Pad with the edge values of the array (repeat boundary values). - 'wrap': Wrap around to the opposite edge (circular/periodic boundary). - 'reflect': Mirror reflection without repeating the edge values. - 'symmetric': Mirror reflection with edge values repeated. It is important to note that when a boundary condition is applied, interpolated values are not systematically invalidated by the output mask. This is because the mask itself is subject to the same boundary condition, ensuring consistency in the handling of edge cases. The `boundary_condition` parameter requires the `standalone` parameter to be set to `True`, which is the default value. We will explore this mode in more detail later. Now, let's apply a 'reflect' boundary condition and examine the effect on the output data. ```python # Lets call the grid resampling with a boundary condition array_out_cropped_zoom_reflect, array_out_cropped_zoom_mask_reflect = array_grid_resampling( interp="cubic", array_in=array_in_cropped, grid_row=yy, grid_col=xx, grid_resolution=(zoom,zoom), array_out=None, win=None, array_in_mask=None, grid_mask=None, array_out_mask=True, nodata_out=None, boundary_condition="reflect", standalone=True, ) ``` ![bc_reflect.png](grid_resampling_001_files/bc_reflect.png) ```python # Left output edge after reflect np.round(array_out_cropped_zoom_reflect[4:9, 0:7], 2), array_out_cropped_zoom_mask_reflect[4:9, 0:7] ``` (array([[127. , 127.4 , 127.31, 124.82, 118. , 102.01, 79.31], [118.77, 119.96, 121.73, 121.37, 116.18, 101.47, 80.01], [107.62, 110.25, 115.36, 118.8 , 116.44, 103.69, 83.93], [ 96.67, 100.52, 108.52, 115.36, 115.73, 105.36, 88.05], [ 89. , 93. , 101.5 , 109.25, 111. , 103.18, 89.31]]), array([[1, 1, 1, 1, 1, 1, 1], [1, 1, 1, 1, 1, 1, 1], [1, 1, 1, 1, 1, 1, 1], [1, 1, 1, 1, 1, 1, 1], [1, 1, 1, 1, 1, 1, 1]], dtype=uint8)) As you can see, all output data is now valid and there are no more `NaN` values. ## The standalone mode Coming in the future... ## Advanced Usage: Pipeline Integration (standalone=False) ### Adjusting the Origin of Input Coordinates The `array_grid_resampling` function includes the parameter `array_in_origin` (default: (0., 0.)) that allows shifting the origin of the input coordinate system. From an implementation standpoint, this shift is applied directly to the computed target grid coordinates. This feature is useful in several scenarios: - Handling different grid conventions without modifying the grid itself. As a reminder, GridR assumes that integer coordinates refer to pixel centers. In other systems, a pixel may be treated as an area, with its center located relatively at (0.5, 0.5). - Working with subregions of a larger input array. In this case, the shift can represent the position of the subregion within the full array’s coordinate system. As an example we will take a previous example using an idendity transformation. ```python # create identity grid if image.ndim == 2: x = np.arange(0, image.shape[0], dtype=grid_dtype) y = np.arange(0, image.shape[1], dtype=grid_dtype) xx, yy = np.meshgrid(x, y) win_center = np.array((58, 175)) # left eye win_shape = np.array((100, 100)) win = np.array((win_center - win_shape//2, win_center - win_shape//2 + win_shape - 1)).T # Call it without specifying array_in_origin array_out_origin_0_0, _ = array_grid_resampling( interp="cubic", array_in=array_in, grid_row=yy, grid_col=xx, grid_resolution=(1,1), array_out=None, win=win, array_in_mask=None, grid_mask=None, array_out_mask=None, nodata_out=None, ) # Call it with specifying array_in_origin and passing standalone=False array_out_origin_30_m10, _ = array_grid_resampling( interp="cubic", array_in=array_in, array_in_origin=(30, -10), # <= shift +30 on rows and -10 on columns grid_row=yy, grid_col=xx, grid_resolution=(1,1), array_out=None, win=win, array_in_mask=None, grid_mask=None, array_out_mask=None, nodata_out=None, standalone=False, ) ``` ![009-array-origin.png](grid_resampling_001_files/009-array-origin.png) As you can observe, the coordinate shift is reflected in the bottom image.