How to transform or distort or skew a rectangular image to fit arbitrary quad

I have a rectangular image which need to transform so the image fits an arbitrary convex quad. Easy? Except I have an important constraint, the features on the original image must maintain their relative area after the transformation.

I've tried a perspective transformation, but of course feature's relative area's change significantly with perspective.

The (hand-drawn) image below illustrates the effect I need. The original is on the left and the transformed is on the right. Notice that all 4 features maintain same relative area.

I would be happy with solution that maintained the areas approximately. • It seems to me, more like a Computer Graphics question than a pure DSP. Jan 24 '15 at 19:45

1. Assume that the center of the square maps to the center of the vertexes.
2. The vertices in the output are known. Assume they are $(x_1,y_1),...,(x_4,y_4)$
3. Their center is $( \frac{x_1+x_2+x_3+x_4}{4},\frac{y_1+y_2+y_3+y_4}{4})$
4. Now you have 4 degrees of freedom of where to put the points where the color change. $( w_1 x_1 +(1-w_1) x_2 ,w_1 y_1 +(1-w_1 ) y_2 )$, and so on 1. Areas of each of the 4 parts equal to the areas of two triangles (Can be calculated analytically)

2. All of these together should give you enough constraints to find $w_1,w_2,w_3,w_4$

3. Once you found the weights, you can create a grid based mapping that transforms squares into the required forms.
• if I am understanding your approach correctly, it only works for the sample image not images generally (which is what I am after)
– Ken
Jan 26 '15 at 10:40
• It will work for 4 patches, just like you drew. If you want more patches, it can be extended. If you want the relative size in every area, I suggest seeing @Nikies answer, which is good IMHO. My solution has the advantage of being not approximate but exact. His solution is not exact, but has this property. Jan 26 '15 at 10:51
• Yes I wanted a smooth transformation. @Nikies look like it might be exactly what I need, but It'll take some work before I understand it.
– Ken
Jan 26 '15 at 10:55

Except I have an important constraint, the features on the original image must maintain their relative area after the transformation.

I'm assuming you're looking for a smooth transformation, where all areas (not just the 4 squares in your image) maintain their relative sizes.

IIRC, this constraint is called "fluid registration". (I'm pretty sure I've read about it in this book, but I can't find my copy right now.) The analogy is of course that if your image was "painted" on an incompressible fluid, and you "moved the fluid around", areas would stay constant.

From the incompressible fluid analogy it seems clear that such a mapping should be possible for any convex quad, since you can "squish" an incompressible fluid to any convex quad as long as the area stays the same. You even have degrees of freedom left that you can use to make the transformation as smooth as possible.

Mathematically, the "constant area" constraint means that you're looking for a mapping, where the determinant of the Jacobian of the map is constant. To get the smoothest possible solution, you could e.g. minimize the squared Laplacian. Unfortunately, I don't think there's a closed form solution for this kind of problem, so you'd have to numerical optimize an array of control point to approximate an optimal mapping.

For instance, in Mathematica, you could write something like:

targetPts = {{0.36, 0.62}, {0.34, 0.75}, {0.57, 0.58}, {0.76, 0.92}};

(* simple bilinear mapping *)
grid = Array[
Function[{x,
y}, {(1 - x)*(1 - y), x*(1 - y), (1 - x)*y, x*y}.targetPts], {10,
10}, {{0., 1.}, {0., 1.}}];

offsetVars = Array[var, Dimensions[grid]];

(* correct the bilinear mapping by variable offsets, optimize offsets \
to get constant area *)
mappingPts = grid + offsetVars;

areaQ[pts_] :=
1/2 Det[{pts[] - pts[], pts[] - pts[]}] -
Det[{pts[] - pts[], pts[] - pts[]}]

(* area of each 4x4 control point quad *)
areas = Map[areaQ[Flatten[#, 1]] &,
Partition[mappingPts, {2, 2}, {1, 1}], {2}];

sameArea =
Total[Flatten[{Differences[areas, {0, 1}],
Differences[areas, {1, 0}]}]^2];

smoothness =
Total[Flatten[{Differences[mappingPts, {2, 0}],
Differences[mappingPts, {0, 2}]}]^2];

optimizationProblem = Flatten@{
smoothness,
Flatten[{Differences[areas, {0, 1}],
Differences[areas, {1, 0}]}] == 0,
offsetVars[[1, 1]] == 0,
offsetVars[[1, -1]] == 0,
offsetVars[[-1, 1]] == 0,
offsetVars[[-1, -1]] == 0,
Table[
offsetVars[[1, i]] ==
borderPos[1, i]*(grid[[1, -1]] - grid[[1, 1]]), {i, 2, 9}],
Table[
offsetVars[[-1, i]] ==
borderPos[2, i]*(grid[[-1, -1]] - grid[[-1, 1]]), {i, 2, 9}],
Table[
offsetVars[[i, 1]] ==
borderPos[3, i]*(grid[[-1, 1]] - grid[[1, 1]]), {i, 2, 9}],
Table[
offsetVars[[i, -1]] ==
borderPos[4, i]*(grid[[-1, -1]] - grid[[1, -1]]), {i, 2, 9}]
};

usedVars =
DeleteDuplicates[
Cases[optimizationProblem,
s : (var | borderPos)[i__] :> s, \[Infinity]]];

sol = FindMinimum[
optimizationProblem,
{#, 0} & /@ usedVars];

solGrid = mappingPts /. sol[];
g = Graphics[{Line /@ solGrid, Line /@ Transpose[solGrid], Red,
PointSize[Large], Point[targetPts]}]

This gives nice, smooth mappings for any quad, like this: 