LCOV - code coverage report
Current view: directory - alg - gdal_crs.c (source / functions) Found Hit Coverage
Test: gdal_filtered.info Lines: 373 324 86.9 %
Date: 2012-12-26 Functions: 16 16 100.0 %

       1                 : /******************************************************************************
       2                 :  * $Id: gdal_crs.c 22670 2011-07-08 17:15:19Z 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                 :     Last Update:  6/29/2011 C. F. Stallmann & R. van den Dool (South African National Space Agency)
      30                 :       GCP refinement added
      31                 :       
      32                 : 
      33                 :     Copyright (c) 1992, Michigan State University
      34                 :    
      35                 :     Permission is hereby granted, free of charge, to any person obtaining a
      36                 :     copy of this software and associated documentation files (the "Software"),
      37                 :     to deal in the Software without restriction, including without limitation
      38                 :     the rights to use, copy, modify, merge, publish, distribute, sublicense,
      39                 :     and/or sell copies of the Software, and to permit persons to whom the
      40                 :     Software is furnished to do so, subject to the following conditions:
      41                 :     
      42                 :     The above copyright notice and this permission notice shall be included
      43                 :     in all copies or substantial portions of the Software.
      44                 :     
      45                 :     THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
      46                 :     IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
      47                 :     FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.  IN NO EVENT SHALL
      48                 :     THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
      49                 :     LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
      50                 :     FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER 
      51                 :     DEALINGS IN THE SOFTWARE.
      52                 : 
      53                 :  ****************************************************************************/
      54                 : 
      55                 : #include "gdal_alg.h"
      56                 : #include "cpl_conv.h"
      57                 : #include "cpl_minixml.h"
      58                 : #include "cpl_string.h"
      59                 : 
      60                 : CPL_CVSID("$Id: gdal_crs.c 22670 2011-07-08 17:15:19Z rouault $");
      61                 : 
      62                 : #define MAXORDER 3
      63                 : 
      64                 : struct Control_Points
      65                 : {
      66                 :     int  count;
      67                 :     double *e1;
      68                 :     double *n1;
      69                 :     double *e2;
      70                 :     double *n2;
      71                 :     int *status;
      72                 : };
      73                 : 
      74                 : typedef struct
      75                 : {
      76                 :     GDALTransformerInfo sTI;
      77                 : 
      78                 :     double adfToGeoX[20];
      79                 :     double adfToGeoY[20];
      80                 :     
      81                 :     double adfFromGeoX[20];
      82                 :     double adfFromGeoY[20];
      83                 : 
      84                 :     int    nOrder;
      85                 :     int    bReversed;
      86                 : 
      87                 :     int       nGCPCount;
      88                 :     GDAL_GCP *pasGCPList;
      89                 :     int    bRefine;
      90                 :     int    nMinimumGcps;
      91                 :     double dfTolerance;
      92                 :     
      93                 : } GCPTransformInfo;
      94                 : 
      95                 : CPL_C_START
      96                 : CPLXMLNode *GDALSerializeGCPTransformer( void *pTransformArg );
      97                 : void *GDALDeserializeGCPTransformer( CPLXMLNode *psTree );
      98                 : CPL_C_END
      99                 : 
     100                 : /* crs.c */
     101                 : static int CRS_georef(double, double, double *, double *, 
     102                 :                               double [], double [], int);
     103                 : static int CRS_compute_georef_equations(struct Control_Points *,
     104                 :     double [], double [], double [], double [], int);
     105                 : static int remove_outliers(GCPTransformInfo *);
     106                 : 
     107                 : static char *CRS_error_message[] = {
     108                 :     "Failed to compute GCP transform: Not enough points available",
     109                 :     "Failed to compute GCP transform: Transform is not solvable",
     110                 :     "Failed to compute GCP transform: Not enough memory",
     111                 :     "Failed to compute GCP transform: Parameter error",
     112                 :     "Failed to compute GCP transform: Internal error"
     113                 : };
     114                 : 
     115                 : 
     116                 : /************************************************************************/
     117                 : /*                      GDALCreateGCPTransformer()                      */
     118                 : /************************************************************************/
     119                 : 
     120                 : /**
     121                 :  * Create GCP based polynomial transformer.
     122                 :  *
     123                 :  * Computes least squares fit polynomials from a provided set of GCPs,
     124                 :  * and stores the coefficients for later transformation of points between
     125                 :  * pixel/line and georeferenced coordinates. 
     126                 :  *
     127                 :  * The return value should be used as a TransformArg in combination with
     128                 :  * the transformation function GDALGCPTransform which fits the 
     129                 :  * GDALTransformerFunc signature.  The returned transform argument should
     130                 :  * be deallocated with GDALDestroyGCPTransformer when no longer needed.
     131                 :  *
     132                 :  * This function may fail (returning NULL) if the provided set of GCPs
     133                 :  * are inadequate for the requested order, the determinate is zero or they
     134                 :  * are otherwise "ill conditioned".  
     135                 :  *
     136                 :  * Note that 2nd order requires at least 6 GCPs, and 3rd order requires at
     137                 :  * least 10 gcps.  If nReqOrder is 0 the highest order possible with the
     138                 :  * provided gcp count will be used.
     139                 :  *
     140                 :  * @param nGCPCount the number of GCPs in pasGCPList.
     141                 :  * @param pasGCPList an array of GCPs to be used as input.
     142                 :  * @param nReqOrder the requested polynomial order.  It should be 1, 2 or 3.
     143                 :  * 
     144                 :  * @return the transform argument or NULL if creation fails. 
     145                 :  */
     146                 : 
     147              27 : void *GDALCreateGCPTransformerEx( int nGCPCount, const GDAL_GCP *pasGCPList, 
     148                 :                                 int nReqOrder, int bReversed, int bRefine, double dfTolerance, int nMinimumGcps)
     149                 : 
     150                 : {
     151                 :     GCPTransformInfo *psInfo;
     152                 :     double *padfGeoX, *padfGeoY, *padfRasterX, *padfRasterY;
     153                 :     int    *panStatus, iGCP;
     154                 :     int    nCRSresult;
     155                 :     struct Control_Points sPoints;
     156                 : 
     157              27 :     if( nReqOrder == 0 )
     158                 :     {
     159              20 :         if( nGCPCount >= 10 )
     160               0 :             nReqOrder = 2; /*for now we avoid 3rd order since it is unstable*/
     161              20 :         else if( nGCPCount >= 6 )
     162               1 :             nReqOrder = 2;
     163                 :         else
     164              19 :             nReqOrder = 1;
     165                 :     }
     166                 :     
     167              27 :     psInfo = (GCPTransformInfo *) CPLCalloc(sizeof(GCPTransformInfo),1);
     168              27 :     psInfo->bReversed = bReversed;
     169              27 :     psInfo->nOrder = nReqOrder;
     170              27 :     psInfo->bRefine = bRefine;
     171              27 :     psInfo->dfTolerance = dfTolerance;
     172              27 :     psInfo->nMinimumGcps = nMinimumGcps;
     173                 : 
     174              27 :     psInfo->pasGCPList = GDALDuplicateGCPs( nGCPCount, pasGCPList );
     175              27 :     psInfo->nGCPCount = nGCPCount;
     176                 : 
     177              27 :     strcpy( psInfo->sTI.szSignature, "GTI" );
     178              27 :     psInfo->sTI.pszClassName = "GDALGCPTransformer";
     179              27 :     psInfo->sTI.pfnTransform = GDALGCPTransform;
     180              27 :     psInfo->sTI.pfnCleanup = GDALDestroyGCPTransformer;
     181              27 :     psInfo->sTI.pfnSerialize = GDALSerializeGCPTransformer;
     182                 :     
     183                 : /* -------------------------------------------------------------------- */
     184                 : /*      Compute the forward and reverse polynomials.                    */
     185                 : /* -------------------------------------------------------------------- */
     186                 : 
     187              27 :     if(bRefine)
     188                 :     {
     189               1 :         nCRSresult = remove_outliers(psInfo);
     190                 :     }
     191                 :     else
     192                 :     {
     193                 :         /* -------------------------------------------------------------------- */
     194                 :         /*      Allocate and initialize the working points list.                */
     195                 :         /* -------------------------------------------------------------------- */
     196              26 :         padfGeoX = (double *) CPLCalloc(sizeof(double),nGCPCount);
     197              26 :         padfGeoY = (double *) CPLCalloc(sizeof(double),nGCPCount);
     198              26 :         padfRasterX = (double *) CPLCalloc(sizeof(double),nGCPCount);
     199              26 :         padfRasterY = (double *) CPLCalloc(sizeof(double),nGCPCount);
     200              26 :         panStatus = (int *) CPLCalloc(sizeof(int),nGCPCount);
     201                 :     
     202             130 :         for( iGCP = 0; iGCP < nGCPCount; iGCP++ )
     203                 :         {
     204             104 :             panStatus[iGCP] = 1;
     205             104 :             padfGeoX[iGCP] = pasGCPList[iGCP].dfGCPX;
     206             104 :             padfGeoY[iGCP] = pasGCPList[iGCP].dfGCPY;
     207             104 :             padfRasterX[iGCP] = pasGCPList[iGCP].dfGCPPixel;
     208             104 :             padfRasterY[iGCP] = pasGCPList[iGCP].dfGCPLine;
     209                 :         }
     210                 : 
     211              26 :         sPoints.count = nGCPCount;
     212              26 :         sPoints.e1 = padfRasterX;
     213              26 :         sPoints.n1 = padfRasterY;
     214              26 :         sPoints.e2 = padfGeoX;
     215              26 :         sPoints.n2 = padfGeoY;
     216              26 :         sPoints.status = panStatus;
     217              26 :         nCRSresult = CRS_compute_georef_equations( &sPoints,
     218                 :                                                 psInfo->adfToGeoX, psInfo->adfToGeoY,
     219                 :                                                 psInfo->adfFromGeoX, psInfo->adfFromGeoY,
     220                 :                                                 nReqOrder );
     221              26 :         CPLFree( padfGeoX );
     222              26 :         CPLFree( padfGeoY );
     223              26 :         CPLFree( padfRasterX );
     224              26 :         CPLFree( padfRasterY );
     225              26 :         CPLFree( panStatus );
     226                 :     }
     227                 : 
     228              27 :     if (nCRSresult != 1)
     229                 :     {
     230               0 :         CPLError( CE_Failure, CPLE_AppDefined, "%s", CRS_error_message[-nCRSresult]);
     231               0 :         GDALDestroyGCPTransformer( psInfo );
     232               0 :         return NULL;
     233                 :     }
     234                 :     else
     235                 :     {
     236              27 :         return psInfo;
     237                 :     }
     238                 : }
     239                 : 
     240              26 : void *GDALCreateGCPTransformer( int nGCPCount, const GDAL_GCP *pasGCPList, 
     241                 :                                 int nReqOrder, int bReversed )
     242                 : 
     243                 : {
     244              26 :     return GDALCreateGCPTransformerEx(nGCPCount, pasGCPList, nReqOrder, bReversed, FALSE, -1, -1);
     245                 : }
     246                 : 
     247               1 : void *GDALCreateGCPRefineTransformer( int nGCPCount, const GDAL_GCP *pasGCPList, 
     248                 :                                 int nReqOrder, int bReversed, double dfTolerance, int nMinimumGcps)
     249                 : 
     250                 : {
     251                 :     //If no minimumGcp parameter was passed, we  use the default value according to the model
     252               1 :     if(nMinimumGcps == -1)
     253                 :     {
     254               0 :         nMinimumGcps = ((nReqOrder+1) * (nReqOrder+2)) / 2 + 1;
     255                 :     }
     256               1 :     return GDALCreateGCPTransformerEx(nGCPCount, pasGCPList, nReqOrder, bReversed, TRUE, dfTolerance, nMinimumGcps);
     257                 : }
     258                 : 
     259                 : 
     260                 : /************************************************************************/
     261                 : /*                     GDALDestroyGCPTransformer()                      */
     262                 : /************************************************************************/
     263                 : 
     264                 : /**
     265                 :  * Destroy GCP transformer.
     266                 :  *
     267                 :  * This function is used to destroy information about a GCP based
     268                 :  * polynomial transformation created with GDALCreateGCPTransformer(). 
     269                 :  *
     270                 :  * @param pTransformArg the transform arg previously returned by 
     271                 :  * GDALCreateGCPTransformer(). 
     272                 :  */
     273                 : 
     274              27 : void GDALDestroyGCPTransformer( void *pTransformArg )
     275                 : 
     276                 : {
     277              27 :     GCPTransformInfo *psInfo = (GCPTransformInfo *) pTransformArg;
     278                 : 
     279              27 :     VALIDATE_POINTER0( pTransformArg, "GDALDestroyGCPTransformer" );
     280                 : 
     281              27 :     GDALDeinitGCPs( psInfo->nGCPCount, psInfo->pasGCPList );
     282              27 :     CPLFree( psInfo->pasGCPList );
     283                 : 
     284              27 :     CPLFree( pTransformArg );
     285                 : }
     286                 : 
     287                 : /************************************************************************/
     288                 : /*                          GDALGCPTransform()                          */
     289                 : /************************************************************************/
     290                 : 
     291                 : /**
     292                 :  * Transforms point based on GCP derived polynomial model.
     293                 :  *
     294                 :  * This function matches the GDALTransformerFunc signature, and can be
     295                 :  * used to transform one or more points from pixel/line coordinates to
     296                 :  * georeferenced coordinates (SrcToDst) or vice versa (DstToSrc). 
     297                 :  *
     298                 :  * @param pTransformArg return value from GDALCreateGCPTransformer(). 
     299                 :  * @param bDstToSrc TRUE if transformation is from the destination 
     300                 :  * (georeferenced) coordinates to pixel/line or FALSE when transforming
     301                 :  * from pixel/line to georeferenced coordinates.
     302                 :  * @param nPointCount the number of values in the x, y and z arrays.
     303                 :  * @param x array containing the X values to be transformed.
     304                 :  * @param y array containing the Y values to be transformed.
     305                 :  * @param z array containing the Z values to be transformed.
     306                 :  * @param panSuccess array in which a flag indicating success (TRUE) or
     307                 :  * failure (FALSE) of the transformation are placed.
     308                 :  *
     309                 :  * @return TRUE.
     310                 :  */
     311                 : 
     312           12223 : int GDALGCPTransform( void *pTransformArg, int bDstToSrc, 
     313                 :                       int nPointCount, 
     314                 :                       double *x, double *y, double *z, 
     315                 :                       int *panSuccess )
     316                 : 
     317                 : {
     318                 :     int    i;
     319           12223 :     GCPTransformInfo *psInfo = (GCPTransformInfo *) pTransformArg;
     320                 : 
     321           12223 :     if( psInfo->bReversed )
     322               0 :         bDstToSrc = !bDstToSrc;
     323                 :     
     324          103106 :     for( i = 0; i < nPointCount; i++ )
     325                 :     {
     326           90883 :         if( x[i] == HUGE_VAL || y[i] == HUGE_VAL )
     327                 :         {
     328               0 :             panSuccess[i] = FALSE;
     329               0 :             continue;
     330                 :         }
     331                 : 
     332           90883 :         if( bDstToSrc )
     333                 :         {
     334           81898 :             CRS_georef( x[i], y[i], x + i, y + i, 
     335                 :                         psInfo->adfFromGeoX, psInfo->adfFromGeoY, 
     336                 :                         psInfo->nOrder );
     337                 :         }
     338                 :         else
     339                 :         {
     340            8985 :             CRS_georef( x[i], y[i], x + i, y + i, 
     341                 :                         psInfo->adfToGeoX, psInfo->adfToGeoY, 
     342                 :                         psInfo->nOrder );
     343                 :         }
     344           90883 :         panSuccess[i] = TRUE;
     345                 :     }
     346                 : 
     347           12223 :     return TRUE;
     348                 : }
     349                 : 
     350                 : /************************************************************************/
     351                 : /*                    GDALSerializeGCPTransformer()                     */
     352                 : /************************************************************************/
     353                 : 
     354               3 : CPLXMLNode *GDALSerializeGCPTransformer( void *pTransformArg )
     355                 : 
     356                 : {
     357                 :     CPLXMLNode *psTree;
     358               3 :     GCPTransformInfo *psInfo = (GCPTransformInfo *) pTransformArg;
     359                 : 
     360               3 :     VALIDATE_POINTER1( pTransformArg, "GDALSerializeGCPTransformer", NULL );
     361                 : 
     362               3 :     psTree = CPLCreateXMLNode( NULL, CXT_Element, "GCPTransformer" );
     363                 : 
     364                 : /* -------------------------------------------------------------------- */
     365                 : /*      Serialize Order and bReversed.                                  */
     366                 : /* -------------------------------------------------------------------- */
     367               3 :     CPLCreateXMLElementAndValue( 
     368                 :         psTree, "Order", 
     369                 :         CPLSPrintf( "%d", psInfo->nOrder ) );
     370                 :                          
     371               3 :     CPLCreateXMLElementAndValue( 
     372                 :         psTree, "Reversed", 
     373                 :         CPLSPrintf( "%d", psInfo->bReversed ) );
     374                 : 
     375               3 :     if( psInfo->bRefine )
     376                 :     {
     377               0 :         CPLCreateXMLElementAndValue(
     378                 :             psTree, "Refine",
     379                 :             CPLSPrintf( "%d", psInfo->bRefine ) );
     380                 : 
     381               0 :         CPLCreateXMLElementAndValue(
     382                 :             psTree, "MinimumGcps",
     383                 :             CPLSPrintf( "%d", psInfo->nMinimumGcps ) );
     384                 : 
     385               0 :         CPLCreateXMLElementAndValue(
     386                 :             psTree, "Tolerance",
     387                 :             CPLSPrintf( "%f", psInfo->dfTolerance ) );
     388                 :     }
     389                 :                                  
     390                 : /* -------------------------------------------------------------------- */
     391                 : /*  Attach GCP List.            */
     392                 : /* -------------------------------------------------------------------- */
     393               3 :     if( psInfo->nGCPCount > 0 )
     394                 :     {
     395                 :         int iGCP;
     396               3 :         CPLXMLNode *psGCPList = CPLCreateXMLNode( psTree, CXT_Element, 
     397               3 :                                                   "GCPList" );
     398                 : 
     399               3 :   if(psInfo->bRefine)
     400                 :   {
     401               0 :     remove_outliers(psInfo);
     402                 :   }
     403                 :   
     404              15 :         for( iGCP = 0; iGCP < psInfo->nGCPCount; iGCP++ )
     405                 :         {
     406                 :             CPLXMLNode *psXMLGCP;
     407              12 :             GDAL_GCP *psGCP = psInfo->pasGCPList + iGCP;
     408                 : 
     409              12 :             psXMLGCP = CPLCreateXMLNode( psGCPList, CXT_Element, "GCP" );
     410                 : 
     411              12 :             CPLSetXMLValue( psXMLGCP, "#Id", psGCP->pszId );
     412                 : 
     413              12 :             if( psGCP->pszInfo != NULL && strlen(psGCP->pszInfo) > 0 )
     414               0 :                 CPLSetXMLValue( psXMLGCP, "Info", psGCP->pszInfo );
     415                 : 
     416              12 :             CPLSetXMLValue( psXMLGCP, "#Pixel", 
     417                 :                             CPLSPrintf( "%.4f", psGCP->dfGCPPixel ) );
     418                 : 
     419              12 :             CPLSetXMLValue( psXMLGCP, "#Line", 
     420                 :                             CPLSPrintf( "%.4f", psGCP->dfGCPLine ) );
     421                 : 
     422              12 :             CPLSetXMLValue( psXMLGCP, "#X", 
     423                 :                             CPLSPrintf( "%.12E", psGCP->dfGCPX ) );
     424                 : 
     425              12 :             CPLSetXMLValue( psXMLGCP, "#Y", 
     426                 :                             CPLSPrintf( "%.12E", psGCP->dfGCPY ) );
     427                 : 
     428              12 :             if( psGCP->dfGCPZ != 0.0 )
     429               0 :                 CPLSetXMLValue( psXMLGCP, "#GCPZ", 
     430                 :                                 CPLSPrintf( "%.12E", psGCP->dfGCPZ ) );
     431                 :         }
     432                 :     }
     433                 : 
     434               3 :     return psTree;
     435                 : }
     436                 : 
     437                 : /************************************************************************/
     438                 : /*               GDALDeserializeReprojectionTransformer()               */
     439                 : /************************************************************************/
     440                 : 
     441               5 : void *GDALDeserializeGCPTransformer( CPLXMLNode *psTree )
     442                 : 
     443                 : {
     444               5 :     GDAL_GCP *pasGCPList = 0;
     445               5 :     int nGCPCount = 0;
     446                 :     void *pResult;
     447                 :     int nReqOrder;
     448                 :     int bReversed;
     449                 :     int bRefine;
     450                 :     int nMinimumGcps;
     451                 :     double dfTolerance;
     452                 : 
     453                 :     /* -------------------------------------------------------------------- */
     454                 :     /*      Check for GCPs.                                                 */
     455                 :     /* -------------------------------------------------------------------- */
     456               5 :     CPLXMLNode *psGCPList = CPLGetXMLNode( psTree, "GCPList" );
     457                 : 
     458               5 :     if( psGCPList != NULL )
     459                 :     {
     460               5 :         int  nGCPMax = 0;
     461                 :         CPLXMLNode *psXMLGCP;
     462                 :          
     463                 :         // Count GCPs.
     464              36 :         for( psXMLGCP = psGCPList->psChild; psXMLGCP != NULL; 
     465              26 :              psXMLGCP = psXMLGCP->psNext )
     466              26 :             nGCPMax++;
     467                 :          
     468               5 :         pasGCPList = (GDAL_GCP *) CPLCalloc(sizeof(GDAL_GCP),nGCPMax);
     469                 : 
     470              36 :         for( psXMLGCP = psGCPList->psChild; psXMLGCP != NULL; 
     471              26 :              psXMLGCP = psXMLGCP->psNext )
     472                 :         {
     473              26 :             GDAL_GCP *psGCP = pasGCPList + nGCPCount;
     474                 : 
     475              52 :             if( !EQUAL(psXMLGCP->pszValue,"GCP") || 
     476              26 :                 psXMLGCP->eType != CXT_Element )
     477               0 :                 continue;
     478                 :              
     479              26 :             GDALInitGCPs( 1, psGCP );
     480                 :              
     481              26 :             CPLFree( psGCP->pszId );
     482              26 :             psGCP->pszId = CPLStrdup(CPLGetXMLValue(psXMLGCP,"Id",""));
     483                 :              
     484              26 :             CPLFree( psGCP->pszInfo );
     485              26 :             psGCP->pszInfo = CPLStrdup(CPLGetXMLValue(psXMLGCP,"Info",""));
     486                 :              
     487              26 :             psGCP->dfGCPPixel = atof(CPLGetXMLValue(psXMLGCP,"Pixel","0.0"));
     488              26 :             psGCP->dfGCPLine = atof(CPLGetXMLValue(psXMLGCP,"Line","0.0"));
     489                 :              
     490              26 :             psGCP->dfGCPX = atof(CPLGetXMLValue(psXMLGCP,"X","0.0"));
     491              26 :             psGCP->dfGCPY = atof(CPLGetXMLValue(psXMLGCP,"Y","0.0"));
     492              26 :             psGCP->dfGCPZ = atof(CPLGetXMLValue(psXMLGCP,"Z","0.0"));
     493              26 :             nGCPCount++;
     494                 :         }
     495                 :     }
     496                 : 
     497                 : /* -------------------------------------------------------------------- */
     498                 : /*      Get other flags.                                                */
     499                 : /* -------------------------------------------------------------------- */
     500               5 :     nReqOrder = atoi(CPLGetXMLValue(psTree,"Order","3"));
     501               5 :     bReversed = atoi(CPLGetXMLValue(psTree,"Reversed","0"));
     502               5 :     bRefine = atoi(CPLGetXMLValue(psTree,"Refine","0"));
     503               5 :     nMinimumGcps = atoi(CPLGetXMLValue(psTree,"MinimumGcps","6"));
     504               5 :     dfTolerance = atof(CPLGetXMLValue(psTree,"Tolerance","1.0"));
     505                 : 
     506                 : /* -------------------------------------------------------------------- */
     507                 : /*      Generate transformation.                                        */
     508                 : /* -------------------------------------------------------------------- */
     509               5 :     if(bRefine)
     510                 :     {
     511               1 :         pResult = GDALCreateGCPRefineTransformer( nGCPCount, pasGCPList, nReqOrder, 
     512                 :                                         bReversed, dfTolerance, nMinimumGcps );
     513                 :     }
     514                 :     else
     515                 :     {
     516               4 :         pResult = GDALCreateGCPTransformer( nGCPCount, pasGCPList, nReqOrder, 
     517                 :                                         bReversed );
     518                 :     }
     519                 :     
     520                 : /* -------------------------------------------------------------------- */
     521                 : /*      Cleanup GCP copy.                                               */
     522                 : /* -------------------------------------------------------------------- */
     523               5 :     GDALDeinitGCPs( nGCPCount, pasGCPList );
     524               5 :     CPLFree( pasGCPList );
     525                 : 
     526               5 :     return pResult;
     527                 : }
     528                 : 
     529                 : /************************************************************************/
     530                 : /* ==================================================================== */
     531                 : /*      Everything below this point derived from the CRS.C from GRASS.  */
     532                 : /* ==================================================================== */
     533                 : /************************************************************************/
     534                 : 
     535                 : 
     536                 : /* STRUCTURE FOR USE INTERNALLY WITH THESE FUNCTIONS.  THESE FUNCTIONS EXPECT
     537                 :    SQUARE MATRICES SO ONLY ONE VARIABLE IS GIVEN (N) FOR THE MATRIX SIZE */
     538                 : 
     539                 : struct MATRIX
     540                 : {
     541                 :     int     n;     /* SIZE OF THIS MATRIX (N x N) */
     542                 :     double *v;
     543                 : };
     544                 : 
     545                 : /* CALCULATE OFFSET INTO ARRAY BASED ON R/C */
     546                 : 
     547                 : #define M(row,col) m->v[(((row)-1)*(m->n))+(col)-1]
     548                 : 
     549                 : 
     550                 : #define MSUCCESS     1 /* SUCCESS */
     551                 : #define MNPTERR      0 /* NOT ENOUGH POINTS */
     552                 : #define MUNSOLVABLE -1 /* NOT SOLVABLE */
     553                 : #define MMEMERR     -2 /* NOT ENOUGH MEMORY */
     554                 : #define MPARMERR    -3 /* PARAMETER ERROR */
     555                 : #define MINTERR     -4 /* INTERNAL ERROR */
     556                 : 
     557                 : /***************************************************************************/
     558                 : /*
     559                 :     FUNCTION PROTOTYPES FOR STATIC (INTERNAL) FUNCTIONS
     560                 : */
     561                 : /***************************************************************************/
     562                 : 
     563                 : static int calccoef(struct Control_Points *,double *,double *,int);
     564                 : static int calcls(struct Control_Points *,struct MATRIX *,
     565                 :                   double *,double *,double *,double *);
     566                 : static int exactdet(struct Control_Points *,struct MATRIX *,
     567                 :                     double *,double *,double *,double *);
     568                 : static int solvemat(struct MATRIX *,double *,double *,double *,double *);
     569                 : static double term(int,double,double);
     570                 : 
     571                 : /***************************************************************************/
     572                 : /*
     573                 :     TRANSFORM A SINGLE COORDINATE PAIR.
     574                 : */
     575                 : /***************************************************************************/
     576                 : 
     577                 : static int 
     578           90883 : CRS_georef (
     579                 :     double e1,  /* EASTINGS TO BE TRANSFORMED */
     580                 :     double n1,  /* NORTHINGS TO BE TRANSFORMED */
     581                 :     double *e,  /* EASTINGS TO BE TRANSFORMED */
     582                 :     double *n,  /* NORTHINGS TO BE TRANSFORMED */
     583                 :     double E[], /* EASTING COEFFICIENTS */
     584                 :     double N[], /* NORTHING COEFFICIENTS */
     585                 :     int order  /* ORDER OF TRANSFORMATION TO BE PERFORMED, MUST MATCH THE
     586                 :                ORDER USED TO CALCULATE THE COEFFICIENTS */
     587                 : )
     588                 :   {
     589                 :   double e3, e2n, en2, n3, e2, en, n2;
     590                 : 
     591           90883 :   switch(order)
     592                 :     {
     593                 :     case 1:
     594                 : 
     595           88739 :       *e = E[0] + E[1] * e1 + E[2] * n1;
     596           88739 :       *n = N[0] + N[1] * e1 + N[2] * n1;
     597           88739 :       break;
     598                 : 
     599                 :     case 2:
     600                 : 
     601            2144 :       e2 = e1 * e1;
     602            2144 :       n2 = n1 * n1;
     603            2144 :       en = e1 * n1;
     604                 : 
     605            4288 :       *e = E[0]      + E[1] * e1 + E[2] * n1 +
     606            2144 :            E[3] * e2 + E[4] * en + E[5] * n2;
     607            4288 :       *n = N[0]      + N[1] * e1 + N[2] * n1 +
     608            2144 :            N[3] * e2 + N[4] * en + N[5] * n2;
     609            2144 :       break;
     610                 : 
     611                 :     case 3:
     612                 : 
     613               0 :       e2  = e1 * e1;
     614               0 :       en  = e1 * n1;
     615               0 :       n2  = n1 * n1;
     616               0 :       e3  = e1 * e2;
     617               0 :       e2n = e2 * n1;
     618               0 :       en2 = e1 * n2;
     619               0 :       n3  = n1 * n2;
     620                 : 
     621               0 :       *e = E[0]      +
     622               0 :            E[1] * e1 + E[2] * n1  +
     623               0 :            E[3] * e2 + E[4] * en  + E[5] * n2  +
     624               0 :            E[6] * e3 + E[7] * e2n + E[8] * en2 + E[9] * n3;
     625               0 :       *n = N[0]      +
     626               0 :            N[1] * e1 + N[2] * n1  +
     627               0 :            N[3] * e2 + N[4] * en  + N[5] * n2  +
     628               0 :            N[6] * e3 + N[7] * e2n + N[8] * en2 + N[9] * n3;
     629               0 :       break;
     630                 : 
     631                 :     default:
     632                 : 
     633               0 :       return(MPARMERR);
     634                 :     }
     635                 : 
     636           90883 :   return(MSUCCESS);
     637                 :   }
     638                 : 
     639                 : /***************************************************************************/
     640                 : /*
     641                 :     COMPUTE THE GEOREFFERENCING COEFFICIENTS BASED ON A SET OF CONTROL POINTS
     642                 : */
     643                 : /***************************************************************************/
     644                 : 
     645                 : static int 
     646              29 : CRS_compute_georef_equations (struct Control_Points *cp, 
     647                 :                                       double E12[], double N12[], 
     648                 :                                       double E21[], double N21[], 
     649                 :                                       int order)
     650                 : {
     651                 :     double *tempptr;
     652                 :     int status;
     653                 : 
     654              29 :     if(order < 1 || order > MAXORDER)
     655               0 :         return(MPARMERR);
     656                 : 
     657                 :     /* CALCULATE THE FORWARD TRANSFORMATION COEFFICIENTS */
     658                 : 
     659              29 :     status = calccoef(cp,E12,N12,order);
     660              29 :     if(status != MSUCCESS)
     661               0 :         return(status);
     662                 : 
     663                 :     /* SWITCH THE 1 AND 2 EASTING AND NORTHING ARRAYS */
     664                 : 
     665              29 :     tempptr = cp->e1;
     666              29 :     cp->e1 = cp->e2;
     667              29 :     cp->e2 = tempptr;
     668              29 :     tempptr = cp->n1;
     669              29 :     cp->n1 = cp->n2;
     670              29 :     cp->n2 = tempptr;
     671                 : 
     672                 :     /* CALCULATE THE BACKWARD TRANSFORMATION COEFFICIENTS */
     673                 : 
     674              29 :     status = calccoef(cp,E21,N21,order);
     675                 : 
     676                 :     /* SWITCH THE 1 AND 2 EASTING AND NORTHING ARRAYS BACK */
     677                 : 
     678              29 :     tempptr = cp->e1;
     679              29 :     cp->e1 = cp->e2;
     680              29 :     cp->e2 = tempptr;
     681              29 :     tempptr = cp->n1;
     682              29 :     cp->n1 = cp->n2;
     683              29 :     cp->n2 = tempptr;
     684                 : 
     685              29 :     return(status);
     686                 : }
     687                 : 
     688                 : /***************************************************************************/
     689                 : /*
     690                 :     COMPUTE THE GEOREFFERENCING COEFFICIENTS BASED ON A SET OF CONTROL POINTS
     691                 : */
     692                 : /***************************************************************************/
     693                 : 
     694                 : static int 
     695              58 : calccoef (struct Control_Points *cp, double E[], double N[], int order)
     696                 : {
     697                 :     struct MATRIX m;
     698                 :     double *a;
     699                 :     double *b;
     700                 :     int numactive;   /* NUMBER OF ACTIVE CONTROL POINTS */
     701                 :     int status, i;
     702                 : 
     703                 :     /* CALCULATE THE NUMBER OF VALID CONTROL POINTS */
     704                 : 
     705             320 :     for(i = numactive = 0 ; i < cp->count ; i++)
     706                 :     {
     707             262 :         if(cp->status[i] > 0)
     708             262 :             numactive++;
     709                 :     }
     710                 : 
     711                 :     /* CALCULATE THE MINIMUM NUMBER OF CONTROL POINTS NEEDED TO DETERMINE
     712                 :        A TRANSFORMATION OF THIS ORDER */
     713                 : 
     714              58 :     m.n = ((order + 1) * (order + 2)) / 2;
     715                 : 
     716              58 :     if(numactive < m.n)
     717               0 :         return(MNPTERR);
     718                 : 
     719                 :     /* INITIALIZE MATRIX */
     720                 : 
     721              58 :     m.v = (double *)CPLCalloc(m.n*m.n,sizeof(double));
     722              58 :     if(m.v == NULL)
     723                 :     {
     724               0 :         return(MMEMERR);
     725                 :     }
     726              58 :     a = (double *)CPLCalloc(m.n,sizeof(double));
     727              58 :     if(a == NULL)
     728                 :     {
     729               0 :         CPLFree((char *)m.v);
     730               0 :         return(MMEMERR);
     731                 :     }
     732              58 :     b = (double *)CPLCalloc(m.n,sizeof(double));
     733              58 :     if(b == NULL)
     734                 :     {
     735               0 :         CPLFree((char *)m.v);
     736               0 :         CPLFree((char *)a);
     737               0 :         return(MMEMERR);
     738                 :     }
     739                 : 
     740              58 :     if(numactive == m.n)
     741               8 :         status = exactdet(cp,&m,a,b,E,N);
     742                 :     else
     743              50 :         status = calcls(cp,&m,a,b,E,N);
     744                 : 
     745              58 :     CPLFree((char *)m.v);
     746              58 :     CPLFree((char *)a);
     747              58 :     CPLFree((char *)b);
     748                 : 
     749              58 :     return(status);
     750                 : }
     751                 : 
     752                 : /***************************************************************************/
     753                 : /*
     754                 :     CALCULATE THE TRANSFORMATION COEFFICIENTS WITH EXACTLY THE MINIMUM
     755                 :     NUMBER OF CONTROL POINTS REQUIRED FOR THIS TRANSFORMATION.
     756                 : */
     757                 : /***************************************************************************/
     758                 : 
     759               8 : static int exactdet (
     760                 :     struct Control_Points *cp,
     761                 :     struct MATRIX *m,
     762                 :     double a[],
     763                 :     double b[],
     764                 :     double E[],     /* EASTING COEFFICIENTS */
     765                 :     double N[]     /* NORTHING COEFFICIENTS */
     766                 : )
     767                 :   {
     768                 :   int pntnow, currow, j;
     769                 : 
     770               8 :   currow = 1;
     771              32 :   for(pntnow = 0 ; pntnow < cp->count ; pntnow++)
     772                 :     {
     773              24 :     if(cp->status[pntnow] > 0)
     774                 :       {
     775                 :       /* POPULATE MATRIX M */
     776                 : 
     777              96 :       for(j = 1 ; j <= m->n ; j++)
     778                 :         {
     779              72 :         M(currow,j) = term(j,cp->e1[pntnow],cp->n1[pntnow]);
     780                 :         }
     781                 : 
     782                 :       /* POPULATE MATRIX A AND B */
     783                 : 
     784              24 :       a[currow-1] = cp->e2[pntnow];
     785              24 :       b[currow-1] = cp->n2[pntnow];
     786                 : 
     787              24 :       currow++;
     788                 :       }
     789                 :     }
     790                 : 
     791               8 :   if(currow - 1 != m->n)
     792               0 :     return(MINTERR);
     793                 : 
     794               8 :   return(solvemat(m,a,b,E,N));
     795                 :   }
     796                 : 
     797                 : /***************************************************************************/
     798                 : /*
     799                 :     CALCULATE THE TRANSFORMATION COEFFICIENTS WITH MORE THAN THE MINIMUM
     800                 :     NUMBER OF CONTROL POINTS REQUIRED FOR THIS TRANSFORMATION.  THIS
     801                 :     ROUTINE USES THE LEAST SQUARES METHOD TO COMPUTE THE COEFFICIENTS.
     802                 : */
     803                 : /***************************************************************************/
     804                 : 
     805              50 : static int calcls (
     806                 :     struct Control_Points *cp,
     807                 :     struct MATRIX *m,
     808                 :     double a[],
     809                 :     double b[],
     810                 :     double E[],     /* EASTING COEFFICIENTS */
     811                 :     double N[]     /* NORTHING COEFFICIENTS */
     812                 : )
     813                 : {
     814              50 :     int i, j, n, numactive = 0;
     815                 : 
     816                 :     /* INITIALIZE THE UPPER HALF OF THE MATRIX AND THE TWO COLUMN VECTORS */
     817                 : 
     818             206 :     for(i = 1 ; i <= m->n ; i++)
     819                 :     {
     820             486 :         for(j = i ; j <= m->n ; j++)
     821             330 :             M(i,j) = 0.0;
     822             156 :         a[i-1] = b[i-1] = 0.0;
     823                 :     }
     824                 : 
     825                 :     /* SUM THE UPPER HALF OF THE MATRIX AND THE COLUMN VECTORS ACCORDING TO
     826                 :        THE LEAST SQUARES METHOD OF SOLVING OVER DETERMINED SYSTEMS */
     827                 : 
     828             288 :     for(n = 0 ; n < cp->count ; n++)
     829                 :     {
     830             238 :         if(cp->status[n] > 0)
     831                 :         {
     832             238 :             numactive++;
     833            1000 :             for(i = 1 ; i <= m->n ; i++)
     834                 :             {
     835            2430 :                 for(j = i ; j <= m->n ; j++)
     836            1668 :                     M(i,j) += term(i,cp->e1[n],cp->n1[n]) * term(j,cp->e1[n],cp->n1[n]);
     837                 : 
     838             762 :                 a[i-1] += cp->e2[n] * term(i,cp->e1[n],cp->n1[n]);
     839             762 :                 b[i-1] += cp->n2[n] * term(i,cp->e1[n],cp->n1[n]);
     840                 :             }
     841                 :         }
     842                 :     }
     843                 : 
     844              50 :     if(numactive <= m->n)
     845               0 :         return(MINTERR);
     846                 : 
     847                 :     /* TRANSPOSE VALUES IN UPPER HALF OF M TO OTHER HALF */
     848                 : 
     849             156 :     for(i = 2 ; i <= m->n ; i++)
     850                 :     {
     851             280 :         for(j = 1 ; j < i ; j++)
     852             174 :             M(i,j) = M(j,i);
     853                 :     }
     854                 : 
     855              50 :     return(solvemat(m,a,b,E,N));
     856                 : }
     857                 : 
     858                 : /***************************************************************************/
     859                 : /*
     860                 :     CALCULATE THE X/Y TERM BASED ON THE TERM NUMBER
     861                 : 
     862                 : ORDER\TERM   1    2    3    4    5    6    7    8    9   10
     863                 :   1        e0n0 e1n0 e0n1
     864                 :   2        e0n0 e1n0 e0n1 e2n0 e1n1 e0n2
     865                 :   3        e0n0 e1n0 e0n1 e2n0 e1n1 e0n2 e3n0 e2n1 e1n2 e0n3
     866                 : */
     867                 : /***************************************************************************/
     868                 : 
     869            4932 : static double term (int term, double e, double n)
     870                 : {
     871            4932 :     switch(term)
     872                 :     {
     873            1500 :       case  1: return((double)1.0);
     874            1500 :       case  2: return((double)e);
     875            1500 :       case  3: return((double)n);
     876             144 :       case  4: return((double)(e*e));
     877             144 :       case  5: return((double)(e*n));
     878             144 :       case  6: return((double)(n*n));
     879               0 :       case  7: return((double)(e*e*e));
     880               0 :       case  8: return((double)(e*e*n));
     881               0 :       case  9: return((double)(e*n*n));
     882               0 :       case 10: return((double)(n*n*n));
     883                 :     }
     884               0 :     return((double)0.0);
     885                 : }
     886                 : 
     887                 : /***************************************************************************/
     888                 : /*
     889                 :     SOLVE FOR THE 'E' AND 'N' COEFFICIENTS BY USING A SOMEWHAT MODIFIED
     890                 :     GAUSSIAN ELIMINATION METHOD.
     891                 : 
     892                 :     | M11 M12 ... M1n | | E0   |   | a0   |
     893                 :     | M21 M22 ... M2n | | E1   | = | a1   |
     894                 :     |  .   .   .   .  | | .    |   | .    |
     895                 :     | Mn1 Mn2 ... Mnn | | En-1 |   | an-1 |
     896                 : 
     897                 :     and
     898                 : 
     899                 :     | M11 M12 ... M1n | | N0   |   | b0   |
     900                 :     | M21 M22 ... M2n | | N1   | = | b1   |
     901                 :     |  .   .   .   .  | | .    |   | .    |
     902                 :     | Mn1 Mn2 ... Mnn | | Nn-1 |   | bn-1 |
     903                 : */
     904                 : /***************************************************************************/
     905                 : 
     906              58 : static int solvemat (struct MATRIX *m,
     907                 :   double a[], double b[], double E[], double N[])
     908                 : {
     909                 :     int i, j, i2, j2, imark;
     910                 :     double factor, temp;
     911                 :     double  pivot;  /* ACTUAL VALUE OF THE LARGEST PIVOT CANDIDATE */
     912                 : 
     913             238 :     for(i = 1 ; i <= m->n ; i++)
     914                 :     {
     915             180 :         j = i;
     916                 : 
     917                 :         /* find row with largest magnitude value for pivot value */
     918                 : 
     919             180 :         pivot = M(i,j);
     920             180 :         imark = i;
     921             378 :         for(i2 = i + 1 ; i2 <= m->n ; i2++)
     922                 :         {
     923             198 :             temp = fabs(M(i2,j));
     924             198 :             if(temp > fabs(pivot))
     925                 :             {
     926             120 :                 pivot = M(i2,j);
     927             120 :                 imark = i2;
     928                 :             }
     929                 :         }
     930                 : 
     931                 :         /* if the pivot is very small then the points are nearly co-linear */
     932                 :         /* co-linear points result in an undefined matrix, and nearly */
     933                 :         /* co-linear points results in a solution with rounding error */
     934                 : 
     935             180 :         if(pivot == 0.0)
     936               0 :             return(MUNSOLVABLE);
     937                 : 
     938                 :         /* if row with highest pivot is not the current row, switch them */
     939                 :  
     940             180 :         if(imark != i)
     941                 :         {
     942             396 :             for(j2 = 1 ; j2 <= m->n ; j2++)
     943                 :             {
     944             303 :                 temp = M(imark,j2);
     945             303 :                 M(imark,j2) = M(i,j2);
     946             303 :                 M(i,j2) = temp;
     947                 :             }
     948                 : 
     949              93 :             temp = a[imark-1];
     950              93 :             a[imark-1] = a[i-1];
     951              93 :             a[i-1] = temp;
     952                 : 
     953              93 :             temp = b[imark-1];
     954              93 :             b[imark-1] = b[i-1];
     955              93 :             b[i-1] = temp;
     956                 :         }
     957                 : 
     958                 :         /* compute zeros above and below the pivot, and compute
     959                 :            values for the rest of the row as well */
     960                 : 
     961             756 :         for(i2 = 1 ; i2 <= m->n ; i2++)
     962                 :         {
     963             576 :             if(i2 != i)
     964                 :             {
     965             396 :                 factor = M(i2,j) / pivot;
     966            1278 :                 for(j2 = j ; j2 <= m->n ; j2++)
     967             882 :                     M(i2,j2) -= factor * M(i,j2);
     968             396 :                 a[i2-1] -= factor * a[i-1];
     969             396 :                 b[i2-1] -= factor * b[i-1];
     970                 :             }
     971                 :         }
     972                 :     }
     973                 : 
     974                 :     /* SINCE ALL OTHER VALUES IN THE MATRIX ARE ZERO NOW, CALCULATE THE
     975                 :        COEFFICIENTS BY DIVIDING THE COLUMN VECTORS BY THE DIAGONAL VALUES. */
     976                 : 
     977             238 :     for(i = 1 ; i <= m->n ; i++)
     978                 :     {
     979             180 :         E[i-1] = a[i-1] / M(i,i);
     980             180 :         N[i-1] = b[i-1] / M(i,i);
     981                 :     }
     982                 : 
     983              58 :     return(MSUCCESS);
     984                 : }
     985                 : 
     986                 : /***************************************************************************/
     987                 : /*
     988                 :   DETECTS THE WORST OUTLIER IN THE GCP LIST AND RETURNS THE INDEX OF THE
     989                 :   OUTLIER.
     990                 :   
     991                 :   THE WORST OUTLIER IS CALCULATED BASED ON THE CONTROL POINTS, COEFFICIENTS
     992                 :   AND THE ALLOWED TOLERANCE:
     993                 :   
     994                 :   sampleAdj = a0 + a1*sample + a2*line + a3*line*sample
     995                 :   lineAdj = b0 + b1*sample + b2*line + b3*line*sample
     996                 :   
     997                 :   WHERE sampleAdj AND lineAdj ARE CORRELATED GCPS
     998                 :   
     999                 :   [residualSample] = [A1][sampleCoefficients] - [b1]
    1000                 :   [residualLine] = [A2][lineCoefficients] - [b2]
    1001                 :   
    1002                 :   sampleResidual^2 = sum( [residualSample]^2 )
    1003                 :   lineResidual^2 = sum( [lineSample]^2 )
    1004                 :   
    1005                 :   residuals(i) = squareRoot( residualSample(i)^2 + residualLine(i)^2 )
    1006                 :   
    1007                 :   THE GCP WITH THE GREATEST DISTANCE residual(i) GREATER THAN THE TOLERANCE WILL
    1008                 :   CONSIDERED THE WORST OUTLIER.
    1009                 :   
    1010                 :   IF NO OUTLIER CAN BE FOUND, -1 WILL BE RETURNED.
    1011                 : */
    1012                 : /***************************************************************************/
    1013               3 : static int worst_outlier(struct Control_Points *cp, double E[], double N[], double dfTolerance)
    1014                 : {
    1015                 :     double *padfResiduals;
    1016                 :     int nI, nIndex;
    1017                 :     double dfDifference, dfSampleResidual, dfLineResidual, dfSampleRes, dfLineRes, dfCurrentDifference;
    1018                 :     double dfE1, dfN1, dfE2, dfN2, dfEn;
    1019                 :   
    1020               3 :     padfResiduals = (double *) CPLCalloc(sizeof(double),cp->count);
    1021               3 :     dfSampleResidual = 0.0;
    1022               3 :     dfLineResidual = 0.0;
    1023                 :   
    1024              30 :     for(nI = 0; nI < cp->count; nI++)
    1025                 :     {
    1026              27 :         dfE1 = cp->e1[nI];
    1027              27 :         dfN1 = cp->n1[nI];
    1028              27 :         dfE2 = dfE1 * dfE1;
    1029              27 :         dfN2 = dfN1 * dfN1;
    1030              27 :         dfEn = dfE1 * dfN1;
    1031                 : 
    1032              27 :         dfSampleRes = E[0] + E[1] * dfE1 + E[2] * dfN1 + E[3] * dfE2 + E[4] * dfEn + E[5] * dfN2 - cp->e2[nI];
    1033              27 :         dfLineRes = N[0] + N[1] * dfE1 + N[2] * dfN1 + N[3] * dfE2 + N[4] * dfEn + N[5] * dfN2 - cp->n2[nI];
    1034                 :     
    1035              27 :         dfSampleResidual += dfSampleRes*dfSampleRes;
    1036              27 :         dfLineResidual += dfLineRes*dfLineRes;
    1037                 :     
    1038              27 :         padfResiduals[nI] = sqrt(dfSampleRes*dfSampleRes + dfLineRes*dfLineRes);
    1039                 :     }
    1040                 :   
    1041               3 :     nIndex = -1;
    1042               3 :     dfDifference = -1.0;
    1043              30 :     for(nI = 0; nI < cp->count; nI++)
    1044                 :     {
    1045              27 :         dfCurrentDifference = padfResiduals[nI];
    1046              27 :         if(fabs(dfCurrentDifference) < 1.19209290E-07F /*FLT_EPSILON*/)
    1047                 :         {
    1048               8 :             dfCurrentDifference = 0.0;
    1049                 :         }
    1050              27 :         if(dfCurrentDifference > dfDifference && dfCurrentDifference >= dfTolerance)
    1051                 :         {
    1052               6 :             dfDifference = dfCurrentDifference;
    1053               6 :             nIndex = nI;
    1054                 :         }
    1055                 :     }
    1056               3 :     CPLFree( padfResiduals );
    1057               3 :     return nIndex;
    1058                 : }
    1059                 : 
    1060                 : /***************************************************************************/
    1061                 : /*
    1062                 :   REMOVES THE WORST OUTLIERS ITERATIVELY UNTIL THE MINIMUM NUMBER OF GCPS
    1063                 :   ARE REACHED OR NO OUTLIERS CAN BE DETECTED.
    1064                 :   
    1065                 :   1. WE CALCULATE THE COEFFICIENTS FOR ALL THE GCPS.
    1066                 :   2. THE GCP LIST WILL BE SCANED TO DETERMINE THE WORST OUTLIER USING
    1067                 :      THE CALCULATED COEFFICIENTS.
    1068                 :   3. THE WORST OUTLIER WILL BE REMOVED FROM THE GCP LIST.
    1069                 :   4. THE COEFFICIENTS WILL BE RECALCULATED WITHOUT THE WORST OUTLIER.
    1070                 :   5. STEP 1 TO 4 ARE EXECUTED UNTIL THE MINIMUM NUMBER OF GCPS WERE REACHED
    1071                 :      OR IF NO GCP IS CONSIDERED AN OUTLIER WITH THE PASSED TOLERANCE.
    1072                 : */
    1073                 : /***************************************************************************/
    1074               1 : static int remove_outliers( GCPTransformInfo *psInfo )
    1075                 : {
    1076                 :     double *padfGeoX, *padfGeoY, *padfRasterX, *padfRasterY;
    1077                 :     int *panStatus;
    1078                 :     int nI, nCRSresult, nGCPCount, nMinimumGcps, nReqOrder;
    1079                 :     double dfTolerance;
    1080                 :     struct Control_Points sPoints;
    1081                 :     
    1082               1 :     nGCPCount = psInfo->nGCPCount;
    1083               1 :     nMinimumGcps = psInfo->nMinimumGcps;
    1084               1 :     nReqOrder = psInfo->nOrder;
    1085               1 :     dfTolerance = psInfo->dfTolerance;
    1086                 :     
    1087               1 :     padfGeoX = (double *) CPLCalloc(sizeof(double),nGCPCount);
    1088               1 :     padfGeoY = (double *) CPLCalloc(sizeof(double),nGCPCount);
    1089               1 :     padfRasterX = (double *) CPLCalloc(sizeof(double),nGCPCount);
    1090               1 :     padfRasterY = (double *) CPLCalloc(sizeof(double),nGCPCount);
    1091               1 :     panStatus = (int *) CPLCalloc(sizeof(int),nGCPCount);
    1092                 :     
    1093              11 :     for( nI = 0; nI < nGCPCount; nI++ )
    1094                 :     {
    1095              10 :         panStatus[nI] = 1;
    1096              10 :         padfGeoX[nI] = psInfo->pasGCPList[nI].dfGCPX;
    1097              10 :         padfGeoY[nI] = psInfo->pasGCPList[nI].dfGCPY;
    1098              10 :         padfRasterX[nI] = psInfo->pasGCPList[nI].dfGCPPixel;
    1099              10 :         padfRasterY[nI] = psInfo->pasGCPList[nI].dfGCPLine;
    1100                 :     }
    1101                 : 
    1102               1 :     sPoints.count = nGCPCount;
    1103               1 :     sPoints.e1 = padfRasterX;
    1104               1 :     sPoints.n1 = padfRasterY;
    1105               1 :     sPoints.e2 = padfGeoX;
    1106               1 :     sPoints.n2 = padfGeoY;
    1107               1 :     sPoints.status = panStatus;
    1108                 :   
    1109               1 :     nCRSresult = CRS_compute_georef_equations( &sPoints,
    1110                 :                                       psInfo->adfToGeoX, psInfo->adfToGeoY,
    1111                 :                                       psInfo->adfFromGeoX, psInfo->adfFromGeoY,
    1112                 :                                       nReqOrder );
    1113                 : 
    1114               4 :     while(sPoints.count > nMinimumGcps)
    1115                 :     {
    1116                 :         int nIndex;
    1117                 : 
    1118               3 :         nIndex = worst_outlier(&sPoints, psInfo->adfFromGeoX, psInfo->adfFromGeoY, dfTolerance);
    1119                 : 
    1120                 :         //If no outliers were detected, stop the GCP elimination
    1121               3 :         if(nIndex == -1)
    1122                 :         {
    1123               1 :             break;
    1124                 :         }
    1125                 : 
    1126               2 :         CPLFree(psInfo->pasGCPList[nIndex].pszId);
    1127               2 :         CPLFree(psInfo->pasGCPList[nIndex].pszInfo);
    1128                 : 
    1129               8 :         for( nI = nIndex; nI < sPoints.count - 1; nI++ )
    1130                 :         {
    1131               6 :             sPoints.e1[nI] = sPoints.e1[nI + 1];
    1132               6 :             sPoints.n1[nI] = sPoints.n1[nI + 1];
    1133               6 :             sPoints.e2[nI] = sPoints.e2[nI + 1];
    1134               6 :             sPoints.n2[nI] = sPoints.n2[nI + 1];
    1135               6 :             psInfo->pasGCPList[nI].pszId = psInfo->pasGCPList[nI + 1].pszId;
    1136               6 :             psInfo->pasGCPList[nI].pszInfo = psInfo->pasGCPList[nI + 1].pszInfo;
    1137                 :         }
    1138                 : 
    1139               2 :         sPoints.count = sPoints.count - 1;
    1140                 : 
    1141               2 :         nCRSresult = CRS_compute_georef_equations( &sPoints,
    1142                 :                                       psInfo->adfToGeoX, psInfo->adfToGeoY,
    1143                 :                                       psInfo->adfFromGeoX, psInfo->adfFromGeoY,
    1144                 :                                       nReqOrder );
    1145                 :     }
    1146                 : 
    1147               9 :     for( nI = 0; nI < sPoints.count; nI++ )
    1148                 :     {
    1149               8 :         psInfo->pasGCPList[nI].dfGCPX = sPoints.e2[nI];
    1150               8 :         psInfo->pasGCPList[nI].dfGCPY = sPoints.n2[nI];
    1151               8 :         psInfo->pasGCPList[nI].dfGCPPixel = sPoints.e1[nI];
    1152               8 :         psInfo->pasGCPList[nI].dfGCPLine = sPoints.n1[nI];
    1153                 :     }
    1154               1 :     psInfo->nGCPCount = sPoints.count;
    1155                 :     
    1156               1 :     CPLFree( sPoints.e1 );
    1157               1 :     CPLFree( sPoints.n1 );
    1158               1 :     CPLFree( sPoints.e2 );
    1159               1 :     CPLFree( sPoints.n2 );
    1160               1 :     CPLFree( sPoints.status );
    1161               1 :     return nCRSresult;
    1162                 : }

Generated by: LCOV version 1.7