GCP Refinement

[goocreations]

Here is a patch for gdalwarp to refine GCPs by using a -refine_gcps flag. GCPs will be detected as outliers by using a Least Squares Fit and not be used during warping.

[goocreations]

SVN patch for the GCP refinement

[Even Rouault]

I've reviewed the patch (but not tested). My general feeling is that it looks good. However here are my remarks :

• in apps/gdalwarp.cpp, the help "[-refine_gcps]" should mention it needs 2 args. --> [-refine_gcps minimumGcps tolerance].
• "The tolerance is passed to adjust when a GCP will be eliminated." --> it doesn't mention in which units tolerance should be expressed : pixels or target SRS units ?
• I think that the -refince_gcps option only works with polynomial interpolation, not -tps option, right ? If so this should be documented.
• the 2 new transformation options, REFINE_MINIMUM and REFINE_TOLERANCE, should be documented in GDALCreateGenImgProjTransformer2() (alg/gdaltransformer.cpp)
• I'm wondering if the minimumGcps parameter should be validated. Particularly if someone passes a too small value and the interpolation doesn't work. There's perhaps validation at lower stages, but I think you should try with small values like 0,1,2 and see what happens.
• GDALSerializeGCPTransformer() and GDALDeserializeGCPTransformer() should be updated with the 2 new params, so that gdalwarp -of VRT works. Make sure that the deserializer still works with old VRTs that wouldn't have the 2 new options.
• it would be good to stick to the reversed polish notations used generally in GDAL for variable naming. E.g. minimumGcps --> nMinimumGcsp (n for integer) and tolerance -> dfTolerance (df for double float). At least in files where it is used in order to keep consistency. (you'll tell me that alg/gdal_crs.c is an exception to this convention... --> it is third-party code originally). Also, we tend to use 4 spaces for each indentation level.
• in worst_outlier(), residuals array is allocated but never freed.
• I think I got the gist of it, but a bit more comment in worst_outlier() and remove_outliers() (especially in the former) would be appreciated to ease the understanding the maths involved.

We'd appreciate if you could attach a sample dataset to the ticket and a gdalwarp command line example that illustrates the benefit of the patch. Ideally your dataset would be small and freely redistribuable so that we can add a new regression test to the GDAL autotest suite.

As far as the target version, this is really a new feature involving new API, options, parameters, etc.., so 1.8.1 isn't really appropriate. So make sure your patch is based against latest trunk (or at least that it applies without rejects on it). Thanks !

[goocreations]

I've attached an example. I have 51 GCPs, one of these is misplaced. The data shows an extreme case. Using gdalwarp will result in a strangely warped image, but -refine_gcps will remove this misplaced GCP. I'm still busy with the comments/improvements as suggested above, but here is the test data. Usage:

gdalwarp -order 3 test.tif out.tif (will result in inaccurate warping) gdalwarp -order 3 -refine_gcps 32 1.9 test.tif out.tif (will remove the misplaced GCP and result in accurate warping)

​http://www.foxhat.org/gdal/test.tif

[goocreations]

Hi rouault, hope you get this message. regarding the minimumGcps that should be validated: Currently if the user specifies any number lower than needed by the interpolation, the process is started, but then the user is presented with an error: "Failed to compute GCP transform: Not enough points available" (as specified in gdal_crs.c). Is this enough, or should there be a more advanced technique to calculate this before the process is started? I'm currently changing the code so that 0 is the minimum and if minimumGcps is higher than the actual number of GCPs, minimumGcps is set to the number of GCPs.

[Even Rouault]

I think the current error message is fine. I just wanted to make sure there's no crash in this situation.

[goocreations]

Updated GCP refinement

[goocreations]

I've updated the patch, build and tested it against revision 22644. Please note that the parameters for -refine_gcps have changed: the tolerance (in SRS units) is now first and the minimum GCPs the second parameter. If the user omits minimum_gcps, the minimum_gcps according to the current model will be chosen. The minimum GCPs are limited to [0, GCPs available]. If too many GCPs were removed and the model cannot be constructed, a default GDAL error of too few points will be shown to the user. I've updated everything and improved the comments. I hope everything is in order. Please let me know if I forgot something.

[Even Rouault]

Hum, are you sure that it works correctly if no explicit minimumGcps is passed ? In alg/gdaltransformer.cpp, I see at line 1108 that you create GDALCreateGCPRefineTransformer() only if nMinimumGcps >= 0. But if no explicit minimumGcps is passed, nMinimumGcps is set to -1 if REFINE_MINIMUM_GCPS is unset or set to default. I think you would need a bRefine flag in that method, that would be set to TRUE if REFINE_TOLERANCE isn't NULL for example.

I also think that the test if(dfCurrentDifference < 0) at line 1046 of alg/gdal_crs.c will always evaluate to false. dfCurrentDifference is set just at the line before to padfResiduals[nI] which is a square root, thus >= 0.

[goocreations]

Updated GCP refinement

[goocreations]

Sorry about that. Hope everything is in order now.

[Even Rouault]

r22655 /trunk/ (6 files in 4 dirs): gdalwarp: add -refine_gcps option to discard outliers GCPs before warping (#4143)

I'v made the following changes w.r.t your latest patch :

• small fix when -refine_gcps is on the last position of the command line and the minimumGcps is not provided
• fix memory leak of removed psInfo->pasGCPList[nIndex].pszId and psInfo->pasGCPList[nIndex].pszInfo in remove_outliers()
• simplify code of remove_outliers() to avoid the creation of a temporary array in the loop
• I've also added a dump autotest to check this new feature

Due to my changes, please check that it still works however.

[goocreations]

There is a major flaw in my code. I've updated my comments, but never my actual code. So here is a small patch, only changes 3 lines of code. Without this the refinement won't behave as expected.

Sorry about this.

[goocreations]

I've also changed the gdalwarp comment to indicate that minimum_gcps is optional

[Even Rouault]

r22670 /trunk/gdal/ (alg/gdal_crs.c apps/gdalwarp.cpp): GCP refinement: compare directly against tolerance value (gcp_refinement_22669.diff, #4143)

[riaanvddool]

Question asked on irc:

[10:34] <riaanvddool> hi guys, just wondering if there is any info on when 1.9.0 would be released? [10:34] <riaanvddool> i am specifically interested in the new gcp_refine functionality [10:35] <riaanvddool> i would like to write a magazine article for a local (ZA) magazine on our contribution to this functionality, but it doesn't make much sense until the functionality is available in a release [10:36] <riaanvddool> obviously the article would also provide good publicity for GDAL [10:36] <riaanvddool> also: if any other test cases etc are needed, would also like to know about it

[Even Rouault]

GDAL releases occur generally on a yearly basis, in late December or January. If people want to try it before, they can for example try the Windows daily builds (the ones tagged -development) : ​http://www.gisinternals.com/sdk/