Find a continuous map between two sets of 2D points

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?