Opened 13 years ago

Closed 3 years ago

#2652 closed defect (wontfix)

GCP Transform with Order 2

Reported by: rjc1 Owned by: warmerdam
Priority: normal Milestone: closed_because_of_github_migration
Component: Utilities Version: svn-trunk
Severity: major Keywords: GCP transformation warp polynomial
Cc: Even Rouault, gislab

Description

I use GDALGCPTransform to transform points with order 1. This works. I then pass this same transform to GDALWarpKernal to transform an image with order 1. This also works. I repeat the process with order = 2. This again works with GDALGCPTransform, but it does not work with GDALWarpKernal (a hole appears in the image center). Attached Results (zip file) are as follows: GDAL Result 1a.JPG: Screen shot (upper left = source image, lower left image = 1st order result, upper image right = 2nd order result, and Spreadsheet on left shows the GCP points). GDAL GCP 1a.JPG: GCP points are graphed. Upper right "blue" plot = source GCPs, Upper right "purple" plot = target GCPs, Lower left plot shows transformed points and target points perfectly overlap after GDALGCPTransform (the spreadsheet above this plot show the numeric results showing the DX and DY differences are zero). Since the GDALGCPTransform at order 2 seems to work, I assume the problem is in GDALWarpKernal. The source code files (ImageWarpGD is attached). Note that in order to minimize our program footprint, I copied some files, then made deletions to the copied files. This was done to remove file GDALWarpKernal's complete dependency on GDALDataSet and DALWarpOperations files). I suggest this be done in version 1.5.4 and 1.6 since this was the intention according to the GDAL documentation (much better for those not having large images or do not need drivers). The source image is in a separate attached Zip file. This image is provided in two forms JPG and RAW. The JPG can be read by a standard driver. The RAW file is simply a 256x256 floating point array - so no driver is needed as it can be passed directly to GDALWarpKernal through the options structure. A final warpGD.zip file contains all the GDAL source files that were used along with the copied modified ones. As can be seen only a few of the files were needed for this application. Although I assigned this to version 1.5.4, I prefer it gets assigned to 1.6 if a fix is faster.

Attachments (1)

ESDI_GCP.zip (591.7 KB) - added by rjc1 13 years ago.
Support Files for Ticket #2652

Download all attachments as: .zip

Change History (14)

Changed 13 years ago by rjc1

Attachment: ESDI_GCP.zip added

Support Files for Ticket #2652

comment:1 Changed 13 years ago by warmerdam

Cc: Even Rouault added
Component: defaultUtilities
Keywords: warp polynomial added
Milestone: 1.5.4
Priority: highestnormal
Status: newassigned
Version: unspecifiedsvn-trunk

Even,

I'm on the road this week with very limited internet access. Would you be interested in looking into this before I am back on Monday?

comment:2 Changed 13 years ago by Even Rouault

I'm going to look a bit but I won't have much time either as I will be off from Friday to next Wednesday.

As far as this ticket is concerned, there are 2 different things here:

  • a possible bug with order 2 warping. My first tests with the provided data give even worse results that those reported by rjc1. I'm not clear how rjc1 produced that result. That would be nice to have simple steps with gdal_translate and gdalwarp command line utilities. With order 2, when using gdalwarp, I get a completely black image with coordinates with negative y. With order 1 or TPS, the result looks normal. What should be noted is that there are just 6 GCPs provided, and 6 is the minimal number for order 2 interpolation. So this is why we get perfect fit on the GCP, and probably junk outside them. I've looked at the "GDAL Result 1a.JPG" screenshot, and I've already seen similar figures to the top right one on completely different datasets of mine. I'm wondering if it's not another form of the non exact invertness of GCP transform (see #2630)
  • an improvement in the way the code is structured. But IMHO, that's definitely not a 1.5.4 target (1.5 branch is a maintenance branch. No big structural changes), and probably too late for 1.6.0. Or the suggested modifications should be better reported to us. Currently I see a bunch of files in the zip archive. Looks like they are against 1.5 branch and not trunk. And I can't see at first sight what should be done. The ideal format for us would be to get a proper patch made against a trunk checkout. Well, I sound demanding, don't I ? ;-)

For my tests, I've produced a SourceImage?.vrt with the GCP in it with

gdal_translate -of VRT SourceImage.jpg SourceImage.vrt \
-gcp 70.07 55.78 161.22 65.31 \
-gcp 162.59 55.78 202.03 126.53 \
-gcp 212.24 119.05 165.99 196.59 \
-gcp 171.43 205.44 95.93 185.71 \
-gcp 91.84 199.32 48.31 140.13 \
-gcp 38.78 151.02 78.92 60.54

and then

gdalwarp -order 2 SourceImage.vrt dest.tif

comment:3 Changed 13 years ago by gislab

