[GRASS5] Re: [bug #3166] (grass) bad target region calculations with i.rectify order=3

[Hamish]

> > this bug's URL: http://intevation.de/rt/webrt?serial_num=3166
> > -------------------------------------------------------------------
> > 
> > Subject: bad target region calculations with i.rectify order=3

[1st and 2nd order work fine]

> > 
> > i.rectify (both 6.0.0, 6.1cvs) is creating rasters with incorrect
> > regions and resolutions when doing a 3rd order transform. It works
> > fine with 1st or second order transforms. 'i.rectify -c' (or lower
> > order) is needed to get a good transform.
> > 
> > Starting with unreferenced simple XY map (PNG imported with
> > r.in.gdal) and re-regioned to positive values with r.region, keeping
> > res=1 pixel.
> > 
> > Target location is lat/lon (both wgs84 and none/intl datums fail)
> > 
> > e.g. eastern boundary ends up being +2 degrees east of what it
> > should be, resolution is calculated from e-w/cols so that ends up
> > very coarse... southern boundary is wrong too.
> > 
> > CRS_georef() in imagery/i.rectify/crs.c ?
> > 
> > more info available upon request.
> 
> Hamish,
> 
> you may compare with GDAL:
> 
>  gdal/alg/gdal_crs.c

ok.
 
> - is the algorithm identical?

I think so.

The CRS_georef() fns are essentially the same (GRASS 6/GDAL 1.2.5).

I have cut the top off of both files to CRS_compute_georef_equations(),
diff -wi attached.. No difference there that I can see.

Note that Frank has some comments in his:

 * $Log: gdal_crs.c,v $

[gdal 1.2.6]
+ * Revision 1.8  2004/12/26 16:12:21  fwarmerdam
+ * thin plate spline support now implemented
+ *
[nice..]

 * Revision 1.7  2004/08/11 19:01:56  warmerda
 * default to 2nd order if we have enough points for 3rd order but 0 requested
 *
 * Revision 1.6  2004/08/09 14:38:27  warmerda
 * added serialize/deserialize support for warpoptions and transformers
 *
 * Revision 1.5  2004/03/19 15:22:02  warmerda
 * Fixed double free of padfRasterX. Submitted by Scott Reynolds.
 *
 * Revision 1.4  2002/12/07 17:09:50  warmerda
 * re-enable 3rd order even though unstable
 *
 * Revision 1.3  2002/12/06 17:58:00  warmerda
 * fix a few bugs
 *
 * Revision 1.2  2002/12/05 05:46:08  warmerda
 * fixed linkage problem
 *
 * Revision 1.1  2002/12/05 05:43:23  warmerda
 * New


which tells us 3rd order is not stable...?

If we can get the inputs out, maybe we can try doing the matrix math
in Matlab/Octave and see what answers it gives... (polyfit() fn.)


> - does gdalwarp produce the same (error)

I am starting with a flat PNG(Tiff) image with no georefencing.
I assume you can use something in GDAL to insert GCPs as defined
in the GRASS group/POINTS file, but not the four corners of the 
map - that's where it all goes wrong in i.rectify.

So I don't know how to try gdalwarp on it.
?

 
> If GDAL performs better for above operation, its improvements
> might be worth backporting.

I think GDAL will have the same.
?


Hamish
-------------- next part --------------
--- irec_crs_crop.c	2005-04-21 16:50:03.000000000 +1200
+++ gdal_crs_crop.c	2005-04-21 16:48:10.000000000 +1200
@@ -4,8 +4,11 @@
 */
 /***************************************************************************/
 
-int 
-CRS_compute_georef_equations (struct Control_Points *cp, double E12[], double N12[], double E21[], double N21[], int order)
+static int 
+CRS_compute_georef_equations (struct Control_Points *cp, 
+                                      double E12[], double N12[], 
+                                      double E21[], double N21[], 
+                                      int order)
   {
   double *tempptr;
   int status;
@@ -77,22 +80,22 @@
 
   /* INITIALIZE MATRIX */
 
-  m.v = (DOUBLE *)G_calloc(m.n*m.n,sizeof(DOUBLE));
+    m.v = (double *)CPLCalloc(m.n*m.n,sizeof(double));
   if(m.v == NULL)
     {
     return(MMEMERR);
     }
-  a = (DOUBLE *)G_calloc(m.n,sizeof(DOUBLE));
+    a = (double *)CPLCalloc(m.n,sizeof(double));
   if(a == NULL)
     {
-    G_free((char *)m.v);
+        CPLFree((char *)m.v);
     return(MMEMERR);
     }
-  b = (DOUBLE *)G_calloc(m.n,sizeof(DOUBLE));
+    b = (double *)CPLCalloc(m.n,sizeof(double));
   if(b == NULL)
     {
-    G_free((char *)m.v);
-    G_free((char *)a);
+        CPLFree((char *)m.v);
+        CPLFree((char *)a);
     return(MMEMERR);
     }
 
@@ -101,9 +104,9 @@
   else
     status = calcls(cp,&m,a,b,E,N);
 
-  G_free((char *)m.v);
-  G_free((char *)a);
-  G_free((char *)b);
+    CPLFree((char *)m.v);
+    CPLFree((char *)a);
+    CPLFree((char *)b);
 
   return(status);
   }