This implementation is wrong. The inverse-mapping computed in the function MLSD2DWarp is the wrong solution to the problem of backward image interpolation.

The paper of [Schaefer etal. 2006] provides the math for the forward-mapping and it cannot be inverted trivially. One possible solution is to swap the point-sets (p and q) and lose the benefits of the pre-computations, but one can use the fast backward-image-warping for interpolation.

One other solution is to compute the fast forward-mapping as in [Schaefer etal. 2006] and use the slow function "griddata" (or the slightly faster "TriScatteredInterp") instead of "interp2" in order to interpolate the image. In any case, it is not that fast...

Inverting the computed forward lookup-coordinates as done in this code is in any case wrong. You can observe this, if you pull on the handles: one cannot reconstruct the results of the papers. The handles one pulls are not followed by the interpolated image...