Given two sets of points in 2D points1
and points2
inside a rectangle, I would like to find a continuous map f
defined on the rectangle that sends points1
to points2
. The set of points are ordered, which makes it considerably easier (points1[[i]]
should be sent to points2[[i]]
).
{points1, points2, trianglesIndices} = Import["https://pastebin.com/raw/n8teGNSM"] // ToExpression;
GraphicsRow[{ListPlot[points1], ListPlot[points2]}]
First, I made a Delaunay mesh and find the affine transformations on each pair of triangles:
triangles1 = points1[[#]] & /@ trianglesOrder;
triangles2 = points2[[#]] & /@ trianglesOrder;
Then, using inPolyQ
from this answer to defined a piecewise function, mapping each triangle to its image:
inPolyQ[poly_, pt_] := Graphics`PolygonUtils`PointWindingNumber[poly, pt] =!= 0
mats = Table[FindGeometricTransform[triangles2[[i]], triangles1[[i]],
TransformationClass -> "Affine"][[2, 1]], {i, Length@triangles1}];
f[{x_, y_}, coef_] := Piecewise[Table[{TransformationFunction[
coef*mats[[i]] + (1 - coef)*IdentityMatrix[3]][{x, y}],
inPolyQ[triangles1[[i]], {x, y}]}, {i, Length@triangles1}], {x, y}]
We can check f
maps points1
to points2
, with the parameter coef
ponderating inbetween:
coef = 0;
Norm[func[#, coef] & /@ points1 - points1]
coef = 1;
Norm[func[#, coef] & /@ points1 - points2]
(* 0. *)
(* 7.50594*10^-12 *)
The problem is, it is very very slow when used with ImageTransform
(code provided below).
I tried other approaches: transforming each triangle independently and assembling everything in the end, but I got to defaults because of the borders of the original image being giving spurious diagonals. I also tried making masks, making an affine transformation of the whoel original image for each triangle and applying the masks, then combining together. It almost works, but not totally, possibly because of some small margins that interfere.
Anyway, I believe I'm missing some simpler approaches, and that's my question.
Here is an example image if need be, consisting of a grid and the triangles whose vertices constitute points1
:
width = 300;
height = 400;
start = Graphics[{White, EdgeForm[Thick], Triangle /@ triangles1,
Thin, Black, Line@Table[{{1, h}, {300, h}}, {h, 1, 400, 10}],
Line@Table[{{w, 1}, {w, height}}, {w, 1, width, 10}]},
PlotRangePadding -> 0, ImageSize -> {width, height}]
Show[start, Graphics[{Red, PointSize[0.02], Point[points1]},
PlotRangePadding -> None]]
Then, the final image can be obtained with ImageTransform
but takes several minutes (!):
funcInv[{x_, y_}, coef_] := Piecewise[ Table[{TransformationFunction[
coef*Inverse[mats[[i]]] + (1 - coef)*IdentityMatrix[3]][{x, y}],
inPolyQ[triangles2[[i]], {x, y}]}, {i, Length@triangles2}], {x, y}]
(* very slow, but works *)
coef = 1;
final = ImageTransformation[start, funcInv[#, coef] &, DataRange -> Full]
Show[final, Graphics[{Blue, PointSize[0.02], Point[points2]},
PlotRangePadding -> None]]
Is there a faster approach?
FindGeometricTransform
does precisely this. $\endgroup$FindGeometricTransform
, that I am using on each triangle, limited to linear (or affine) mappings? $\endgroup$