Skip to content

quatica.qslst

quatica.qslst

QSLST: Quaternion Special Least Squares with Tikhonov Regularization

This module provides a practical implementation of Algorithm 2 (QSLST) from [1] for quaternion-valued image restoration. It follows the formulation

(A^T A + λ I) X = A^T B,

where A is a real-valued blurring operator, X and B are quaternion images. For convolutional A (Gaussian/motion blur with periodic boundary), we exploit the diagonalization in the Fourier domain to obtain the closed-form Tikhonov solution. For a generic real matrix A, we provide a faithful pinv-based path that implements Algorithm 2 exactly without explicitly forming A(T) = A(A^T A + λ I).

Quaternion representation

We represent a quaternion image Q as a float array of shape (H, W, 4) whose last dimension stores components [q0, q1, q2, q3]. For RGB color, the common choice is q0 = 0, q1 = R, q2 = G, q3 = B. Use rgb_to_quat / quat_to_rgb helpers to convert.

Core API

  • qslst_restore_fft(Bq, psf, lam, boundary="periodic") -> Xq Efficient path for convolutional blur with periodic boundary (BCCB).

  • qslst_restore_matrix(Bq, A_mat, lam) -> Xq Faithful Algorithm 2 using T = A^T A + lam I and pseudo-inverse of T.

  • build_psf_gaussian(radius, sigma) -> psf

  • build_psf_motion(length, angle_deg) -> psf
  • apply_blur_fft(Q, psf, boundary="periodic") -> blurred quaternion
  • add_awgn_snr(Q, snr_db, rng=None) -> noisy quaternion
  • psnr(x, x_ref), relative_error(x, x_ref)

All functions are fully numpy-based.

References

[1] Fei, W., Tang, J., & Shan, M. Quaternion special least squares with Tikhonov regularization method in image restoration. Numerical Algorithms, 1-20. (2025) https://doi.org/10.1007/s11075-025-02187-6

add_awgn_snr(Q, snr_db, rng=None)

Add white Gaussian noise to reach a target SNR (in dB) per quaternion image.

Adds calibrated noise to achieve a specific signal-to-noise ratio. SNR definition: SNR_dB = 10 * log10( ||signal||_F^2 / ||noise||_F^2 )

Parameters:

Q : np.ndarray Clean quaternion image of shape (H, W, 4) snr_db : float Target signal-to-noise ratio in decibels rng : np.random.Generator, optional Random number generator (default: None, uses default_rng())

Returns:

np.ndarray Noisy quaternion image of shape (H, W, 4)

Notes:

The noise variance is computed to achieve the exact target SNR. Returns the original image unchanged if signal power is zero.

apply_blur_fft(Q, psf, boundary='periodic')

Convolve quaternion image with PSF using FFT (per channel).

Applies blurring to a quaternion image by convolving each quaternion component independently with the given point spread function using FFT.

Parameters:

Q : np.ndarray Quaternion image array of shape (H, W, 4) psf : np.ndarray Point spread function kernel of shape (kH, kW) boundary : str, optional Boundary condition, only "periodic" is supported for BCCB (default: "periodic")

Returns:

np.ndarray Blurred quaternion image of shape (H, W, 4)

Notes:

Uses FFT-based convolution for computational efficiency. Each quaternion component is convolved independently, preserving quaternion structure.

build_psf_gaussian(radius, sigma)

Isotropic Gaussian PSF, truncated to a (2*radius+1) window and normalized.

Creates a 2D Gaussian point spread function for image blurring operations. The PSF is truncated to a square window and normalized to unit sum.

Parameters:

radius : int Truncation radius r (PSF size will be 2*radius + 1) sigma : float Standard deviation of the Gaussian distribution

Returns:

np.ndarray Gaussian PSF array of shape (K, K) where K = 2*radius + 1, with sum(psf) = 1

Notes:

Uses the 2D Gaussian formula: exp(-(x² + y²) / (2σ²)) / (2πσ²) The PSF is centered and normalized for convolution operations.

build_psf_motion(length, angle_deg)

Simple linear motion blur PSF of given length and angle, normalized.

Creates a motion blur kernel representing linear movement during exposure. The PSF models uniform motion along a straight line.

Parameters:

length : int Number of pixels in the motion kernel (>= 1) angle_deg : float Motion angle in degrees measured counter-clockwise from +x axis

Returns:

np.ndarray Motion blur PSF array of shape (K, K) with sum = 1, where K is chosen to tightly contain the motion line

Notes:

The kernel size K is automatically determined to contain the motion path. Points along the line are sampled uniformly and rounded to pixel positions.

psnr(x, x_ref, data_range=None)

Peak Signal-to-Noise Ratio for arrays with same shape.

Computes PSNR = 10 * log10(data_range² / MSE) where MSE is the mean squared error. Higher values indicate better image quality/reconstruction.

Parameters:

x : np.ndarray Estimated/reconstructed array x_ref : np.ndarray Reference/ground truth array data_range : float, optional Dynamic range of the data (default: None, uses max(x_ref) - min(x_ref))

