2d signal module ---------------- .. automodule:: lazylinop.signal2d This module provides lazy linear operators associated to fast linear transforms commonly used in image processing. List of transforms ~~~~~~~~~~~~~~~~~~ 1. :meth:`lazylinop.signal2d.convolve2d` Convolution with a given (2D) filter. 2. :meth:`lazylinop.signal2d.dwt2d` Discrete Wavelet Transform 3. :meth:`lazylinop.signal2d.idwt2d` Inverse of Discrete Wavelet transform 4. :meth:`lazylinop.signal2d.dft2d` Discrete Fourier Transform 5. :meth:`lazylinop.signal2d.dct2d` Discrete Cosine Transform of types I to IV 6. :meth:`lazylinop.signal2d.dst2d` Discrete Sine Transform of types I to IV 7. :meth:`lazylinop.signal2d.mdct2d` Modified DCT Input (and output) vector shape :octicon:`alert-fill;1em;sd-text-danger` ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ Applying any of these linear operators ``L`` to a single input image ``X`` provided as a 2D-array (or a batch of such 2D-arrays) requires prior *flattening* of images. A typical usage of 2d signal module would be: .. code-block:: python >>> shape2d = X.shape[0], X.shape[1] # X is an image or a batch of images >>> L = foo2d(shape2d, ...) >>> y = L @ colvec(X) # Apply L to the flattened version of x >>> Y = uncolvec(y, shape2d) # Un-flattened version of y This requirement is due to the need to ensure compatibility with the LazyLinOp API, where the result of ``L @ V`` is specified as the concatenation of the column vectors obtained by applying ``L`` to each column of ``V``. Do not try to compute ``y = L @ X``, it would not give the desired result (and would even simply generate an error due to non-matching dimensions of the columns of X with the expected input dimension of the operator ``L``)! However, note that `signal2d` functions usually take an `in_shape` as first argument, where :math:`\text{in_shape}=\left(M,~N\right)` the two dimensions of the input image(s). As a concrete example: - if ``X`` is a single 2D image, it is needed to perform ``y = L @ X.ravel()``. - if instead ``X`` is a batch of 2D images, where the first two dimension are the 2D image dimension and the third is the batch dimension (as usual with LazyLinOp API), it would need: ``y = L @ X.reshape(-1, X.shape[2])`` For convenience, LazyLinOp provides ``colvec(X)`` utility function (:py:func:`colvec`) that performs the required reshaping in both -single or batch of- image situations (internally it implements ``X.reshape(-1, *(X.shape[2:]))``. :octicon:`alert-fill;1em;sd-text-danger` For operators where it can be expected than the output returns an transformed image (e.g. ``convolve2d``, ``padder2d``,...), the output `y` is also *flattened*. For the same reason, ``y`` requires *reshaping* to be manipulated as an image (or an image batch). We also provide an ``uncolvec(y, in_shape)`` utility function (:py:func:`uncolvec`) for that purpose. Padding or cropping ~~~~~~~~~~~~~~~~~~~ Traditional implementations of image processing transforms offer the possibility to either crop or zero-pad the analyzed image. We choose not to include this in our implementation, as mimicing this feature is simple with the generic lazylinop interface. .. [FIXME: the exception is the DWT2D, do we want to keep this exception? Same question in the 1d case in fact] Consider for the example ``F2 = dft2d(in_shape=(N, N))``, the operator associated to the 2D DFT matrix on :math:`N\times N` images. To apply it to an image ``X`` of size :math:`N\times N`, we can define ``G = eye(N, n)``, ``G2 = kron(G, G)`` and observe that ``F2 @ G2`` does exactly what we need: - padding: if :math:`nN`, similarly ``G2 @ colvec(X)`` is a (flattened) cropped version of ``X`` of size N x N, and ``H = F2 @ G2`` is again exactly what you need. .. [FIXME: to check !] For convenience, we provide an operator ``padder2d(in_shape, width, mode)`` to pad a 2D array (provided "of course" in a flattened form) with specific boundary conditions ``zero``, ``periodic``, ``symmetric``, ``antisymmetric`` and ``reflect``. ``padder2d()`` is obtained from the Kronecker product ``kron(L1, L2)`` of two operators ``L1 = padder(in_shape[0], width[0], mode)`` and ``L2 = padder(in_shape[1], width[1], mode)``. .. warning:: Do not confuse ``padder2d`` with ``pad`` (see :ref:`basicops-construc` for more details). Inverse transforms ~~~~~~~~~~~~~~~~~~ As in the signal module :doc:`api_signal`, most of our 2D transforms are *orthonormal*, and their inverse is thus their adjoint, e.g., with ``F`` the DFT, the adjoint ``F.H`` is the lazy linear operator associated to the inverse transform (appropriately taking into account all parameters of the original transform). The main exceptions are :py:func:`convolve2d`, :py:func:`mdct2d` (this operator satisfies ``L @ L.T = L @ L.H = Id`` but not ``L.T @ L = Id``) and :py:func:`dwt2d` (see its documentation for details on parameters ensuring that ``L = dwt2d(...)`` satisfies ``L @ L.T = L.T @ L = Id`` or only ``L.T @ L = Id`` when ``L`` is rectangular). .. [FIXME idwt2d??] .. note:: This is because of most our 2D transforms are Kronecker product :math:`K=L_1\otimes L_2` of two 1D transforms :math:`L_1` and :math:`L_2`. The matrix product of the adjoint :math:`K^H` with :math:`K` is given by :math:`K^HK=\left(L_1^H\otimes L_2^H\right)\left(L_1\otimes L_2\right)` and we have :math:`K^HK=\left(L_1^HL\right)\otimes\left(L_2^HL_2\right)` (mixed-product property). If :math:`L_1` and :math:`L_2` are orthogonal, :math:`K^HK` is equal to the identity matrix. Therefore, the 2D transform associated to :math:`K` is orthogonal. Non-orthonormal versions of these transforms correspond to alternative normalizations described below. Various normalizations ~~~~~~~~~~~~~~~~~~~~~~ Traditional implementations of image processing transforms offer various normalization (e.g., with division by :math:`N`, :math:`\sqrt{N}`, or no division). They can all be mimicked (if really needed) by pre- or post-composing the transform ``F``. As an illustration, SciPy’s :py:func:`fft2d` and :py:func:`ifft2d` with the default normalization is mimicked as follows. .. code-block:: python >>> import numpy as np >>> from scipy.fft import fft2 as sp_fft2d >>> from scipy.fft import ifft2 as sp_ifft2d >>> from lazylinop.signal2d import dft2d as lz_dft2d >>> from lazylinop.signal2d import colvec >>> from lazylinop.signal2d import uncolvec >>> N = 32 >>> X = np.random.randn(N, N) >>> F2 = lz_dft2d((N, N)) >>> scale = N >>> y = scale * F2 @ colvec(X) >>> Y = sp_fft2d(X) >>> np.allclose(colvec(y, (N, N)), Y) >>> True >>> x_ = F2.H @ y / scale >>> X_ = sp_ifft2d(Y) >>> np.allclose(colvec(x_, (N, N)), X_) >>> True To mimick SciPy's DCT/DST called with ``orthogonalize=True``, the same trick holds where ``scale`` depends on the transform's type (I,II,III,IV) and the choice of ``norm`` (``'backward'`` or ``'forward'``; NB: ``scale = 1`` if ``norm = 'ortho'``). Mimicking the DCT/DST with ``orthogonalize=False`` requires pre- and/or post-composing by diagonal operators that depend on the type. For example, the default DCT-II behavior (``norm = 'ortho'``): .. code-block:: python >>> from lazylinop.signal2d import dst2d as lz_dst2d >>> from scipy.fft import dstn as sp_dstn >>> import numpy as np >>> M, N = 32, 32 >>> X = np.random.randn(M, N) >>> L = lz_dst2d(X.shape) >>> Y = L @ colvec(X) >>> from lazylinop.basicops import diag >>> v = np.full(N, 1.0) >>> v[-1] = np.sqrt(2.0) >>> D = diag(v) >>> Z = sp_dstn(X, 2, (M, N), (0, 1), 'ortho', False, 1, orthogonalize=False) >>> np.allclose(D @ Y.reshape(M, N) @ D, Z) True Transforms ~~~~~~~~~~ .. autofunction:: lazylinop.signal2d.convolve2d .. autofunction:: lazylinop.signal2d.dwt2d .. autofunction:: lazylinop.signal2d.idwt2d .. autofunction:: lazylinop.signal2d.dft2d .. autofunction:: lazylinop.signal2d.fft2d .. autofunction:: lazylinop.signal2d.dct2d .. autofunction:: lazylinop.signal2d.dst2d .. autofunction:: lazylinop.signal2d.mdct2d .. autofunction:: lazylinop.signal2d.padder2d Utility functions ~~~~~~~~~~~~~~~~~ 1. :meth:`lazylinop.flatten` 2. :meth:`lazylinop.unflatten` 3. :meth:`lazylinop.signal2d.dwt2d_to_pywt_coeffs` 4. :meth:`lazylinop.signal2d.dwt2d_coeffs_shapes` 5. :meth:`lazylinop.colvec` 6. :meth:`lazylinop.uncolvec` .. autofunction:: lazylinop.signal2d.flatten .. autofunction:: lazylinop.signal2d.unflatten .. autofunction:: lazylinop.signal2d.dwt2d_to_pywt_coeffs .. autofunction:: lazylinop.signal2d.dwt2d_coeffs_shapes .. autofunction:: lazylinop.signal2d.colvec .. autofunction:: lazylinop.signal2d.uncolvec