probably related, I reported it half a year ago through the list http://article.gmane.org/gmane.comp.gis.gdal.devel/14911

comment:4 Changed 13 years ago by Even Rouault

It looks like a problem of instability of the algorithm. I've done the following experiment. I've used gdaltransform on the above VRT (with the 6 original GCP), with the TPS algorithm to get the transformed coordinates for (100, 100). The result is 139.17471873512 100.14623756591. I've added this new GCP to the list of 6 and used again order 2. Then the result of gdalwarp looks fine.

I've gone on with that game. And modified the Y transformed component of the additionnal GCP to make it smaller. For example 100 100 139.1747 90 still gives a reasonable result. But with 100 100 139.1747 80 it begins to look very weird with -order 2, but still ok with -tps.

I'm not sure which conclusion to draw on that. Maybe it makes no sense to use just 6 GCPs with order 2 and we must guess an additionnal one ?

gislab, in the thread you're refering too, you mention other software that gives good result whereas GDAL fails. Any chance that some of those software are Open Source ?

comment:5 Changed 13 years ago by gislab

I've tried in ERDAS and ArcGIS, I'll also try GRASS. Here is some more info, if you take R and solve for my points, it solves and predicts correctly, I don't see why GDAL is not doing that, or - it IS doing that, but the image comes up of black. I'm not sure.

-Maxim

######## x = c(448,1761,444,3072,1774,3104) y = c(942,938,3871,928,3870,3858) x2 = c(103,103.5,103,104,103.5,104) y2 = c(53.333333333,53.333333333,52.66666667,53.33333333,52.66666667,52.66666667) an = solve(mat, x2) an bn = solve(mat, y2) bn testpoint = c(500, 1500) xpred = an[1] + an[2]* testpoint[1] + an[3]* testpoint[2] + an[4]* testpoint[1]2 + an[5]* testpoint[1]*testpoint[2] + an[6]*testpoint[2]2 xpred ypred = bn[1] + bn[2]* testpoint[1] + bn[3]* testpoint[2] + bn[4]* testpoint[1]2 + bn[5]* testpoint[1]*testpoint[2] + bn[6]*testpoint[2]2 ypred #102.9993,53.17749

comment:6 Changed 13 years ago by gislab

sorry, trac messed everything up, removed <cr> I hope it is decodable :) Let me know if I need to repost it.

comment:7 Changed 13 years ago by warmerdam

Cc: gislab added

I have done some testing with the provided sample GCPs (per Even's suggestion) and I get the oval output reported by rjc1 (if I use an output file with extents set by an order 1 or tps transformation). It would appear that order 2 polynomial solving is rather sensitive to GCPs in a circle as this does not constrain the middle-to-edge scaling and thus weird things can happen (as is the case here). I'm not sure if this represents a numerical weakness in the gdal_crs.cpp code or if it is inherent in 2nd order polynomials.

I will note that the gdal_crs.cpp code was derived from GRASS so I'm interested in evidence that it produces noticably different results than GRASS.

At this point I'm sliding out of my depth (mathematically speaking) and I'm not sure what to do. It should be possible to try plugging in alternate transformation code if the perspection is that it is the particular implementation here that is the issue.

In the past I had used some numerical recipes code and it would sometimes give noticably different results. However, it was license encumbered in ways that made it impractical to include in GDAL.

PS. Maxim - I was unable to decode the equations you provided though I'm not clear on what they will prove.

comment:8 Changed 13 years ago by gislab

We've finally digged through GDAL/GRASS code to get the part which is calculating transformation matrix. We've checked coefficients matrix calculated using this code and in R - they match perfectly. This means there is no error in calculations of the matrix itself. We'll do some more testing on rasters now and will report here.

Maxim

comment:9 Changed 7 years ago by Jukka Rahkonen

You three (gislab, rouault, warmerdam) did quite a lot of work for solving this issue. Did you ever reach the goal?

comment:10 Changed 7 years ago by Jukka Rahkonen

I repeated the Even's VRT test step by step with GDAL 2.0-dev. The result was a totally black image. Warping with -order 1 of -tps succeeds. The issue has remained the same.

comment:11 Changed 6 years ago by Even Rouault

Severity: blockermajor

Removing blocker criticality

comment:12 Changed 4 years ago by Jukka Rahkonen

Even's VRT test creates a black image still with GDAL 2.3dev.

comment:13 Changed 3 years ago by Even Rouault

Milestone: closed_because_of_github_migration
Resolution: wontfix
Status: assignedclosed

This ticket has been automatically closed because Trac is no longer used for GDAL bug tracking, since the project has migrated to GitHub?. If you believe this ticket is still valid, you may file it to https://github.com/OSGeo/gdal/issues if it is not already reported there.

Note: See TracTickets for help on using tickets.