WARP_TRI implementation?

[Arcky87]

Any chances to implement warp_tri.pro?
https://www.nv5geospatialsoftware.com/docs/WARP_TRI.html

[alaingdl]

Thanks you for your message.
It is the first time I listen about this procedure !

At first order, in fact this function is a few-lines interface to TRIANGULATE, INTERPOLATE and TRIGRID, but tricks in details :

• Currently : we don't have /quintic keyword in TRIGRID
• We don't have the GRID_TPS() function if WARP_TRI called with /TPS keyword activated

I have a problem : running the IDL example, I got :

GDL> help, WARPED_IMAGE
WARPED_IMAGE    BYTE      = Array[509, 51]
IDL> help, WARPED_IMAGE
WARPED_IMAGE    BYTE      = Array[509, 509]

then we do have another basic problem :(

What I can suggest : we look at the problem of size which should not happen, and we/you test a basic interface without /quintic

/quintic clearly will need more work, and GRID_TPS() also :(

[Arcky87]

• Currently : we don't have /quintic keyword in TRIGRID
• We don't have the GRID_TPS() function if WARP_TRI called with /TPS keyword activated

Looking at my data, I'd say I don't really need /quintic as much as grid_tps(), but I can overcome the lack of both.
The real struggle is this:

GDL> help, WARPED_IMAGE
WARPED_IMAGE    BYTE      = Array[509, 51]
IDL> help, WARPED_IMAGE
WARPED_IMAGE    BYTE      = Array[509, 509]

At the moment, I cannot figure out where exactly the two algorithms diverge.

[alaingdl]

hum hum.
After some tests : TRIANGULATE() gives different values in IDL & GDL, and also the outputs of TRIGRID()
(If I put the IDL values in TRIGRID(), result is as wrong as with the GDL inputs ...)

The field covered by "test_triangulate.pro" is very limited, I don't know id wider tests on TRIANGULATE() & TRIGRID() have bee made ...

[Arcky87]

Unfortunately, I'm not familiar with current GDL tests coverage. I cannot reproduce the examples of TRIGRID() and TRIANGULATE() from NV5 manuals because they're mixing function and procedure syntax which is not implemented in GDL yet.

[alaingdl]

OK, thanks, I go it, we have a bug/feature in line 1228 in triangulation.cpp
if (GS != NULL ) //&& NyDefined)
I can run to the end and the output make sense.

Concerning the modern IDL plotting utilities (isurface & co, plot() functions) since we have close to zero requests on that, we don't move :)

[GillesDuvert]

any chances this can be done with poly_2d ?
look for "warp" occurrences in gdl code

[Arcky87]

any chances this can be done with poly_2d ?
I guess I don't have another option but to rewrite the code with another functions.
The thing is that I have the tricky IDL pipeline for correcting the geometry of spectropolarimetric data. Warp_tri has proven itself as the reliable method for straightening spectra and I've never tested another ways. But as long as it's a combination of existing functionality I can reconstruct the functionality, I suppose

[alaingdl]

I think 99.9% of the code is OK, we have to fix an index problem somewhere, I send to @GillesDuvert a simple test which should help to find the culprit. (code is mostly TRIANGULATE + TRIGRID + INTERPOLATE, I think the problem is in TRIGRID and only for the Y direction). I can put this test case here if someone is ready to have a look.

[alaingdl]

In #2225 I corrected line 1228 in triangulation.cpp
if (GS != NULL ) //&& NyDefined)
It would be great if @GillesDuvert could give an eye on that (why it was like that)

I also add into the Git a very preliminary version of WARP_TRI which is working fine is some cases. I noticed that it is fine for

x0 = [ 10.0, 240.0, 110.0,  30.0]
y0 = [ 20.0,  40.0, 220.0, 170.0]

but it is wrong with

x0=[5,250,480,250]
y0=[250,5,250,480]

I have no idea why :(
(but if I make a shift by 1 on the second, it do work fine ...)

[alaingdl]

@Arcky87 @GillesDuvert
I think I have good news, and obviously bad news !

The not-working case is what is call a "Delaunay co-circular" configuration (the 4 points are on a circle) and may gives degeneracy.
(I know nothing on that topic 😿 )

If I move, eg, y0=[250,5,250,480] to y0=[252,5,250,480] everything is fine 😄

Clearly Gilles wrote a lot of comments in the code (about determinant, signes ..),
we may have to go back to see if we can "fix" something for such specific cases.
We can also add a detector of such configuration.
And I have some ideas on tests

But following my tests, when we are in a "normal" case, our current version of TRIANGULATE + TRIGRID is working fine with the very basic WARP_TRI we have now in PR #2225

[Arcky87]

@alaingdl
Wow, thanks for clarification. Indeed delaunay works great with irregular grid but fails when the grid points are on the same circle OR colinear. In IDL, TRIANGULATE has a TOLERANCE keyword to set the threshold on the test for colinear points, which is useful when you have (almost) regular grid. The default is zero and works pretty well for non-colinear case, but maybe it's not zero by default for "co-circle"? I don't know.

The IDL manual recommends to use value of MAX(X,Y)*10^-12 for double precision as tolerance. Could you test if the x0=[5,250,480,250], y0=[250,5,250,480] will work fine if you add random tolerance to the points?

[alaingdl]

Thanks @Arcky87 to point this out, but I read in the code and the keyword Tolerance is not enforced ...

(line ~95 in triangulation.cpp

    DDouble dtol = isDouble ? 1e-12 : 1e-6;
    //Tol is irrelevant in our implementation, as (s)tripack work with tol=machine precision. 

For sure we will have to look at that.

And also I have new ideas for extending the test suite. (important for me for long term support)

At the moment, are you sure your case is concerned by exceptions ? Aligned or co-circular points ? If not, I just released PR #2225 then please test !

I suggest not closing this issue as long as no clarification for Tol and exceptions, but I have more confident than a week ago : the core of the code works fine (@GillesDuvert 👍 ), we just need to better manage exceptions.

[GillesDuvert]

@Arcky87 @alaingdl Sorry if I was not reactive (I'm in a too complex GDL code update to switch interest). But you seem to handle this just perfect.
For the record, t if the triangulation is on the plane and not on the sphere, GDL uses the "delaunator" engine that triangulates incredibly faster than the code used in IDL. Consequently the triangles order are not those of IDL. In spherical triangulation the 60's stripack fortran code is used (as in IDL: results then should be identical). Both delaunator and stripack use a tolerance in their code, but it is 'hardwired' and the code will need to be slightly hacked to get access and modify this tolerance. I suggest to create an issue with just this topic.