I am trying to create an algorithm that creates smooth color gradient functions, given control points in the red, green, and blue components.
Mathematically, each curve would have a domain [0, 1] and range [0, 255]; every control point would be included in the function, and the interpolation between each point (x0, y0) -> (x1, y1) would not leave the range of [y0, y1]. Moreover, the interpolation should have a continuous first derivative.
Ideally, the result would look something like this visually (ignoring the red lines):
In all of my research, cubic spline interpolation came closest to my requirements, but it did not satisfy the conditions constraining the range of the function. For example, between the 4th and 5th control points (counting from the left) on the bottom function, a cubic spline interpolation would allow the curve to descend below y=0.
As a programmer, I am obviously missing out on some obvious mathematical construct, as these types of curves are seen everywhere in the graphics industry; is anyone aware of a rigorous mathematical definition?
Best Answer
An interpolation by cubic splines has enough degrees of freedom to do what you need. You decompose your domain $[0, 1]$ into intervals $[x_i, x_{i+1}]$, corresponding to your control points. Now, on each interval you define a cubic function $$s_i(x) = a_i x^3+b_i x^2+c_i x+d_i.$$ To compute the coefficients, you have to apply your constraints: $$\begin{align} s_i(x_i) &= y_i \\ s_i(x_{i+1}) &= y_{i+1}\\ s_i'(x_{i+1}) &= s_{i+1}'(x_{i+1})\end{align}$$ This gives you two conditions per cubic for the position and one for the tangent continuity, except for the first and last, where there's no tangent continuity condition.
These aren't enough conditions yet to give a unique solution, as you need 4 conditions for each cubic to determine the 4 coefficients. Ignoring your requirement to stay inside [0, 255] for a second, a common way to interpolate is the "natural" cubic spline, which additionally requires continuity of the second derivative and sets the second derivative to zero at the end-points to give a complete equation system. You would then have to solve that linear system of equations. An alternative are Catmull-Rom splines, which are local and thus don't require solving a linear system.
This does not address your "no overshooting" requirement. You can ignore the second derivatives from the cubic spline and try to work this requirement into the system instead. However, this does not uniquely define your function. One solution to fulfill your requirements is to set the tangents to zero at each control point:
$$\begin{align} s_i(x_i) &= y_i \\ s_i(x_{i+1}) &= y_{i+1}\\ s_i'(x_{i}) &= 0\\ s_i'(x_{i+1}) &= 0\end{align}$$
This will give you sigmoidal shapes for each interval and prevent overshooting. However, it will not look very pretty, it has way too much tension (The natural spline will have minimal tension).
A better solution might be to start with a "good" interpolant like the Catmull-Rom spline and then compute its representation in B-Spline form. B-Splines (like Bezier curves) have the convex hull property. This means, the curve will never be outside of the convex hull of the control points. You can use this to prevent over-shooting. Just ensure that the y-coordinate of your B-Spline control points are in [0, 255]. If the control points you got from computing the spline are outside this range, just set them to 0 or 255 accordingly.
In many vector graphic tools like Inkscape you can draw cubic splines and manipulate the corresponding control points, so you can toy with it there to get a feeling what you need.