Notes:

The Wikipedia article defines the bilinear blended Coons patch $C\left(s,t\right)$ bounded by the four curves $c_0\left(s\right)$, $c_1\left(s\right)$, $d_0\left(t\right)$, and $d_1\left(t\right)$ as

For SVG mesh gradients the four curves are cubic Bézier curves.

The inkscape wiki states that color blending inside the patch is bilinear. (I have not found where this is specified in Advanced Gradients for SVG ) . So presumably given colors $H_{ij}$ for the four corners we can compute a color within the patch using the formula

My goal is to develop automated software for creating a mesh gradient which best matches a masked source image. There has already been research in this area (Image Vectorization using Optimized Gradients for Ferguson patches) .

A reasonable goal is to minimize the error

which might be an oversimplification since $I\left(x,y\right)$ is sampling an RGB image. (What if we tried to optimize in HSV space?)

The Jian Sun et al paper identifies this as a non-linear least squares problem and the paramters include color parameters $H_{ij}$ , and the space parameters are the control points for $c_0$, $c_1$, $d_0$, and $d_1$ (which I will label $P_{ij}$). They suggest a Levenberg-Marquardt algorithm which depends upon computing a Jacobian (with which I am not very familiar). Maybe it is a Jacobian instead of a gradient because the image has red/green/blue channels.

In any case, I am going to have to figure out how to compute formulas for things like

and evaluate them over the entire image. (I use $\chi$ to represent a chroma element like $r$, $g$, or $b$; and $\gamma$ stands in for a geometric component like $x$ or $y$.)

Computing these partial differentials is relatively straightforward if we iterate over $s$ and $t$.

Now that we have an example of the partial differential with respect to one of the chroma parameters we can ask ourselves: do we really want to calculate the error as the sum over $s$ and $t$?

I feel it more appropriate to calculate it as the sum over $x$ and $y$. That then raises the question: should our algorithm iterate over the relevant pixels of the image and calculate the correct $s$ and $t$ to give each $x$ and $y$? Or should it instead iterate over $s$ and $t$ and scale each element of the sum by $\frac{\partial x}{\partial s}$, $\frac{\partial x}{\partial t}$, $\frac{\partial y}{\partial s}$, and $\frac{\partial y}{\partial t}$ as appropriate? The nature of the formula $I\left(C\left(s,t\right)\right)$ especially concerns me since $⟨x,y⟩=C\left(s,t\right)$ and I’m not absolutely confident of the proper formulation of $\frac{\partial E}{\partial P_{\gamma ij}}$ for a specific $x,y$.

(Now that I have reacquainted myself with integration by substitution for multiple variables) One option is to sum over $s$ and $t$ and use the determinant of the Jacobian to scale the result to $x$ & $y$ space. This gives us a new formula for $E$.

Let’s find out what that Jacobian determinant will look like:

And after all that I need a vacation. $\frac{\partial }{\partial t}C\left(s,t\right)$ will be just as tortuous and have a similar structure, although the structure of the formulas that fall out of $L_d$ and $L_c$ will be swapped.

The source image will be a discrete matrix of RGB triples. $C\left(s,t\right)$ will often not give integer $x$ and $y$ values, so some form of interpolation will be necessary. In images with noise (any photograph, or even a JPEG of a smooth computer painting) it is reasonable to perform some smoothing which might synergize with the interpolation nicely.

Computing geometric component $\frac{\partial E}{\partial P_{\gamma ij}}$ is somewhat more complicated than the chroma components.

That leaves $\frac{\partial }{\partial \gamma }I\left(x,y\right)$. One tactic for computing this would be to use the derivative of the bilinear interpolation of the samples on the boundary of the cell that $x,y$ falls within, but that definition of the derivative is discontinous at cell boundaries (where $x=\lfloor x\rfloor$ or $y=\lfloor y\rfloor$). Upgrading to a bicubic might be an option when I stop being lazy and figure out the related math. Then again, let’s be super-lazy and apply a Gaussian blur to the image!

Let us call the original image $J\left(u,v\right)$ over integer $u$ and $v$ coordinates. Then