Interpolating (hue) angles

In image processing, retouchers may want to apply a saturation boost on specific hues only. Typically, uniform saturation corrections follow a basic linear transfer function $sat_{out} = gain \cdot sat_{in}$, where $gain$ is a real positive constant. To target specific hues, we simply rewrite $sat_{out}(hue) = gain(hue) \cdot sat_{in}$ where $gain$ is then a function. The most common way to declare $gain$ is to provide users with 8 discrete saturation nodes for 8 evenly-spaced hues, and interpolate smoothly between nodes.

The problem is hues are angles, and as far as angles are concerned, $0° \equiv 360° \pm 360°$, meaning they are not your regular 1D euclidean space, but they are actually a ring. This will need extra care not only to connect the curve between 0° and 360°, but also to ensure the continuity of the slope and curvature (first- and second-order derivatives).

How to not do it

I have seen a smart and talented researcher fall into the trap of using regular Catmull-Rom  spline to perform this kind of interpolation. Why is that bad ? Well…

I used here a cubic spline, not the Catmull-Rom, just because it was already available in the Python library Scipy, but the principle holds (although the Catmull-Rom would not overshoot like this). Because I put the hue angle on an euclidean axis, it may not look obviously wrong straight away, but graphs can be deceiving and misleading. Since it’s an angle, let’s display it in polar coordinates.

It should be immediately clear here that we have a problem : hues of 0.1° and 359.9° are expected to get a very similar correction, since they are conceptually so close, but our interpolation doesn’t know it and doesn’t care.

Therefore, we need to enforce a periodic condition, that is, ensure that the interpolated shape will close 0° to 360° as a proper ring and ensure at least $C^1$ continuity (continuity of the first-order derivative), such that the blending of the effects is flawless. None of the usual cubic/Bézier/etc. splines will do it there.

First approach

It makes sense, intuitively, to use trigonometric functions as interpolant here. In fact, Lagrange and Dirichlet have done it already in the 19th century . Since we already know that we will use 8 nodes evenly-spaced, there is even a nice simplification already done for us:

$$\DeclareMathOperator{\sinc}{sinc} y(x) = \sum_{k=1}^{n} y_k \cdot \frac{\sinc(n X / \pi) \cdot \cos(X)}{\sinc(X / \pi)}$$

where $X = (x - x_k) / 2$, $(x_k, y_k)$ are respectively the angular and radial coordinates of the hue nodes, $n$ is the number of nodes (8 for us here), and $\sinc$ is the cardinal sine  defined as $\sinc(x) = \sin(x) / x$ and set to $1$ when $x = 0$ by limit analysis.

Here is the Python code:

 1import numpy as np
 2
 3def get_x_nodes(N):
 4    step_nodes = 2 * np.pi / (N)
 5    result = []
 6    for i in range(N):
 7        result.append(step_nodes * i)
 8    return result
 9
10def trigo_interp(x:np.array, p:np.array):
11    # Interpolate x between N nodes evenly spaced
12    # N needs to be even
13    # p are the user-defined parameters for the radial coordinate of the nodes
14    N = p.size
15    result = np.zeros_like(x)
16    x_nodes = get_x_nodes(N)
17    for i in range(N):
18        X = 0.5 * (x - x_nodes[i])
19        D = np.sinc(N * X / np.pi) * np.cos(X) / np.sinc(X / np.pi) # numpy sinc function is normalized by pi, we don't want that
20        result += D * p[i]
21
22    return result

And let’s display straight away the polar graph :

