Opened 6 years ago

Closed 6 years ago

Last modified 5 years ago

#4143 closed enhancement (fixed)

GCP Refinement

Reported by: goocreations Owned by: warmerdam
Priority: normal Milestone: 1.9.0
Component: default Version: svn-trunk
Severity: normal Keywords: gcp refine refinement
Cc:

Description

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.

Attachments (4)

gcp_refinement.diff (14.3 KB) - added by goocreations 6 years ago.
SVN patch for the GCP refinement
gcp_refinement_22644.diff (24.7 KB) - added by goocreations 6 years ago.
Updated GCP refinement
gcp_refinement_22652.diff (24.9 KB) - added by goocreations 6 years ago.
Updated GCP refinement
gcp_refinement_22669.diff (2.4 KB) - added by goocreations 6 years ago.

Download all attachments as: .zip

Change History (19)

Changed 6 years ago by goocreations

SVN patch for the GCP refinement

comment:1 Changed 6 years ago by 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 !

comment:2 Changed 6 years ago by 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

comment:3 Changed 6 years ago by goocreations

  • Milestone changed from 1.8.1 to 1.9.0

comment:4 Changed 6 years ago by 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.

comment:5 Changed 6 years ago by rouault

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

Changed 6 years ago by goocreations

Updated GCP refinement

comment:6 Changed 6 years ago by 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.

comment:7 Changed 6 years ago by 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.

Changed 6 years ago by goocreations

Updated GCP refinement

comment:8 Changed 6 years ago by goocreations

Sorry about that. Hope everything is in order now.

comment:9 Changed 6 years ago by rouault

  • Resolution set to fixed
  • Status changed from new to closed

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.

comment:10 Changed 6 years ago by 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.

Changed 6 years ago by goocreations

comment:11 Changed 6 years ago by goocreations

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

comment:12 Changed 6 years ago by riaanvddool

  • Resolution fixed deleted
  • Status changed from closed to reopened

comment:13 Changed 6 years ago by rouault

  • Resolution set to fixed
  • Status changed from reopened to closed

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

comment:14 Changed 5 years ago by 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

comment:15 Changed 5 years ago by 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/

Note: See TracTickets for help on using tickets.