Ticket #4143: gcp_refinement_22652.diff

File gcp_refinement_22652.diff, 24.9 KB (added by goocreations, 5 years ago)

Updated GCP refinement

  • apps/gdalwarp.cpp

     
    7878    [-r resampling_method] [-wm memory_in_mb] [-multi] [-q]
    7979    [-cutline datasource] [-cl layer] [-cwhere expression]
    8080    [-csql statement] [-cblend dist_in_pixels] [-crop_to_cutline]
    81     [-of format] [-co "NAME=VALUE"]* [-overwrite]
     81    [-of format] [-co "NAME=VALUE"]* [-overwrite] [-refine_gcps tolerance minimum_gcps]
    8282    srcfile* dstfile
    8383\endverbatim
    8484
     
    163163the output format driver. Multiple <b>-co</b> options may be listed. See
    164164format specific documentation for legal creation options for each format.
    165165</dd>
     166<dt> <b>-refine_gcps</b> <em>tolerance minimum_gcps</em>:</dt><dd> refines the GCPs by automatically eliminating outliers.
     167Outliers will be eliminated until minimum_gcps are left or when no outliers can be detected.
     168The tolerance is passed to adjust when a GCP will be eliminated.
     169Not that GCP refinement only works with polynomial interpolation.
     170The tolerance is in pixel units if no projection is available, otherwise it is in SRS units.
     171If minimum_gcps is not provided, the minimum GCPs according to the polynomial model is used.</dd>
    166172
    167173<dt> <b>-cutline</b> <em>datasource</em>:</dt><dd>Enable use of a blend cutline from the name OGR support datasource.</dd>
    168174<dt> <b>-cl</b> <em>layername</em>:</dt><dd>Select the named layer from the
     
    231237        "    [-r resampling_method] [-wm memory_in_mb] [-multi] [-q]\n"
    232238        "    [-cutline datasource] [-cl layer] [-cwhere expression]\n"
    233239        "    [-csql statement] [-cblend dist_in_pixels] [-crop_to_cutline]\n"
    234         "    [-of format] [-co \"NAME=VALUE\"]* [-overwrite]\n"
     240        "    [-of format] [-co \"NAME=VALUE\"]* [-overwrite] [-refine_gcps tolerance minimum_gcps]\n"
    235241        "    srcfile* dstfile\n"
    236242        "\n"
    237243        "Available resampling methods:\n"
     
    388394        {
    389395            papszTO = CSLSetNameValue( papszTO, "MAX_GCP_ORDER", argv[++i] );
    390396        }
     397        else if( EQUAL(argv[i],"-refine_gcps") && i < argc-2 )
     398        {
     399            papszTO = CSLSetNameValue( papszTO, "REFINE_TOLERANCE", argv[++i] );
     400            if(atof(argv[i]) < 0)
     401            {
     402                printf( "The tolerance for -refine_gcps may not be negative\n");
     403                Usage();
     404            }
     405            if(atoi(argv[i+1]) >= 0 && isdigit(argv[i+1][0]))
     406            {
     407                papszTO = CSLSetNameValue( papszTO, "REFINE_MINIMUM_GCPS", argv[++i] );
     408            }
     409            else
     410            {
     411                papszTO = CSLSetNameValue( papszTO, "REFINE_MINIMUM_GCPS", "-1" );
     412            }
     413        }
    391414        else if( EQUAL(argv[i],"-tps") )
    392415        {
    393416            papszTO = CSLSetNameValue( papszTO, "METHOD", "GCP_TPS" );
     
    774797        CSLDestroy( papszCreateOptions );
    775798        papszCreateOptions = NULL;
    776799    }
    777 
     800 
    778801    if( hDstDS == NULL )
    779802        exit( 1 );
    780803
     
    799822/*      Check that there's at least one raster band                     */
    800823/* -------------------------------------------------------------------- */
    801824        if ( GDALGetRasterCount(hSrcDS) == 0 )
    802         {
     825        {     
    803826            fprintf(stderr, "Input file %s has no raster bands.\n", papszSrcFiles[iSrc] );
    804827            exit( 1 );
    805828        }
     
    811834/*      Warns if the file has a color table and something more          */
    812835/*      complicated than nearest neighbour resampling is asked          */
    813836/* -------------------------------------------------------------------- */
    814 
     837 
    815838        if ( eResampleAlg != GRA_NearestNeighbour &&
    816839             GDALGetRasterColorTable(GDALGetRasterBand(hSrcDS, 1)) != NULL)
    817840        {
  • alg/gdaltransformer.cpp

     
    10001000 * <li> SRC_SRS: WKT SRS to be used as an override for hSrcDS.
    10011001 * <li> DST_SRS: WKT SRS to be used as an override for hDstDS.
    10021002 * <li> GCPS_OK: If false, GCPs will not be used, default is TRUE.
     1003 * <li> REFINE_MINIMUM_GCPS: The minimum amount of GCPs that should be available after the refinement.
     1004 * <li> REFINE_TOLERANCE: The tolernace that specifies when a GCP will be eliminated.
    10031005 * <li> MAX_GCP_ORDER: the maximum order to use for GCP derived polynomials if
    10041006 * possible.  The default is to autoselect based on the number of GCPs. 
    10051007 * A value of -1 triggers use of Thin Plate Spline instead of polynomials.
     
    10301032    GDALRPCInfo sRPCInfo;
    10311033    const char *pszMethod = CSLFetchNameValue( papszOptions, "METHOD" );
    10321034    const char *pszValue;
    1033     int nOrder = 0, bGCPUseOK = TRUE;
     1035    int nOrder = 0, bGCPUseOK = TRUE, nMinimumGcps = -1, bRefine = FALSE;
     1036    double dfTolerance = 0.0;
    10341037    const char *pszSrcWKT = CSLFetchNameValue( papszOptions, "SRC_SRS" );
    10351038    const char *pszDstWKT = CSLFetchNameValue( papszOptions, "DST_SRS" );
    10361039
     
    10421045    if( pszValue )
    10431046        bGCPUseOK = CSLTestBoolean(pszValue);
    10441047
     1048    pszValue = CSLFetchNameValue( papszOptions, "REFINE_MINIMUM_GCPS" );
     1049    if( pszValue )
     1050    {
     1051        if( atoi(pszValue) != -1)
     1052            nMinimumGcps = atoi(pszValue);
     1053    }
     1054
     1055    pszValue = CSLFetchNameValue( papszOptions, "REFINE_TOLERANCE" );
     1056    if( pszValue )
     1057    {
     1058        dfTolerance = atof(pszValue);
     1059        bRefine = TRUE;
     1060    }
     1061
    10451062/* -------------------------------------------------------------------- */
    10461063/*      Initialize the transform info.                                  */
    10471064/* -------------------------------------------------------------------- */
     
    10891106             && (pszMethod == NULL || EQUAL(pszMethod,"GCP_POLYNOMIAL") )
    10901107             && GDALGetGCPCount( hSrcDS ) > 0 && nOrder >= 0 )
    10911108    {
    1092         psInfo->pSrcGCPTransformArg =
    1093             GDALCreateGCPTransformer( GDALGetGCPCount( hSrcDS ),
    1094                                       GDALGetGCPs( hSrcDS ), nOrder,
    1095                                       FALSE );
     1109        if(bRefine)
     1110        {
     1111                psInfo->pSrcGCPTransformArg =
     1112                    GDALCreateGCPRefineTransformer( GDALGetGCPCount( hSrcDS ),
     1113                                                    GDALGetGCPs( hSrcDS ), nOrder,
     1114                                                    FALSE, dfTolerance, nMinimumGcps );
     1115        }
     1116        else
     1117        {
     1118            psInfo->pSrcGCPTransformArg =
     1119                GDALCreateGCPTransformer( GDALGetGCPCount( hSrcDS ),
     1120                                          GDALGetGCPs( hSrcDS ), nOrder,
     1121                                          FALSE );
     1122        }
    10961123
    10971124        if( psInfo->pSrcGCPTransformArg == NULL )
    10981125        {
  • alg/gdal_alg.h

     
    165165void CPL_DLL *
    166166GDALCreateGCPTransformer( int nGCPCount, const GDAL_GCP *pasGCPList,
    167167                          int nReqOrder, int bReversed );
     168                         
     169/* GCP based transformer with refinement of the GCPs ... forward is to georef coordinates */
     170void CPL_DLL *
     171GDALCreateGCPRefineTransformer( int nGCPCount, const GDAL_GCP *pasGCPList,
     172                                int nReqOrder, int bReversed, double tolerance, int minimumGcps);
     173                         
    168174void CPL_DLL GDALDestroyGCPTransformer( void *pTransformArg );
    169175int CPL_DLL GDALGCPTransform(
    170176    void *pTransformArg, int bDstToSrc, int nPointCount,
  • alg/gdal_crs.c

     
    2626      Added printout of trnfile. Triggered by BDEBUG.
    2727    Last Update:  1/27/92 Brian J. Buckley
    2828      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     
    2932
    30 
    3133    Copyright (c) 1992, Michigan State University
    3234   
    3335    Permission is hereby granted, free of charge, to any person obtaining a
     
    5557#include "cpl_minixml.h"
    5658#include "cpl_string.h"
    5759
     60
    5861CPL_CVSID("$Id$");
    5962
    6063#define MAXORDER 3
     
    6972    int *status;
    7073};
    7174
     75typedef struct
     76{
     77    GDALTransformerInfo sTI;
     78
     79    double adfToGeoX[20];
     80    double adfToGeoY[20];
     81   
     82    double adfFromGeoX[20];
     83    double adfFromGeoY[20];
     84
     85    int    nOrder;
     86    int    bReversed;
     87
     88    int       nGCPCount;
     89    GDAL_GCP *pasGCPList;
     90    int    bRefine;
     91    int    nMinimumGcps;
     92    double dfTolerance;
     93   
     94} GCPTransformInfo;
     95
    7296CPL_C_START
    7397CPLXMLNode *GDALSerializeGCPTransformer( void *pTransformArg );
    7498void *GDALDeserializeGCPTransformer( CPLXMLNode *psTree );
     
    79103                              double [], double [], int);
    80104static int CRS_compute_georef_equations(struct Control_Points *,
    81105    double [], double [], double [], double [], int);
     106static int remove_outliers(GCPTransformInfo *);
    82107
    83 
    84108static char *CRS_error_message[] = {
    85109    "Failed to compute GCP transform: Not enough points available",
    86110    "Failed to compute GCP transform: Transform is not solvable",
     
    89113    "Failed to compute GCP transform: Internal error"
    90114};
    91115
    92 typedef struct
    93 {
    94     GDALTransformerInfo sTI;
    95116
    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 
    111117/************************************************************************/
    112118/*                      GDALCreateGCPTransformer()                      */
    113119/************************************************************************/
     
    139145 * @return the transform argument or NULL if creation fails.
    140146 */
    141147
    142 void *GDALCreateGCPTransformer( int nGCPCount, const GDAL_GCP *pasGCPList,
    143                                 int nReqOrder, int bReversed )
     148void *GDALCreateGCPTransformerEx( int nGCPCount, const GDAL_GCP *pasGCPList,
     149                                int nReqOrder, int bReversed, int bRefine, double dfTolerance, int nMinimumGcps)
    144150
    145151{
    146152    GCPTransformInfo *psInfo;
     
    162168    psInfo = (GCPTransformInfo *) CPLCalloc(sizeof(GCPTransformInfo),1);
    163169    psInfo->bReversed = bReversed;
    164170    psInfo->nOrder = nReqOrder;
     171    psInfo->bRefine = bRefine;
     172    psInfo->dfTolerance = dfTolerance;
     173    psInfo->nMinimumGcps = nMinimumGcps;
    165174
    166175    psInfo->pasGCPList = GDALDuplicateGCPs( nGCPCount, pasGCPList );
    167176    psInfo->nGCPCount = nGCPCount;
     
    171180    psInfo->sTI.pfnTransform = GDALGCPTransform;
    172181    psInfo->sTI.pfnCleanup = GDALDestroyGCPTransformer;
    173182    psInfo->sTI.pfnSerialize = GDALSerializeGCPTransformer;
    174 
     183   
    175184/* -------------------------------------------------------------------- */
    176 /*      Allocate and initialize the working points list.                */
     185/*      Compute the forward and reverse polynomials.                    */
    177186/* -------------------------------------------------------------------- */
    178     padfGeoX = (double *) CPLCalloc(sizeof(double),nGCPCount);
    179     padfGeoY = (double *) CPLCalloc(sizeof(double),nGCPCount);
    180     padfRasterX = (double *) CPLCalloc(sizeof(double),nGCPCount);
    181     padfRasterY = (double *) CPLCalloc(sizeof(double),nGCPCount);
    182     panStatus = (int *) CPLCalloc(sizeof(int),nGCPCount);
    183    
    184     for( iGCP = 0; iGCP < nGCPCount; iGCP++ )
     187
     188    if(bRefine)
    185189    {
    186         panStatus[iGCP] = 1;
    187         padfGeoX[iGCP] = pasGCPList[iGCP].dfGCPX;
    188         padfGeoY[iGCP] = pasGCPList[iGCP].dfGCPY;
    189         padfRasterX[iGCP] = pasGCPList[iGCP].dfGCPPixel;
    190         padfRasterY[iGCP] = pasGCPList[iGCP].dfGCPLine;
     190        nCRSresult = remove_outliers(psInfo);
    191191    }
    192 
    193     sPoints.count = nGCPCount;
    194     sPoints.e1 = padfRasterX;
    195     sPoints.n1 = padfRasterY;
    196     sPoints.e2 = padfGeoX;
    197     sPoints.n2 = padfGeoY;
    198     sPoints.status = panStatus;
     192    else
     193    {
     194        /* -------------------------------------------------------------------- */
     195        /*      Allocate and initialize the working points list.                */
     196        /* -------------------------------------------------------------------- */
     197        padfGeoX = (double *) CPLCalloc(sizeof(double),nGCPCount);
     198        padfGeoY = (double *) CPLCalloc(sizeof(double),nGCPCount);
     199        padfRasterX = (double *) CPLCalloc(sizeof(double),nGCPCount);
     200        padfRasterY = (double *) CPLCalloc(sizeof(double),nGCPCount);
     201        panStatus = (int *) CPLCalloc(sizeof(int),nGCPCount);
    199202   
    200 /* -------------------------------------------------------------------- */
    201 /*      Compute the forward and reverse polynomials.                    */
    202 /* -------------------------------------------------------------------- */
    203     nCRSresult = CRS_compute_georef_equations( &sPoints,
    204                                       psInfo->adfToGeoX, psInfo->adfToGeoY,
    205                                       psInfo->adfFromGeoX, psInfo->adfFromGeoY,
    206                                       nReqOrder );
     203        for( iGCP = 0; iGCP < nGCPCount; iGCP++ )
     204        {
     205            panStatus[iGCP] = 1;
     206            padfGeoX[iGCP] = pasGCPList[iGCP].dfGCPX;
     207            padfGeoY[iGCP] = pasGCPList[iGCP].dfGCPY;
     208            padfRasterX[iGCP] = pasGCPList[iGCP].dfGCPPixel;
     209            padfRasterY[iGCP] = pasGCPList[iGCP].dfGCPLine;
     210        }
    207211
    208     CPLFree( padfGeoX );
    209     CPLFree( padfGeoY );
    210     CPLFree( padfRasterX );
    211     CPLFree( padfRasterY );
    212     CPLFree( panStatus );
     212        sPoints.count = nGCPCount;
     213        sPoints.e1 = padfRasterX;
     214        sPoints.n1 = padfRasterY;
     215        sPoints.e2 = padfGeoX;
     216        sPoints.n2 = padfGeoY;
     217        sPoints.status = panStatus;
     218        nCRSresult = CRS_compute_georef_equations( &sPoints,
     219                                                psInfo->adfToGeoX, psInfo->adfToGeoY,
     220                                                psInfo->adfFromGeoX, psInfo->adfFromGeoY,
     221                                                nReqOrder );
     222        CPLFree( padfGeoX );
     223        CPLFree( padfGeoY );
     224        CPLFree( padfRasterX );
     225        CPLFree( padfRasterY );
     226        CPLFree( panStatus );
     227    }
    213228
    214229    if (nCRSresult != 1)
    215230    {
     
    223238    }
    224239}
    225240
     241void *GDALCreateGCPTransformer( int nGCPCount, const GDAL_GCP *pasGCPList,
     242                                int nReqOrder, int bReversed )
     243
     244{
     245    return GDALCreateGCPTransformerEx(nGCPCount, pasGCPList, nReqOrder, bReversed, FALSE, -1, -1);
     246}
     247
     248void *GDALCreateGCPRefineTransformer( int nGCPCount, const GDAL_GCP *pasGCPList,
     249                                int nReqOrder, int bReversed, double dfTolerance, int nMinimumGcps)
     250
     251{
     252    //If no minimumGcp parameter was passed, we  use the default value according to the model
     253    if(nMinimumGcps == -1)
     254    {
     255        nMinimumGcps = ((nReqOrder+1) * (nReqOrder+2)) / 2 + 1;
     256    }
     257    return GDALCreateGCPTransformerEx(nGCPCount, pasGCPList, nReqOrder, bReversed, TRUE, dfTolerance, nMinimumGcps);
     258}
     259
     260
    226261/************************************************************************/
    227262/*                     GDALDestroyGCPTransformer()                      */
    228263/************************************************************************/
     
    337372    CPLCreateXMLElementAndValue(
    338373        psTree, "Reversed",
    339374        CPLSPrintf( "%d", psInfo->bReversed ) );
     375
     376    CPLCreateXMLElementAndValue(
     377        psTree, "Refine",
     378        CPLSPrintf( "%d", psInfo->bRefine ) );
     379
     380    CPLCreateXMLElementAndValue(
     381        psTree, "MinimumGcps",
     382        CPLSPrintf( "%d", psInfo->nMinimumGcps ) );
     383
     384    CPLCreateXMLElementAndValue(
     385        psTree, "Tolerance",
     386        CPLSPrintf( "%f", psInfo->dfTolerance ) );
    340387                                 
    341388/* -------------------------------------------------------------------- */
    342389/*      Attach GCP List.                                                */
     
    347394        CPLXMLNode *psGCPList = CPLCreateXMLNode( psTree, CXT_Element,
    348395                                                  "GCPList" );
    349396
     397        if(psInfo->bRefine)
     398        {
     399          remove_outliers(psInfo);
     400        }
     401       
    350402        for( iGCP = 0; iGCP < psInfo->nGCPCount; iGCP++ )
    351403        {
    352404            CPLXMLNode *psXMLGCP;
     
    392444    void *pResult;
    393445    int nReqOrder;
    394446    int bReversed;
     447    int bRefine;
     448    int nMinimumGcps;
     449    double dfTolerance;
    395450
    396451    /* -------------------------------------------------------------------- */
    397452    /*      Check for GCPs.                                                 */
     
    442497/* -------------------------------------------------------------------- */
    443498    nReqOrder = atoi(CPLGetXMLValue(psTree,"Order","3"));
    444499    bReversed = atoi(CPLGetXMLValue(psTree,"Reversed","0"));
     500    bRefine = atoi(CPLGetXMLValue(psTree,"Refine","0"));
     501    nMinimumGcps = atoi(CPLGetXMLValue(psTree,"MinimumGcps","6"));
     502    dfTolerance = atof(CPLGetXMLValue(psTree,"Tolerance","1.0"));
    445503
    446504/* -------------------------------------------------------------------- */
    447505/*      Generate transformation.                                        */
    448506/* -------------------------------------------------------------------- */
    449     pResult = GDALCreateGCPTransformer( nGCPCount, pasGCPList, nReqOrder,
     507    if(bRefine)
     508    {
     509        pResult = GDALCreateGCPRefineTransformer( nGCPCount, pasGCPList, nReqOrder,
     510                                        bReversed, dfTolerance, nMinimumGcps );
     511    }
     512    else
     513    {
     514        pResult = GDALCreateGCPTransformer( nGCPCount, pasGCPList, nReqOrder,
    450515                                        bReversed );
     516    }
    451517   
    452518/* -------------------------------------------------------------------- */
    453519/*      Cleanup GCP copy.                                               */
     
    914980
    915981    return(MSUCCESS);
    916982}
     983
     984/***************************************************************************/
     985/*
     986  DETECTS THE WORST OUTLIER IN THE GCP LIST AND RETURNS THE INDEX OF THE
     987  OUTLIER.
     988 
     989  THE WORST OUTLIER IS CALCULATED BASED ON THE CONTROL POINTS, COEFFICIENTS
     990  AND THE ALLOWED TOLERANCE:
     991 
     992  sampleAdj = a0 + a1*sample + a2*line + a3*line*sample
     993  lineAdj = b0 + b1*sample + b2*line + b3*line*sample
     994 
     995  WHERE sampleAdj AND lineAdj ARE CORRELATED GCPS
     996 
     997  [residualSample] = [A1][sampleCoefficients] - [b1]
     998  [residualLine] = [A2][lineCoefficients] - [b2]
     999 
     1000  sampleResidual^2 = sum( [residualSample]^2 )
     1001  lineResidual^2 = sum( [lineSample]^2 )
     1002 
     1003  residuals(i) = squareRoot( residualSample(i)^2 + residualLine(i)^2 )
     1004 
     1005  THE GCP WITH THE GREATEST DISTANCE residual(i) GREATER THAN THE TOLERANCE WILL
     1006  CONSIDERED THE WORST OUTLIER.
     1007 
     1008  IF NO OUTLIER CAN BE FOUND, -1 WILL BE RETURNED.
     1009*/
     1010/***************************************************************************/
     1011static int worst_outlier(struct Control_Points *cp, double E[], double N[], double dfTolerance)
     1012{
     1013    double *padfResiduals;
     1014    int nI, nIndex;
     1015    double dfThreshold, dfDifference, dfSampleResidual, dfLineResidual, dfSampleRes, dfLineRes, dfCurrentDifference;
     1016    double dfE1, dfN1, dfE2, dfN2, dfEn;
     1017 
     1018    padfResiduals = (double *) CPLCalloc(sizeof(double),cp->count);
     1019    dfSampleResidual = 0.0;
     1020    dfLineResidual = 0.0;
     1021 
     1022    for(nI = 0; nI < cp->count; nI++)
     1023    {
     1024        dfE1 = cp->e1[nI];
     1025        dfN1 = cp->n1[nI];
     1026        dfE2 = dfE1 * dfE1;
     1027        dfN2 = dfN1 * dfN1;
     1028        dfEn = dfE1 * dfN1;
     1029
     1030        dfSampleRes = E[0] + E[1] * dfE1 + E[2] * dfN1 + E[3] * dfE2 + E[4] * dfEn + E[5] * dfN2 - cp->e2[nI];
     1031        dfLineRes = N[0] + N[1] * dfE1 + N[2] * dfN1 + N[3] * dfE2 + N[4] * dfEn + N[5] * dfN2 - cp->n2[nI];
     1032   
     1033        dfSampleResidual += dfSampleRes*dfSampleRes;
     1034        dfLineResidual += dfLineRes*dfLineRes;
     1035   
     1036        padfResiduals[nI] = sqrt(dfSampleRes*dfSampleRes + dfLineRes*dfLineRes);
     1037    }
     1038 
     1039    dfThreshold = dfTolerance * sqrt( (dfSampleResidual + dfLineResidual) / (double) cp->count );
     1040 
     1041    nIndex = -1;
     1042    dfDifference = -1.0;
     1043    for(nI = 0; nI < cp->count; nI++)
     1044    {
     1045        dfCurrentDifference = padfResiduals[nI];
     1046        if(fabs(dfCurrentDifference) < 1.19209290E-07F /*FLT_EPSILON*/)
     1047        {
     1048            dfCurrentDifference = 0.0;
     1049        }
     1050        if(dfCurrentDifference > dfDifference && dfCurrentDifference >= dfThreshold)
     1051        {
     1052            dfDifference = dfCurrentDifference;
     1053            nIndex = nI;
     1054        }
     1055    }
     1056    CPLFree( padfResiduals );
     1057    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/***************************************************************************/
     1074static 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    nGCPCount = psInfo->nGCPCount;
     1083    nMinimumGcps = psInfo->nMinimumGcps;
     1084    nReqOrder = psInfo->nOrder;
     1085    dfTolerance = psInfo->dfTolerance;
     1086   
     1087    padfGeoX = (double *) CPLCalloc(sizeof(double),nGCPCount);
     1088    padfGeoY = (double *) CPLCalloc(sizeof(double),nGCPCount);
     1089    padfRasterX = (double *) CPLCalloc(sizeof(double),nGCPCount);
     1090    padfRasterY = (double *) CPLCalloc(sizeof(double),nGCPCount);
     1091    panStatus = (int *) CPLCalloc(sizeof(int),nGCPCount);
     1092   
     1093    for( nI = 0; nI < nGCPCount; nI++ )
     1094    {
     1095        panStatus[nI] = 1;
     1096        padfGeoX[nI] = psInfo->pasGCPList[nI].dfGCPX;
     1097        padfGeoY[nI] = psInfo->pasGCPList[nI].dfGCPY;
     1098        padfRasterX[nI] = psInfo->pasGCPList[nI].dfGCPPixel;
     1099        padfRasterY[nI] = psInfo->pasGCPList[nI].dfGCPLine;
     1100    }
     1101
     1102    sPoints.count = nGCPCount;
     1103    sPoints.e1 = padfRasterX;
     1104    sPoints.n1 = padfRasterY;
     1105    sPoints.e2 = padfGeoX;
     1106    sPoints.n2 = padfGeoY;
     1107    sPoints.status = panStatus;
     1108 
     1109    nCRSresult = CRS_compute_georef_equations( &sPoints,
     1110                                      psInfo->adfToGeoX, psInfo->adfToGeoY,
     1111                                      psInfo->adfFromGeoX, psInfo->adfFromGeoY,
     1112                                      nReqOrder );
     1113
     1114    while(sPoints.count > nMinimumGcps)
     1115    {
     1116        int nGcpCounter, nIndex;
     1117        double *padfGeoXNew, *padfGeoYNew, *padfRasterXNew, *padfRasterYNew;
     1118        int *panStatusNew;
     1119   
     1120        nIndex = worst_outlier(&sPoints, psInfo->adfFromGeoX, psInfo->adfFromGeoY, dfTolerance);
     1121
     1122        //If no outliers were detected, stop the GCP elimination
     1123        if(nIndex == -1)
     1124        {
     1125            break;
     1126        }
     1127
     1128        padfGeoXNew = (double *) CPLCalloc(sizeof(double),sPoints.count - 1);
     1129        padfGeoYNew = (double *) CPLCalloc(sizeof(double),sPoints.count - 1);
     1130        padfRasterXNew = (double *) CPLCalloc(sizeof(double),sPoints.count - 1);
     1131        padfRasterYNew = (double *) CPLCalloc(sizeof(double),sPoints.count - 1);
     1132        panStatusNew = (int *) CPLCalloc(sizeof(int),sPoints.count - 1);
     1133
     1134        nGcpCounter = 0;
     1135        for(nI = 0; nI < sPoints.count; nI++)
     1136        {
     1137            if(nI != nIndex)
     1138            {
     1139                panStatusNew[nGcpCounter] = 1;
     1140                padfGeoXNew[nGcpCounter] = sPoints.e2[nI];
     1141                padfGeoYNew[nGcpCounter] = sPoints.n2[nI];
     1142                padfRasterXNew[nGcpCounter] = sPoints.e1[nI];
     1143                padfRasterYNew[nGcpCounter] = sPoints.n1[nI];
     1144                nGcpCounter++;
     1145            }
     1146        }
     1147
     1148        CPLFree( sPoints.e1 );
     1149        CPLFree( sPoints.n1 );
     1150        CPLFree( sPoints.e2 );
     1151        CPLFree( sPoints.n2 );
     1152        CPLFree( sPoints.status );
     1153
     1154        sPoints.count = sPoints.count - 1;
     1155        sPoints.e1 = padfRasterXNew;
     1156        sPoints.n1 = padfRasterYNew;
     1157        sPoints.e2 = padfGeoXNew;
     1158        sPoints.n2 = padfGeoYNew;
     1159        sPoints.status = panStatusNew;
     1160   
     1161        nCRSresult = CRS_compute_georef_equations( &sPoints,
     1162                                      psInfo->adfToGeoX, psInfo->adfToGeoY,
     1163                                      psInfo->adfFromGeoX, psInfo->adfFromGeoY,
     1164                                      nReqOrder );
     1165    }
     1166   
     1167    for( nI = 0; nI < sPoints.count; nI++ )
     1168    {
     1169        psInfo->pasGCPList[nI].dfGCPX = sPoints.e2[nI];
     1170        psInfo->pasGCPList[nI].dfGCPY = sPoints.n2[nI];
     1171        psInfo->pasGCPList[nI].dfGCPPixel = sPoints.e1[nI];
     1172        psInfo->pasGCPList[nI].dfGCPLine = sPoints.n1[nI];
     1173    }
     1174    psInfo->nGCPCount = sPoints.count;
     1175   
     1176    CPLFree( sPoints.e1 );
     1177    CPLFree( sPoints.n1 );
     1178    CPLFree( sPoints.e2 );
     1179    CPLFree( sPoints.n2 );
     1180    CPLFree( sPoints.status );
     1181    return nCRSresult;
     1182}