It’s much better, but there is a problem : the sinc function over- and under-shoot (it’s a known issue and the main flaw of the Lanczos resampling , and the worst is the cusps (inflection points) are not located on the nodes but in-between, which will make that tedious and counter-intuitive to control for the user.

Refined approach

Light and Cheney¹ have published in 1992 a paper about precisely what we want to do here. The problem is this paper is pure, theoretical maths that will be unreadable for most people, and the only numerical example they give is actually to show how it doesn’t work². It took me 6 hours and 3 pages of calculus to digest this pile of hieroglyphs before I finally got some code working. I will try to give you here the primer.

Radial basis functions are basically functions that decay symmetrically as we move away from their center. You may already be familiar with one of them : the Gaussian function . Its center is the mean of the Gaussian distribution. The general form of a radial basis function (RBF) is $f(x) = g(|x - c|)$, where the $|\cdot|$ operator is a scalar norm and $c$ is the coordinate of the center.

The whole idea of the interpolation is to express a set of data as a linear combination of arbitrary continuous base functions. Think of it as making up arbitrary base functions (typically, polynomials), and mix-and-match them. The most simple way is probably the Lagrange interpolation . Solving the interpolation problem therefore falls back to simply determining how much of these functions you need to add in the mix (aka computing the weighting coefficients) such that the continuous interpolation goes through all your discrete data points. The properties of the base functions will decide the properties of the interpolation, such as its monotony and its “adherence” to the nodes to interpolate, and in certain cases, when properly chosen, they can also be used for extrapolation (guessing the value of the function outside of the range where the base data are defined).

RBF interpolation yields an equation like $y(x) = \sum_{k = 0}^{n} \lambda_k f_n(|x - x_k|)$. A linear combination of RBF is called a radiad basis network . RBF are very powerful because the centers $x_k$ can be arbitrarily set at runtime depending on data. But here, we will force them uniformly-spaced on the $[0;2 \pi[$ ring, so $x_k = 2k\pi / N$, where $N$ is the number of nodes and $k \in [0; N[$, and this will allow some nice computational shortcuts and speed-ups.

To compute the $\lambda_k$, we need to write the Gram matrix  $[A]$ of the interpolation, such that we can solve the problem :

$$[A] \cdot \vec{x_k} = \vec{y_k}$$

where $\vec{x_k}$ and $\vec{y_k}$ are respectively the vectors of the angular and radial coordinates four all our control nodes. The square matrix $[A]$ of dimension $N×N$ will be populated at each index $(i, j)$ with the value $f_n(|x_i - x_j|)$, where $x_i$ and $x_j$ are the coordinates of the control nodes. Solving this equation will give us a new vector containing the eigenvalues $\lambda_k$, aka the coefficients for the linear combination of base functions.

Now, if you remember your calculus classes, this will require to diagonalize (or decompose) the matrix, which will fail in 2 cases :

• if the matrix determinant is zero (or, for numerical applications, very close to zero),
• if the matrix condition number is too high.

In addition, since we are interpolating a radius, aka a distance, we can’t output negative values, so we need the matrix to be definite positive.

All these properties of the matrix $[A]$ ultimately depend on the properties of the chosen interpolation base function $f_n$, and the paper derivates very elegantly the necessary conditions on $f_n$ that make the problem solvable all the time, that is, make the $[A]$ matrix robustly diagonalisable.

To shorten, we need $f_n(t) = \sum_{l = 0}^{m \rightarrow \infty} a_l \cos(lt) > 0$, with all $a_l > 0$ and finite, $t = |x - x_k|$, and $f_n(t)$ should be derivable at least twice. In other words, we need to be able to express the function as a cosine series , which is a particular, simplified, case of Fourier series. Which means it needs to have a definite Fourier transform in order to have an expression for $a_l$.

Note that I have taken distances here with the notations of the paper, because the authors express $a_k$ the coefficients of the cosine decomposition, where $k$ denotes the number of terms in the cosine series. Since this overlaps with their previous $k$, the index of the control nodes, I denote it $l$ here³. $m$ is a natural integer arbitrarily chosen depending on the desired precision for the approximation of $f_n$ by the cosine series.

I don’t know for you, but for me this is screaming $\exp(-\alpha x^\beta)$. After some tuning, I converged to $f(t) = \sqrt{\pi} \exp(-x^2)$. From there, getting the cosine transform is easy with Sympy :

1from sympy import *
2from sympy.abc import x, l, a
3cosine_transform(sqrt(pi) * exp(-x**2), x, l)

which returns $a_l = \sqrt{\pi / 2} \cdot \exp(-l^2 / 4) > 0 \forall l < m \in \mathbb{N}$. If you don’t use this function, be sure to understand the constraints lying on $f_n$ to make the interpolation problem solvable all the time.

Note that using a spatial function and computing its cosine transform is only one way of doing it. You can directly declare $a_l$ coefficients the way you like, provided they decay to “almost zero” when $l \rightarrow \infty$. This $a_l$ is actually equivalent to a discrete high-pass filter, so its rate of decay will allow or filter more high frequencies harmonics, allowing more or less oscillations. You could even make the decay tunable, such that $a_l = \exp(-l^2 / c)$, with $c$ user-defined.

Now, at the end of the paper (theorem 26), they add another $\exp$ around our previous $f(t)$ and say “Interpolation at arbitrary distinct nodes is possible with this function”. I would be lying if I said I understand where that comes from, but they seem clever, so I will trust them. Also, I tested without that exponential, and the interpolation actually falls back to the previous trigonometric approach in that case. Being trained in mechanical engineering, I understand this term as a kind of dampener to prevent brutal oscillations. Which is why I love exponentials : they make everything smooth. The shame is $\exp$ also stands for expensive, as far as computers are concerned.

With $f_n$ in our hands, we need to fill the $[A]$ matrix, such that:

$$[A]_{i, j} = \exp\left(\sum_{l = 0}^{m} a_l \cos(l|x_i - x_j|)\right)$$

with $x_i$ and $x_j$, the angular coordinates of the nodes, $i$ being the vertical index, $j$ the horizontal index. Then pseudo-solve through least-squares or whatever to get $\lambda_k$.

Then, the interpolation is expressed by:

$$y(x) = \sum_{k = 0}^{n}\lambda_k \exp\left(\sum_{l = 0}^{m} a_l \cos(l |x - x_k|)\right)$$

As it is customary in maths, the hieroglyphs look worse than the code to make them real:

 1def cosine_coeffs(l, c):
 2    return np.exp(-l**2 / c)
 3
 4def RBF_interpolate(x, p, c):
 5    # Interpolate x between N nodes evenly spaced
 6    # p are the user-defined parameters for the radial coordinate of the nodes
 7    # c is the smoothing parameter
 8
 9    N = p.size
10
11    result = np.zeros_like(x)
12    x_nodes = get_x_nodes(N)
13
14    A = np.zeros((N, N))
15
16    # Number of terms for the cosine series
17    m = int(np.ceil(3 * np.sqrt(c)))
18
19    # Build the A matrix with nodes
20    for i in range(N):
21        for j in range(N):
22            for l in range(m):
23                A[i][j] += cosine_coeffs(l, c) * np.cos(l * np.abs(x_nodes[i] - x_nodes[j]))
24            A[i, j] = np.exp(A[i, j])
25
26    # Solve A * x = y for lambdas
27    lambdas = np.linalg.solve(A, p)
28
29    # Interpolate data for all x
30    for k in range(N):
31        rresult = np.zeros_like(x)
32        for l in range(m):
33            rresult += cosine_coeffs(l, c) * np.cos(l * np.abs(x - x_nodes[k]))
34
35        result += lambdas[k] * np.exp(rresult)
36
37    return result

As you see, the periodic RBF is less prone to over/under-shooting and has the cusps mostly taped to the nodes. This should be less frustrating to control on the user side as well as better blended visually, in terms of the effect.

If we define $a_l(c) = \exp(-l^2 / c)$, which is the un-normalized Gaussian function, with $c$ an user-defined smoothness parameter, we can perform a parametric study of the effect of the decaying rate of $a_l$:

This shows how $c$ relates to the number of harmonics $m$ chosen for the cosine series. Using the well-known properties of the Gaussian distribution , we can deduce that setting $m = 3 \sqrt{c}$ rounded to the nearest upper integer will cover at least 99.6% of the filter’s transmittance × bandwidth while discarding the vanishing coefficients, therefore avoiding useless terms evaluation.

We see that, for $c = 1$, the periodic RBF falls very close to the simple trigonometric approach. As $c$ increases, the decay gets slower, which introduces more harmonic high-frequencies, making the interpolation pattern degrade into a leaf pattern. The interesting property here, however strange the pattern gets when $c > 8$, is that the nodes are all local maxima, which is a desirable property for user-defined control nodes. An optimum seems to be found for $c = \pi$.

Note that the matrix $[A]$ is a Gram matrix, which is a particular Hermitian matrix, which can be solved efficiently by the Choleski decomposition . If you are not using Python and Numpy, I have implemented a Choleski-based matrix solver from scratch in C for darktable  (licensed as GPL v3). Using the properties of the Hermitian matrix, it solves it using half as many computations as the typical Gauss-Jordan elimination method.

Conclusion

The radial basis functions can be fine-tuned for more or less damping, by allowing more or fewer harmonics, through the $c$ smoothing parameter. But they might be gay since they usually have an hard time conforming to straight patterns:

Here, the oscillations have an amplitude of 0.12 %, which is probably acceptable, but the trigonometric interpolation matches the straight line flawlessly. Reducing the magnitude of $c$ cuts more harmonics and better smoothens the interpolation, but at the expense of allowing more overshoot.

All in all, the $a_l$ coefficients can be tweaked to display particular behaviours, like vanishing for even or odd harmonics, allowing only the fundamental frequency, etc. Provided all $a_l$ are positive, the interpolation problem is guaranteed to be solvable.

I find this kind of application wholesome and oddly satisfying, as it bridges clean, elegant, fundamental maths with end-user interface design, in order to provide an intuitive and organically-feeling control over effects applied in the software. As an engineering student, I struggled to pass all my maths exams as all this theory, unconnected to real-life matters, didn’t made any sense for me. I wish I had access to this kind of (fun) application earlier, since starting with a practical problem made unrolling the maths much easier and more interesting, and it’s only by messing around with maths in Python and for practical applications that I began understanding notions I had known for years. Lots of years lost to purely theoretical teaching and dumb exams repetition training…