Returns:

float PSNR value in decibels, or infinity if MSE is zero

Notes:

Standard metric for image quality assessment. Values typically range from 20-50 dB for natural images, with higher values indicating better quality.

qslst_restore_fft(Bq, psf, lam, boundary='periodic')

QSLST (Algorithm 2) specialized to convolutional A with periodic BC.

Implements the FFT-based Tikhonov solution for quaternion image restoration with convolutional blur operators. Per-channel closed-form solution: X_hat = conj(H_hat) * B_hat / (|H_hat|^2 + lam)

Parameters:

Bq : np.ndarray Observed blurred+noisy quaternion image of shape (H, W, 4) psf : np.ndarray Point spread function of shape (kH, kW) lam : float Tikhonov regularization parameter (>= 0) boundary : str, optional Boundary condition, only "periodic" supported (default: "periodic")

Returns:

np.ndarray Restored quaternion image of shape (H, W, 4)

Notes:

Efficient FFT-based implementation for BCCB (block-circulant with circulant blocks) operators. Each quaternion component is restored independently using the same frequency domain filter.

qslst_restore_matrix(Bq, A_mat, lam)

Faithful Algorithm 2 (QSLST) for a generic real matrix A.

Implements the matrix-based QSLST algorithm for non-convolutional operators

T = A^T A + lam I E = A^T B X = T^+ E (Moore-Penrose pseudoinverse)

Since A is real, T is real; then A(T) = I_4 ⊗ T, and A(T)^+ = I_4 ⊗ T^+. Therefore we can solve per quaternion component independently.

Parameters:

Bq : np.ndarray Observed quaternion image of shape (H, W, 4) A_mat : np.ndarray Real blur matrix of shape (N, N) where N = H*W lam : float Tikhonov regularization parameter (>= 0)

Returns:

np.ndarray Restored quaternion image of shape (H, W, 4)

Notes:

General implementation that works with any real matrix A, not just convolutional operators. Uses pseudoinverse for robust solution even when the system is ill-conditioned.

quat_to_rgb(q, clip=True)

Convert quaternion image to RGB by taking imaginary parts [q1, q2, q3].

Extracts RGB color channels from quaternion imaginary components, with optional intelligent clipping for normalized values.

Parameters:

q : np.ndarray Quaternion image array of shape (H, W, 4) clip : bool, optional Whether to clip RGB to [0, 1] if input appears normalized (default: True)

Returns:

np.ndarray RGB image of shape (H, W, 3)

Notes:

Inverse operation of rgb_to_quat. The clipping heuristic checks if values are in range [-0.5, 1.5] to determine if normalization is appropriate.

relative_error(x, x_ref)

Relative Frobenius error: ||x - x_ref||_F / ||x_ref||_F.

Computes the normalized error between two arrays using the Frobenius norm. This metric is scale-invariant and commonly used for restoration quality assessment.

Parameters:

x : np.ndarray Estimated/reconstructed array x_ref : np.ndarray Reference/ground truth array

Returns:

float Relative error value, or infinity if reference norm is zero

Notes:

Values closer to 0 indicate better reconstruction quality. This metric is particularly useful for comparing restoration algorithms.

rgb_to_quat(rgb, real_part=0.0)

Convert an RGB image to quaternion form (H, W, 4).

Maps RGB color channels to quaternion imaginary parts with configurable real component. Common choice: q0=real_part, q1=R, q2=G, q3=B.

Parameters:

rgb : np.ndarray RGB image array of shape (H, W, 3) with values in [0, 1] or [0, 255] real_part : float, optional Value for quaternion real component q0 (default: 0.0)

Returns:

np.ndarray Quaternion image of shape (H, W, 4) with [q0, q1, q2, q3] = [real_part, R, G, B]

Notes:

This conversion enables processing RGB images as quaternion matrices for quaternion-based image restoration algorithms.

split_quat_channels(q)

Split quaternion image into individual component channels.

Separates a quaternion image into its four scalar components for independent processing or analysis.

Parameters:

q : np.ndarray Quaternion image array of shape (H, W, 4)

Returns:

tuple[np.ndarray, np.ndarray, np.ndarray, np.ndarray] Four channel arrays (q0, q1, q2, q3), each of shape (H, W)

Notes:

Useful for component-wise operations where quaternion channels need to be processed separately.

stack_quat_channels(q0, q1, q2, q3)

Stack scalar channels into quaternion image (H, W, 4).

Combines four separate scalar channel arrays into a single quaternion image for quaternion matrix operations.

Parameters:

q0 : np.ndarray Real component array of shape (H, W) q1 : np.ndarray First imaginary component array of shape (H, W) q2 : np.ndarray Second imaginary component array of shape (H, W) q3 : np.ndarray Third imaginary component array of shape (H, W)

Returns:

np.ndarray Quaternion image of shape (H, W, 4)

Notes:

Inverse operation of split_quat_channels. All input arrays must have the same spatial dimensions (H, W).