Source code for pylops.signalprocessing.bilinear

__all__ = ["Bilinear"]

import logging

import numpy as np
import numpy.typing as npt

from pylops import LinearOperator
from pylops.utils.backend import get_add_at, get_array_module, to_numpy
from pylops.utils.decorators import reshaped
from pylops.utils.typing import DTypeLike, InputDimsLike, IntNDArray, NDArray

logging.basicConfig(format="%(levelname)s: %(message)s", level=logging.WARNING)


def _checkunique(iava: npt.ArrayLike) -> None:
    """Check that vector as only unique values"""
    _, count = np.unique(iava, axis=1, return_counts=True)
    if np.any(count > 1):
        raise ValueError("Repeated values in iava array")


[docs]class Bilinear(LinearOperator): r"""Bilinear interpolation operator. Apply bilinear interpolation onto fractionary positions ``iava`` along the first two axes of a n-dimensional array. .. note:: The vector ``iava`` should contain unique pais. If the same pair is repeated twice an error will be raised. Parameters ---------- iava : :obj:`list` or :obj:`numpy.ndarray` Array of size :math:`[2 \times n_\text{ava}]` containing pairs of floating indices of locations of available samples for interpolation. dims : :obj:`list` Number of samples for each dimension forceflat : :obj:`bool`, optional .. versionadded:: 2.2.0 Force an array to be flattened after rmatvec. Note that this is only required when `len(dims)=2`, otherwise pylops will detect whether to return a 1d or nd array. dtype : :obj:`str`, optional Type of elements in input array. name : :obj:`str`, optional .. versionadded:: 2.0.0 Name of operator (to be used by :func:`pylops.utils.describe.describe`) Attributes ---------- shape : :obj:`tuple` Operator shape explicit : :obj:`bool` Operator contains a matrix that can be solved explicitly (``True``) or not (``False``) Raises ------ ValueError If the vector ``iava`` contains repeated values. Notes ----- Bilinear interpolation of a subset of :math:`N` values at locations ``iava`` from an input n-dimensional vector :math:`\mathbf{x}` of size :math:`[m_1 \times m_2 \times ... \times m_{ndim}]` can be expressed as: .. math:: y_{\mathbf{i}} = (1-w^0_{i}) (1-w^1_{i}) x_{l^{l,0}_i, l^{l,1}_i} + w^0_{i} (1-w^1_{i}) x_{l^{r,0}_i, l^{l,1}_i} + (1-w^0_{i}) w^1_{i} x_{l^{l,0}_i, l^{r,1}_i} + w^0_{i} w^1_{i} x_{l^{r,0}_i, l^{r,1}_i} \quad \forall i=1,2,\ldots,M where :math:`\mathbf{l^{l,0}}=[\lfloor l_1^0 \rfloor, \lfloor l_2^0 \rfloor, ..., \lfloor l_N^0 \rfloor]`, :math:`\mathbf{l^{l,1}}=[\lfloor l_1^1 \rfloor, \lfloor l_2^1 \rfloor, ..., \lfloor l_N^1 \rfloor]`, :math:`\mathbf{l^{r,0}}=[\lfloor l_1^0 \rfloor + 1, \lfloor l_2^0 \rfloor + 1, ..., \lfloor l_N^0 \rfloor + 1]`, :math:`\mathbf{l^{r,1}}=[\lfloor l_1^1 \rfloor + 1, \lfloor l_2^1 \rfloor + 1, ..., \lfloor l_N^1 \rfloor + 1]`, are vectors containing the indices of the original array at which samples are taken, and :math:`\mathbf{w^j}=[l_1^i - \lfloor l_1^i \rfloor, l_2^i - \lfloor l_2^i \rfloor, ..., l_N^i - \lfloor l_N^i \rfloor]` (:math:`\forall j=0,1`) are the bilinear interpolation weights. """ def __init__( self, iava: IntNDArray, dims: InputDimsLike, forceflat: bool = None, dtype: DTypeLike = "float64", name: str = "B", ) -> None: # define dimension of data ndims = len(dims) dimsd = [len(iava[1])] + list(dims[2:]) # check if forceflat is needed and set it back to None otherwise if ndims > 2: if forceflat is not None: logging.warning( f"setting forceflat=None since len(dims)={len(dims)}>2. " f"PyLops will automatically detect whether to return " f"a 1d or nd array based on the shape of the input" f"array." ) forceflat = None super().__init__( dtype=np.dtype(dtype), dims=dims, dimsd=dimsd, forceflat=forceflat, name=name, ) ncp = get_array_module(iava) # check non-unique pairs (works only with numpy arrays) _checkunique(to_numpy(iava)) # find indices and weights self.iava_t = ncp.floor(iava[0]).astype(int) self.iava_b = self.iava_t + 1 self.weights_tb = iava[0] - self.iava_t self.iava_l = ncp.floor(iava[1]).astype(int) self.iava_r = self.iava_l + 1 self.weights_lr = iava[1] - self.iava_l # expand dims to weights for nd-arrays if ndims > 2: for _ in range(ndims - 2): self.weights_tb = ncp.expand_dims(self.weights_tb, axis=-1) self.weights_lr = ncp.expand_dims(self.weights_lr, axis=-1) @reshaped def _matvec(self, x: NDArray) -> NDArray: return ( x[self.iava_t, self.iava_l] * (1 - self.weights_tb) * (1 - self.weights_lr) + x[self.iava_t, self.iava_r] * (1 - self.weights_tb) * self.weights_lr + x[self.iava_b, self.iava_l] * self.weights_tb * (1 - self.weights_lr) + x[self.iava_b, self.iava_r] * self.weights_tb * self.weights_lr ) @reshaped def _rmatvec(self, x: NDArray) -> NDArray: ncp = get_array_module(x) ncp_add_at = get_add_at(x) y = ncp.zeros(self.dims, dtype=self.dtype) ncp_add_at( y, tuple([self.iava_t, self.iava_l]), x * (1 - self.weights_tb) * (1 - self.weights_lr), ) ncp_add_at( y, tuple([self.iava_t, self.iava_r]), x * (1 - self.weights_tb) * self.weights_lr, ) ncp_add_at( y, tuple([self.iava_b, self.iava_l]), x * self.weights_tb * (1 - self.weights_lr), ) ncp_add_at( y, tuple([self.iava_b, self.iava_r]), x * self.weights_tb * self.weights_lr ) return y