GIMP LoHalo Pixel Interpolation Mathematics
This article details the mathematical operations occurring within the GEGL (Generic Graphics Library) engine of GIMP when using the LoHalo (Low Halo) interpolation algorithm. LoHalo addresses the ringing artifacts (Gibbs phenomenon) commonly produced by high-order interpolation filters during spatial transformations. It achieves this by adaptively combining high-order interpolation kernels with lower-order, halo-free methods using localized gradient analysis, matrix transformations, and dynamic clamping.
1. Coordinate Mapping and Jacobian Transformation
When an image undergoes scaling, rotation, or affine transformation, GIMP calculates the relationship between the target output grid and the source pixel coordinate space.
- Jacobian Matrix: For any sampling point \((x, y)\), the algorithm computes the \(2 \times 2\) Jacobian matrix of the transformation function: \[J = \begin{bmatrix} \frac{\partial u}{\partial x} & \frac{\partial u}{\partial y} \\ \frac{\partial v}{\partial x} & \frac{\partial v}{\partial y} \end{bmatrix}\]
- Footprint Geometry: The eigenvalues and determinant of \(J\) define an elliptical footprint in the source image space. This footprint determines whether the operation acts primarily as an upsampling (magnification) or downsampling (minification) operation along each principal axis, establishing the width and orientation of the sampling kernel.
2. High-Order Kernel Weighting
For reconstructing continuous intensity surfaces, LoHalo employs higher-order weighting functions—typically cubic convolution or windowed sinc functions (Lanczos).
- Distance Metrics: For each neighbor pixel \(P_i\) located at \((u_i, v_i)\) relative to the continuous target position \((u_0, v_0)\), a normalized Euclidean or elliptical distance \(d_i\) is evaluated: \[d_i = \sqrt{(\Delta u)^2 + (\Delta v)^2}\]
- Polynomial Weight Calculation: Weights \(w_i\) are assigned using a piecewise cubic kernel: \[W(d) = \begin{cases} (a + 2)|d|^3 - (a + 3)|d|^2 + 1 & \text{for } |d| \le 1 \\ a|d|^3 - 5a|d|^2 + 8a|d| - 4a & \text{for } 1 < |d| < 2 \\ 0 & \text{otherwise} \end{cases}\] (where \(a\) typically equals \(-0.5\) or \(-0.75\)).
- Normalization: The preliminary high-order sample \(I_{\text{high}}\) is computed as the normalized dot product of the weights and pixel values: \[I_{\text{high}} = \frac{\sum_{i} w_i \cdot P_i}{\sum_{i} w_i}\]
Negative lobes in these cubic and sinc kernels introduce sharpening, but they also cause negative and excessive values (halos) near sharp transitions.
3. Local Extrema Evaluation
To counter ringing, LoHalo inspects the immediate \(2 \times 2\) pixel neighborhood surrounding the target coordinate \((u_0, v_0)\).
- Bound Constraints: The minimum and maximum pixel values within this local neighborhood are extracted: \[V_{\min} = \min(P_{0,0}, P_{0,1}, P_{1,0}, P_{1,1})\] \[V_{\max} = \max(P_{0,0}, P_{0,1}, P_{1,0}, P_{1,1})\]
These bounds represent the allowable dynamic range for a monotonic reconstruction. Any value outside this range constitutes an overshoot or undershoot (a halo).
4. Bilinear and Area-Based Fallback
Concurrently, LoHalo computes a low-order estimate \(I_{\text{low}}\) using bilinear interpolation:
\[I_{\text{low}} = (1 - s)(1 - t)P_{0,0} + s(1 - t)P_{1,0} + (1 - s)tP_{0,1} + stP_{1,1}\]
where \(s = u_0 - \lfloor u_0 \rfloor\) and \(t = v_0 - \lfloor v_0 \rfloor\). Because bilinear interpolation weights are strictly positive and sum to one, \(I_{\text{low}}\) is mathematically guaranteed to fall entirely within the interval \([V_{\min}, V_{\max}]\), producing zero halos.
5. Adaptive Blending and Clamping
The final pixel value is determined through dynamic blending and clamping functions:
- Hard Clamping Mode: In strict edge conditions, the high-order result is bounded using a simple median/clamping function: \[I_{\text{clamped}} = \max(V_{\min}, \min(V_{\max}, I_{\text{high}}))\]
- Gradient-Weighted Convex Combination: To avoid harsh boundary artifacts from hard clamping, LoHalo calculates a blending parameter \(\alpha \in [0, 1]\) derived from the local contrast (gradient magnitude \(\|\nabla I\|\)). The final interpolated intensity \(I_{\text{final}}\) is determined by a linear blend: \[I_{\text{final}} = \alpha \cdot I_{\text{clamped}} + (1 - \alpha) \cdot I_{\text{low}}\]
This mathematical pipeline ensures that smooth gradients retain the high sharpness and frequency retention of higher-order cubic/sinc sampling, while high-contrast edges automatically collapse to non-oscillating, monotonic bounds to eliminate visual halos.