1 : /******************************************************************************
2 : * $Id: gdal_crs.c 16614 2009-03-17 23:46:39Z rouault $
3 : *
4 : * Project: Mapinfo Image Warper
5 : * Purpose: Implemention of the GDALTransformer wrapper around CRS.C functions
6 : * to build a polynomial transformation based on ground control
7 : * points.
8 : * Author: Frank Warmerdam, warmerdam@pobox.com
9 : *
10 : ***************************************************************************
11 :
12 : CRS.C - Center for Remote Sensing rectification routines
13 :
14 : Written By: Brian J. Buckley
15 :
16 : At: The Center for Remote Sensing
17 : Michigan State University
18 : 302 Berkey Hall
19 : East Lansing, MI 48824
20 : (517)353-7195
21 :
22 : Written: 12/19/91
23 :
24 : Last Update: 12/26/91 Brian J. Buckley
25 : Last Update: 1/24/92 Brian J. Buckley
26 : Added printout of trnfile. Triggered by BDEBUG.
27 : Last Update: 1/27/92 Brian J. Buckley
28 : Fixed bug so that only the active control points were used.
29 :
30 :
31 : Copyright (c) 1992, Michigan State University
32 :
33 : Permission is hereby granted, free of charge, to any person obtaining a
34 : copy of this software and associated documentation files (the "Software"),
35 : to deal in the Software without restriction, including without limitation
36 : the rights to use, copy, modify, merge, publish, distribute, sublicense,
37 : and/or sell copies of the Software, and to permit persons to whom the
38 : Software is furnished to do so, subject to the following conditions:
39 :
40 : The above copyright notice and this permission notice shall be included
41 : in all copies or substantial portions of the Software.
42 :
43 : THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
44 : IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
45 : FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
46 : THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
47 : LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
48 : FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
49 : DEALINGS IN THE SOFTWARE.
50 :
51 : ****************************************************************************/
52 :
53 : #include "gdal_alg.h"
54 : #include "cpl_conv.h"
55 : #include "cpl_minixml.h"
56 : #include "cpl_string.h"
57 :
58 : CPL_CVSID("$Id: gdal_crs.c 16614 2009-03-17 23:46:39Z rouault $");
59 :
60 : #define MAXORDER 3
61 :
62 : struct Control_Points
63 : {
64 : int count;
65 : double *e1;
66 : double *n1;
67 : double *e2;
68 : double *n2;
69 : int *status;
70 : };
71 :
72 : CPL_C_START
73 : CPLXMLNode *GDALSerializeGCPTransformer( void *pTransformArg );
74 : void *GDALDeserializeGCPTransformer( CPLXMLNode *psTree );
75 : CPL_C_END
76 :
77 : /* crs.c */
78 : static int CRS_georef(double, double, double *, double *,
79 : double [], double [], int);
80 : static int CRS_compute_georef_equations(struct Control_Points *,
81 : double [], double [], double [], double [], int);
82 :
83 :
84 : static char *CRS_error_message[] = {
85 : "Failed to compute GCP transform: Not enough points available",
86 : "Failed to compute GCP transform: Transform is not solvable",
87 : "Failed to compute GCP transform: Not enough memory"
88 : "Failed to compute GCP transform: Parameter error",
89 : "Failed to compute GCP transform: Internal error"
90 : };
91 :
92 : typedef struct
93 : {
94 : GDALTransformerInfo sTI;
95 :
96 : double adfToGeoX[20];
97 : double adfToGeoY[20];
98 :
99 : double adfFromGeoX[20];
100 : double adfFromGeoY[20];
101 :
102 : int nOrder;
103 : int bReversed;
104 :
105 : int nGCPCount;
106 : GDAL_GCP *pasGCPList;
107 :
108 : } GCPTransformInfo;
109 :
110 :
111 : /************************************************************************/
112 : /* GDALCreateGCPTransformer() */
113 : /************************************************************************/
114 :
115 : /**
116 : * Create GCP based polynomial transformer.
117 : *
118 : * Computes least squares fit polynomials from a provided set of GCPs,
119 : * and stores the coefficients for later transformation of points between
120 : * pixel/line and georeferenced coordinates.
121 : *
122 : * The return value should be used as a TransformArg in combination with
123 : * the transformation function GDALGCPTransform which fits the
124 : * GDALTransformerFunc signature. The returned transform argument should
125 : * be deallocated with GDALDestroyGCPTransformer when no longer needed.
126 : *
127 : * This function may fail (returning NULL) if the provided set of GCPs
128 : * are inadequate for the requested order, the determinate is zero or they
129 : * are otherwise "ill conditioned".
130 : *
131 : * Note that 2nd order requires at least 6 GCPs, and 3rd order requires at
132 : * least 10 gcps. If nReqOrder is 0 the highest order possible with the
133 : * provided gcp count will be used.
134 : *
135 : * @param nGCPCount the number of GCPs in pasGCPList.
136 : * @param pasGCPList an array of GCPs to be used as input.
137 : * @param nReqOrder the requested polynomial order. It should be 1, 2 or 3.
138 : *
139 : * @return the transform argument or NULL if creation fails.
140 : */
141 :
142 : void *GDALCreateGCPTransformer( int nGCPCount, const GDAL_GCP *pasGCPList,
143 : int nReqOrder, int bReversed )
144 :
145 37 : {
146 : GCPTransformInfo *psInfo;
147 : double *padfGeoX, *padfGeoY, *padfRasterX, *padfRasterY;
148 : int *panStatus, iGCP;
149 : int nCRSresult;
150 : struct Control_Points sPoints;
151 :
152 37 : if( nReqOrder == 0 )
153 : {
154 32 : if( nGCPCount >= 10 )
155 0 : nReqOrder = 2; /*for now we avoid 3rd order since it is unstable*/
156 32 : else if( nGCPCount >= 6 )
157 3 : nReqOrder = 2;
158 : else
159 29 : nReqOrder = 1;
160 : }
161 :
162 37 : psInfo = (GCPTransformInfo *) CPLCalloc(sizeof(GCPTransformInfo),1);
163 37 : psInfo->bReversed = bReversed;
164 37 : psInfo->nOrder = nReqOrder;
165 :
166 37 : psInfo->pasGCPList = GDALDuplicateGCPs( nGCPCount, pasGCPList );
167 37 : psInfo->nGCPCount = nGCPCount;
168 :
169 37 : strcpy( psInfo->sTI.szSignature, "GTI" );
170 37 : psInfo->sTI.pszClassName = "GDALGCPTransformer";
171 37 : psInfo->sTI.pfnTransform = GDALGCPTransform;
172 37 : psInfo->sTI.pfnCleanup = GDALDestroyGCPTransformer;
173 37 : psInfo->sTI.pfnSerialize = GDALSerializeGCPTransformer;
174 :
175 : /* -------------------------------------------------------------------- */
176 : /* Allocate and initialize the working points list. */
177 : /* -------------------------------------------------------------------- */
178 37 : padfGeoX = (double *) CPLCalloc(sizeof(double),nGCPCount);
179 37 : padfGeoY = (double *) CPLCalloc(sizeof(double),nGCPCount);
180 37 : padfRasterX = (double *) CPLCalloc(sizeof(double),nGCPCount);
181 37 : padfRasterY = (double *) CPLCalloc(sizeof(double),nGCPCount);
182 37 : panStatus = (int *) CPLCalloc(sizeof(int),nGCPCount);
183 :
184 197 : for( iGCP = 0; iGCP < nGCPCount; iGCP++ )
185 : {
186 160 : panStatus[iGCP] = 1;
187 160 : padfGeoX[iGCP] = pasGCPList[iGCP].dfGCPX;
188 160 : padfGeoY[iGCP] = pasGCPList[iGCP].dfGCPY;
189 160 : padfRasterX[iGCP] = pasGCPList[iGCP].dfGCPPixel;
190 160 : padfRasterY[iGCP] = pasGCPList[iGCP].dfGCPLine;
191 : }
192 :
193 37 : sPoints.count = nGCPCount;
194 37 : sPoints.e1 = padfRasterX;
195 37 : sPoints.n1 = padfRasterY;
196 37 : sPoints.e2 = padfGeoX;
197 37 : sPoints.n2 = padfGeoY;
198 37 : sPoints.status = panStatus;
199 :
200 : /* -------------------------------------------------------------------- */
201 : /* Compute the forward and reverse polynomials. */
202 : /* -------------------------------------------------------------------- */
203 37 : nCRSresult = CRS_compute_georef_equations( &sPoints,
204 : psInfo->adfToGeoX, psInfo->adfToGeoY,
205 : psInfo->adfFromGeoX, psInfo->adfFromGeoY,
206 : nReqOrder );
207 :
208 37 : CPLFree( padfGeoX );
209 37 : CPLFree( padfGeoY );
210 37 : CPLFree( padfRasterX );
211 37 : CPLFree( padfRasterY );
212 37 : CPLFree( panStatus );
213 :
214 37 : if (nCRSresult != 1)
215 : {
216 0 : CPLError( CE_Failure, CPLE_AppDefined, "%s", CRS_error_message[-nCRSresult]);
217 0 : GDALDestroyGCPTransformer( psInfo );
218 0 : return NULL;
219 : }
220 : else
221 : {
222 37 : return psInfo;
223 : }
224 : }
225 :
226 : /************************************************************************/
227 : /* GDALDestroyGCPTransformer() */
228 : /************************************************************************/
229 :
230 : /**
231 : * Destroy GCP transformer.
232 : *
233 : * This function is used to destroy information about a GCP based
234 : * polynomial transformation created with GDALCreateGCPTransformer().
235 : *
236 : * @param pTransformArg the transform arg previously returned by
237 : * GDALCreateGCPTransformer().
238 : */
239 :
240 : void GDALDestroyGCPTransformer( void *pTransformArg )
241 :
242 37 : {
243 37 : GCPTransformInfo *psInfo = (GCPTransformInfo *) pTransformArg;
244 :
245 37 : VALIDATE_POINTER0( pTransformArg, "GDALDestroyGCPTransformer" );
246 :
247 37 : GDALDeinitGCPs( psInfo->nGCPCount, psInfo->pasGCPList );
248 37 : CPLFree( psInfo->pasGCPList );
249 :
250 37 : CPLFree( pTransformArg );
251 : }
252 :
253 : /************************************************************************/
254 : /* GDALGCPTransform() */
255 : /************************************************************************/
256 :
257 : /**
258 : * Transforms point based on GCP derived polynomial model.
259 : *
260 : * This function matches the GDALTransformerFunc signature, and can be
261 : * used to transform one or more points from pixel/line coordinates to
262 : * georeferenced coordinates (SrcToDst) or vice versa (DstToSrc).
263 : *
264 : * @param pTransformArg return value from GDALCreateGCPTransformer().
265 : * @param bDstToSrc TRUE if transformation is from the destination
266 : * (georeferenced) coordinates to pixel/line or FALSE when transforming
267 : * from pixel/line to georeferenced coordinates.
268 : * @param nPointCount the number of values in the x, y and z arrays.
269 : * @param x array containing the X values to be transformed.
270 : * @param y array containing the Y values to be transformed.
271 : * @param z array containing the Z values to be transformed.
272 : * @param panSuccess array in which a flag indicating success (TRUE) or
273 : * failure (FALSE) of the transformation are placed.
274 : *
275 : * @return TRUE.
276 : */
277 :
278 : int GDALGCPTransform( void *pTransformArg, int bDstToSrc,
279 : int nPointCount,
280 : double *x, double *y, double *z,
281 : int *panSuccess )
282 :
283 12566 : {
284 : int i;
285 12566 : GCPTransformInfo *psInfo = (GCPTransformInfo *) pTransformArg;
286 :
287 12566 : if( psInfo->bReversed )
288 0 : bDstToSrc = !bDstToSrc;
289 :
290 294829 : for( i = 0; i < nPointCount; i++ )
291 : {
292 282263 : if( x[i] == HUGE_VAL || y[i] == HUGE_VAL )
293 : {
294 0 : panSuccess[i] = FALSE;
295 0 : continue;
296 : }
297 :
298 282263 : if( bDstToSrc )
299 : {
300 273286 : CRS_georef( x[i], y[i], x + i, y + i,
301 : psInfo->adfFromGeoX, psInfo->adfFromGeoY,
302 : psInfo->nOrder );
303 : }
304 : else
305 : {
306 8977 : CRS_georef( x[i], y[i], x + i, y + i,
307 : psInfo->adfToGeoX, psInfo->adfToGeoY,
308 : psInfo->nOrder );
309 : }
310 282263 : panSuccess[i] = TRUE;
311 : }
312 :
313 12566 : return TRUE;
314 : }
315 :
316 : /************************************************************************/
317 : /* GDALSerializeGCPTransformer() */
318 : /************************************************************************/
319 :
320 : CPLXMLNode *GDALSerializeGCPTransformer( void *pTransformArg )
321 :
322 3 : {
323 : CPLXMLNode *psTree;
324 3 : GCPTransformInfo *psInfo = (GCPTransformInfo *) pTransformArg;
325 :
326 3 : VALIDATE_POINTER1( pTransformArg, "GDALSerializeGCPTransformer", NULL );
327 :
328 3 : psTree = CPLCreateXMLNode( NULL, CXT_Element, "GCPTransformer" );
329 :
330 : /* -------------------------------------------------------------------- */
331 : /* Serialize Order and bReversed. */
332 : /* -------------------------------------------------------------------- */
333 3 : CPLCreateXMLElementAndValue(
334 : psTree, "Order",
335 : CPLSPrintf( "%d", psInfo->nOrder ) );
336 :
337 3 : CPLCreateXMLElementAndValue(
338 : psTree, "Reversed",
339 : CPLSPrintf( "%d", psInfo->bReversed ) );
340 :
341 : /* -------------------------------------------------------------------- */
342 : /* Attach GCP List. */
343 : /* -------------------------------------------------------------------- */
344 3 : if( psInfo->nGCPCount > 0 )
345 : {
346 : int iGCP;
347 : CPLXMLNode *psGCPList = CPLCreateXMLNode( psTree, CXT_Element,
348 3 : "GCPList" );
349 :
350 15 : for( iGCP = 0; iGCP < psInfo->nGCPCount; iGCP++ )
351 : {
352 : CPLXMLNode *psXMLGCP;
353 12 : GDAL_GCP *psGCP = psInfo->pasGCPList + iGCP;
354 :
355 12 : psXMLGCP = CPLCreateXMLNode( psGCPList, CXT_Element, "GCP" );
356 :
357 12 : CPLSetXMLValue( psXMLGCP, "#Id", psGCP->pszId );
358 :
359 12 : if( psGCP->pszInfo != NULL && strlen(psGCP->pszInfo) > 0 )
360 0 : CPLSetXMLValue( psXMLGCP, "Info", psGCP->pszInfo );
361 :
362 12 : CPLSetXMLValue( psXMLGCP, "#Pixel",
363 : CPLSPrintf( "%.4f", psGCP->dfGCPPixel ) );
364 :
365 12 : CPLSetXMLValue( psXMLGCP, "#Line",
366 : CPLSPrintf( "%.4f", psGCP->dfGCPLine ) );
367 :
368 12 : CPLSetXMLValue( psXMLGCP, "#X",
369 : CPLSPrintf( "%.12E", psGCP->dfGCPX ) );
370 :
371 12 : CPLSetXMLValue( psXMLGCP, "#Y",
372 : CPLSPrintf( "%.12E", psGCP->dfGCPY ) );
373 :
374 12 : if( psGCP->dfGCPZ != 0.0 )
375 0 : CPLSetXMLValue( psXMLGCP, "#GCPZ",
376 : CPLSPrintf( "%.12E", psGCP->dfGCPZ ) );
377 : }
378 : }
379 :
380 3 : return psTree;
381 : }
382 :
383 : /************************************************************************/
384 : /* GDALDeserializeReprojectionTransformer() */
385 : /************************************************************************/
386 :
387 : void *GDALDeserializeGCPTransformer( CPLXMLNode *psTree )
388 :
389 4 : {
390 4 : GDAL_GCP *pasGCPList = 0;
391 4 : int nGCPCount = 0;
392 : void *pResult;
393 : int nReqOrder;
394 : int bReversed;
395 :
396 : /* -------------------------------------------------------------------- */
397 : /* Check for GCPs. */
398 : /* -------------------------------------------------------------------- */
399 4 : CPLXMLNode *psGCPList = CPLGetXMLNode( psTree, "GCPList" );
400 :
401 4 : if( psGCPList != NULL )
402 : {
403 4 : int nGCPMax = 0;
404 : CPLXMLNode *psXMLGCP;
405 :
406 : // Count GCPs.
407 24 : for( psXMLGCP = psGCPList->psChild; psXMLGCP != NULL;
408 16 : psXMLGCP = psXMLGCP->psNext )
409 16 : nGCPMax++;
410 :
411 4 : pasGCPList = (GDAL_GCP *) CPLCalloc(sizeof(GDAL_GCP),nGCPMax);
412 :
413 24 : for( psXMLGCP = psGCPList->psChild; psXMLGCP != NULL;
414 16 : psXMLGCP = psXMLGCP->psNext )
415 : {
416 16 : GDAL_GCP *psGCP = pasGCPList + nGCPCount;
417 :
418 16 : if( !EQUAL(psXMLGCP->pszValue,"GCP") ||
419 : psXMLGCP->eType != CXT_Element )
420 : continue;
421 :
422 16 : GDALInitGCPs( 1, psGCP );
423 :
424 16 : CPLFree( psGCP->pszId );
425 16 : psGCP->pszId = CPLStrdup(CPLGetXMLValue(psXMLGCP,"Id",""));
426 :
427 16 : CPLFree( psGCP->pszInfo );
428 16 : psGCP->pszInfo = CPLStrdup(CPLGetXMLValue(psXMLGCP,"Info",""));
429 :
430 16 : psGCP->dfGCPPixel = atof(CPLGetXMLValue(psXMLGCP,"Pixel","0.0"));
431 16 : psGCP->dfGCPLine = atof(CPLGetXMLValue(psXMLGCP,"Line","0.0"));
432 :
433 16 : psGCP->dfGCPX = atof(CPLGetXMLValue(psXMLGCP,"X","0.0"));
434 16 : psGCP->dfGCPY = atof(CPLGetXMLValue(psXMLGCP,"Y","0.0"));
435 16 : psGCP->dfGCPZ = atof(CPLGetXMLValue(psXMLGCP,"Z","0.0"));
436 16 : nGCPCount++;
437 : }
438 : }
439 :
440 : /* -------------------------------------------------------------------- */
441 : /* Get other flags. */
442 : /* -------------------------------------------------------------------- */
443 4 : nReqOrder = atoi(CPLGetXMLValue(psTree,"Order","3"));
444 4 : bReversed = atoi(CPLGetXMLValue(psTree,"Reversed","0"));
445 :
446 : /* -------------------------------------------------------------------- */
447 : /* Generate transformation. */
448 : /* -------------------------------------------------------------------- */
449 4 : pResult = GDALCreateGCPTransformer( nGCPCount, pasGCPList, nReqOrder,
450 : bReversed );
451 :
452 : /* -------------------------------------------------------------------- */
453 : /* Cleanup GCP copy. */
454 : /* -------------------------------------------------------------------- */
455 4 : GDALDeinitGCPs( nGCPCount, pasGCPList );
456 4 : CPLFree( pasGCPList );
457 :
458 4 : return pResult;
459 : }
460 :
461 : /************************************************************************/
462 : /* ==================================================================== */
463 : /* Everything below this point derived from the CRS.C from GRASS. */
464 : /* ==================================================================== */
465 : /************************************************************************/
466 :
467 :
468 : /* STRUCTURE FOR USE INTERNALLY WITH THESE FUNCTIONS. THESE FUNCTIONS EXPECT
469 : SQUARE MATRICES SO ONLY ONE VARIABLE IS GIVEN (N) FOR THE MATRIX SIZE */
470 :
471 : struct MATRIX
472 : {
473 : int n; /* SIZE OF THIS MATRIX (N x N) */
474 : double *v;
475 : };
476 :
477 : /* CALCULATE OFFSET INTO ARRAY BASED ON R/C */
478 :
479 : #define M(row,col) m->v[(((row)-1)*(m->n))+(col)-1]
480 :
481 :
482 : #define MSUCCESS 1 /* SUCCESS */
483 : #define MNPTERR 0 /* NOT ENOUGH POINTS */
484 : #define MUNSOLVABLE -1 /* NOT SOLVABLE */
485 : #define MMEMERR -2 /* NOT ENOUGH MEMORY */
486 : #define MPARMERR -3 /* PARAMETER ERROR */
487 : #define MINTERR -4 /* INTERNAL ERROR */
488 :
489 : /***************************************************************************/
490 : /*
491 : FUNCTION PROTOTYPES FOR STATIC (INTERNAL) FUNCTIONS
492 : */
493 : /***************************************************************************/
494 :
495 : static int calccoef(struct Control_Points *,double *,double *,int);
496 : static int calcls(struct Control_Points *,struct MATRIX *,
497 : double *,double *,double *,double *);
498 : static int exactdet(struct Control_Points *,struct MATRIX *,
499 : double *,double *,double *,double *);
500 : static int solvemat(struct MATRIX *,double *,double *,double *,double *);
501 : static double term(int,double,double);
502 :
503 : /***************************************************************************/
504 : /*
505 : TRANSFORM A SINGLE COORDINATE PAIR.
506 : */
507 : /***************************************************************************/
508 :
509 : static int
510 : CRS_georef (
511 : double e1, /* EASTINGS TO BE TRANSFORMED */
512 : double n1, /* NORTHINGS TO BE TRANSFORMED */
513 : double *e, /* EASTINGS TO BE TRANSFORMED */
514 : double *n, /* NORTHINGS TO BE TRANSFORMED */
515 : double E[], /* EASTING COEFFICIENTS */
516 : double N[], /* NORTHING COEFFICIENTS */
517 : int order /* ORDER OF TRANSFORMATION TO BE PERFORMED, MUST MATCH THE
518 : ORDER USED TO CALCULATE THE COEFFICIENTS */
519 : )
520 282263 : {
521 : double e3, e2n, en2, n3, e2, en, n2;
522 :
523 282263 : switch(order)
524 : {
525 : case 1:
526 :
527 280119 : *e = E[0] + E[1] * e1 + E[2] * n1;
528 280119 : *n = N[0] + N[1] * e1 + N[2] * n1;
529 280119 : break;
530 :
531 : case 2:
532 :
533 2144 : e2 = e1 * e1;
534 2144 : n2 = n1 * n1;
535 2144 : en = e1 * n1;
536 :
537 2144 : *e = E[0] + E[1] * e1 + E[2] * n1 +
538 : E[3] * e2 + E[4] * en + E[5] * n2;
539 2144 : *n = N[0] + N[1] * e1 + N[2] * n1 +
540 : N[3] * e2 + N[4] * en + N[5] * n2;
541 2144 : break;
542 :
543 : case 3:
544 :
545 0 : e2 = e1 * e1;
546 0 : en = e1 * n1;
547 0 : n2 = n1 * n1;
548 0 : e3 = e1 * e2;
549 0 : e2n = e2 * n1;
550 0 : en2 = e1 * n2;
551 0 : n3 = n1 * n2;
552 :
553 0 : *e = E[0] +
554 : E[1] * e1 + E[2] * n1 +
555 : E[3] * e2 + E[4] * en + E[5] * n2 +
556 : E[6] * e3 + E[7] * e2n + E[8] * en2 + E[9] * n3;
557 0 : *n = N[0] +
558 : N[1] * e1 + N[2] * n1 +
559 : N[3] * e2 + N[4] * en + N[5] * n2 +
560 : N[6] * e3 + N[7] * e2n + N[8] * en2 + N[9] * n3;
561 0 : break;
562 :
563 : default:
564 :
565 0 : return(MPARMERR);
566 : }
567 :
568 282263 : return(MSUCCESS);
569 : }
570 :
571 : /***************************************************************************/
572 : /*
573 : COMPUTE THE GEOREFFERENCING COEFFICIENTS BASED ON A SET OF CONTROL POINTS
574 : */
575 : /***************************************************************************/
576 :
577 : static int
578 : CRS_compute_georef_equations (struct Control_Points *cp,
579 : double E12[], double N12[],
580 : double E21[], double N21[],
581 : int order)
582 37 : {
583 : double *tempptr;
584 : int status;
585 :
586 37 : if(order < 1 || order > MAXORDER)
587 0 : return(MPARMERR);
588 :
589 : /* CALCULATE THE FORWARD TRANSFORMATION COEFFICIENTS */
590 :
591 37 : status = calccoef(cp,E12,N12,order);
592 37 : if(status != MSUCCESS)
593 0 : return(status);
594 :
595 : /* SWITCH THE 1 AND 2 EASTING AND NORTHING ARRAYS */
596 :
597 37 : tempptr = cp->e1;
598 37 : cp->e1 = cp->e2;
599 37 : cp->e2 = tempptr;
600 37 : tempptr = cp->n1;
601 37 : cp->n1 = cp->n2;
602 37 : cp->n2 = tempptr;
603 :
604 : /* CALCULATE THE BACKWARD TRANSFORMATION COEFFICIENTS */
605 :
606 37 : status = calccoef(cp,E21,N21,order);
607 :
608 : /* SWITCH THE 1 AND 2 EASTING AND NORTHING ARRAYS BACK */
609 :
610 37 : tempptr = cp->e1;
611 37 : cp->e1 = cp->e2;
612 37 : cp->e2 = tempptr;
613 37 : tempptr = cp->n1;
614 37 : cp->n1 = cp->n2;
615 37 : cp->n2 = tempptr;
616 :
617 37 : return(status);
618 : }
619 :
620 : /***************************************************************************/
621 : /*
622 : COMPUTE THE GEOREFFERENCING COEFFICIENTS BASED ON A SET OF CONTROL POINTS
623 : */
624 : /***************************************************************************/
625 :
626 : static int
627 : calccoef (struct Control_Points *cp, double E[], double N[], int order)
628 74 : {
629 : struct MATRIX m;
630 : double *a;
631 : double *b;
632 : int numactive; /* NUMBER OF ACTIVE CONTROL POINTS */
633 : int status, i;
634 :
635 : /* CALCULATE THE NUMBER OF VALID CONTROL POINTS */
636 :
637 394 : for(i = numactive = 0 ; i < cp->count ; i++)
638 : {
639 320 : if(cp->status[i] > 0)
640 320 : numactive++;
641 : }
642 :
643 : /* CALCULATE THE MINIMUM NUMBER OF CONTROL POINTS NEEDED TO DETERMINE
644 : A TRANSFORMATION OF THIS ORDER */
645 :
646 74 : m.n = ((order + 1) * (order + 2)) / 2;
647 :
648 74 : if(numactive < m.n)
649 0 : return(MNPTERR);
650 :
651 : /* INITIALIZE MATRIX */
652 :
653 74 : m.v = (double *)CPLCalloc(m.n*m.n,sizeof(double));
654 74 : if(m.v == NULL)
655 : {
656 0 : return(MMEMERR);
657 : }
658 74 : a = (double *)CPLCalloc(m.n,sizeof(double));
659 74 : if(a == NULL)
660 : {
661 0 : CPLFree((char *)m.v);
662 0 : return(MMEMERR);
663 : }
664 74 : b = (double *)CPLCalloc(m.n,sizeof(double));
665 74 : if(b == NULL)
666 : {
667 0 : CPLFree((char *)m.v);
668 0 : CPLFree((char *)a);
669 0 : return(MMEMERR);
670 : }
671 :
672 74 : if(numactive == m.n)
673 0 : status = exactdet(cp,&m,a,b,E,N);
674 : else
675 74 : status = calcls(cp,&m,a,b,E,N);
676 :
677 74 : CPLFree((char *)m.v);
678 74 : CPLFree((char *)a);
679 74 : CPLFree((char *)b);
680 :
681 74 : return(status);
682 : }
683 :
684 : /***************************************************************************/
685 : /*
686 : CALCULATE THE TRANSFORMATION COEFFICIENTS WITH EXACTLY THE MINIMUM
687 : NUMBER OF CONTROL POINTS REQUIRED FOR THIS TRANSFORMATION.
688 : */
689 : /***************************************************************************/
690 :
691 : static int exactdet (
692 : struct Control_Points *cp,
693 : struct MATRIX *m,
694 : double a[],
695 : double b[],
696 : double E[], /* EASTING COEFFICIENTS */
697 : double N[] /* NORTHING COEFFICIENTS */
698 : )
699 0 : {
700 : int pntnow, currow, j;
701 :
702 0 : currow = 1;
703 0 : for(pntnow = 0 ; pntnow < cp->count ; pntnow++)
704 : {
705 0 : if(cp->status[pntnow] > 0)
706 : {
707 : /* POPULATE MATRIX M */
708 :
709 0 : for(j = 1 ; j <= m->n ; j++)
710 : {
711 0 : M(currow,j) = term(j,cp->e1[pntnow],cp->n1[pntnow]);
712 : }
713 :
714 : /* POPULATE MATRIX A AND B */
715 :
716 0 : a[currow-1] = cp->e2[pntnow];
717 0 : b[currow-1] = cp->n2[pntnow];
718 :
719 0 : currow++;
720 : }
721 : }
722 :
723 0 : if(currow - 1 != m->n)
724 0 : return(MINTERR);
725 :
726 0 : return(solvemat(m,a,b,E,N));
727 : }
728 :
729 : /***************************************************************************/
730 : /*
731 : CALCULATE THE TRANSFORMATION COEFFICIENTS WITH MORE THAN THE MINIMUM
732 : NUMBER OF CONTROL POINTS REQUIRED FOR THIS TRANSFORMATION. THIS
733 : ROUTINE USES THE LEAST SQUARES METHOD TO COMPUTE THE COEFFICIENTS.
734 : */
735 : /***************************************************************************/
736 :
737 : static int calcls (
738 : struct Control_Points *cp,
739 : struct MATRIX *m,
740 : double a[],
741 : double b[],
742 : double E[], /* EASTING COEFFICIENTS */
743 : double N[] /* NORTHING COEFFICIENTS */
744 : )
745 74 : {
746 74 : int i, j, n, numactive = 0;
747 :
748 : /* INITIALIZE THE UPPER HALF OF THE MATRIX AND THE TWO COLUMN VECTORS */
749 :
750 314 : for(i = 1 ; i <= m->n ; i++)
751 : {
752 774 : for(j = i ; j <= m->n ; j++)
753 534 : M(i,j) = 0.0;
754 240 : a[i-1] = b[i-1] = 0.0;
755 : }
756 :
757 : /* SUM THE UPPER HALF OF THE MATRIX AND THE COLUMN VECTORS ACCORDING TO
758 : THE LEAST SQUARES METHOD OF SOLVING OVER DETERMINED SYSTEMS */
759 :
760 394 : for(n = 0 ; n < cp->count ; n++)
761 : {
762 320 : if(cp->status[n] > 0)
763 : {
764 320 : numactive++;
765 1424 : for(i = 1 ; i <= m->n ; i++)
766 : {
767 3744 : for(j = i ; j <= m->n ; j++)
768 2640 : M(i,j) += term(i,cp->e1[n],cp->n1[n]) * term(j,cp->e1[n],cp->n1[n]);
769 :
770 1104 : a[i-1] += cp->e2[n] * term(i,cp->e1[n],cp->n1[n]);
771 1104 : b[i-1] += cp->n2[n] * term(i,cp->e1[n],cp->n1[n]);
772 : }
773 : }
774 : }
775 :
776 74 : if(numactive <= m->n)
777 0 : return(MINTERR);
778 :
779 : /* TRANSPOSE VALUES IN UPPER HALF OF M TO OTHER HALF */
780 :
781 240 : for(i = 2 ; i <= m->n ; i++)
782 : {
783 460 : for(j = 1 ; j < i ; j++)
784 294 : M(i,j) = M(j,i);
785 : }
786 :
787 74 : return(solvemat(m,a,b,E,N));
788 : }
789 :
790 : /***************************************************************************/
791 : /*
792 : CALCULATE THE X/Y TERM BASED ON THE TERM NUMBER
793 :
794 : ORDER\TERM 1 2 3 4 5 6 7 8 9 10
795 : 1 e0n0 e1n0 e0n1
796 : 2 e0n0 e1n0 e0n1 e2n0 e1n1 e0n2
797 : 3 e0n0 e1n0 e0n1 e2n0 e1n1 e0n2 e3n0 e2n1 e1n2 e0n3
798 : */
799 : /***************************************************************************/
800 :
801 : static double term (int term, double e, double n)
802 7488 : {
803 7488 : switch(term)
804 : {
805 2064 : case 1: return((double)1.0);
806 2064 : case 2: return((double)e);
807 2064 : case 3: return((double)n);
808 432 : case 4: return((double)(e*e));
809 432 : case 5: return((double)(e*n));
810 432 : case 6: return((double)(n*n));
811 0 : case 7: return((double)(e*e*e));
812 0 : case 8: return((double)(e*e*n));
813 0 : case 9: return((double)(e*n*n));
814 0 : case 10: return((double)(n*n*n));
815 : }
816 0 : return((double)0.0);
817 : }
818 :
819 : /***************************************************************************/
820 : /*
821 : SOLVE FOR THE 'E' AND 'N' COEFFICIENTS BY USING A SOMEWHAT MODIFIED
822 : GAUSSIAN ELIMINATION METHOD.
823 :
824 : | M11 M12 ... M1n | | E0 | | a0 |
825 : | M21 M22 ... M2n | | E1 | = | a1 |
826 : | . . . . | | . | | . |
827 : | Mn1 Mn2 ... Mnn | | En-1 | | an-1 |
828 :
829 : and
830 :
831 : | M11 M12 ... M1n | | N0 | | b0 |
832 : | M21 M22 ... M2n | | N1 | = | b1 |
833 : | . . . . | | . | | . |
834 : | Mn1 Mn2 ... Mnn | | Nn-1 | | bn-1 |
835 : */
836 : /***************************************************************************/
837 :
838 : static int solvemat (struct MATRIX *m,
839 : double a[], double b[], double E[], double N[])
840 74 : {
841 : int i, j, i2, j2, imark;
842 : double factor, temp;
843 : double pivot; /* ACTUAL VALUE OF THE LARGEST PIVOT CANDIDATE */
844 :
845 314 : for(i = 1 ; i <= m->n ; i++)
846 : {
847 240 : j = i;
848 :
849 : /* find row with largest magnitude value for pivot value */
850 :
851 240 : pivot = M(i,j);
852 240 : imark = i;
853 534 : for(i2 = i + 1 ; i2 <= m->n ; i2++)
854 : {
855 294 : temp = fabs(M(i2,j));
856 294 : if(temp > fabs(pivot))
857 : {
858 184 : pivot = M(i2,j);
859 184 : imark = i2;
860 : }
861 : }
862 :
863 : /* if the pivot is very small then the points are nearly co-linear */
864 : /* co-linear points result in an undefined matrix, and nearly */
865 : /* co-linear points results in a solution with rounding error */
866 :
867 240 : if(pivot == 0.0)
868 0 : return(MUNSOLVABLE);
869 :
870 : /* if row with highest pivot is not the current row, switch them */
871 :
872 240 : if(imark != i)
873 : {
874 588 : for(j2 = 1 ; j2 <= m->n ; j2++)
875 : {
876 459 : temp = M(imark,j2);
877 459 : M(imark,j2) = M(i,j2);
878 459 : M(i,j2) = temp;
879 : }
880 :
881 129 : temp = a[imark-1];
882 129 : a[imark-1] = a[i-1];
883 129 : a[i-1] = temp;
884 :
885 129 : temp = b[imark-1];
886 129 : b[imark-1] = b[i-1];
887 129 : b[i-1] = temp;
888 : }
889 :
890 : /* compute zeros above and below the pivot, and compute
891 : values for the rest of the row as well */
892 :
893 1068 : for(i2 = 1 ; i2 <= m->n ; i2++)
894 : {
895 828 : if(i2 != i)
896 : {
897 588 : factor = M(i2,j) / pivot;
898 2034 : for(j2 = j ; j2 <= m->n ; j2++)
899 1446 : M(i2,j2) -= factor * M(i,j2);
900 588 : a[i2-1] -= factor * a[i-1];
901 588 : b[i2-1] -= factor * b[i-1];
902 : }
903 : }
904 : }
905 :
906 : /* SINCE ALL OTHER VALUES IN THE MATRIX ARE ZERO NOW, CALCULATE THE
907 : COEFFICIENTS BY DIVIDING THE COLUMN VECTORS BY THE DIAGONAL VALUES. */
908 :
909 314 : for(i = 1 ; i <= m->n ; i++)
910 : {
911 240 : E[i-1] = a[i-1] / M(i,i);
912 240 : N[i-1] = b[i-1] / M(i,i);
913 : }
914 :
915 74 : return(MSUCCESS);
916 : }
|
