Index: apps/gdalwarp.cpp
===================================================================
--- apps/gdalwarp.cpp (revision 22652)
+++ apps/gdalwarp.cpp (working copy)
@@ -78,7 +78,7 @@
[-r resampling_method] [-wm memory_in_mb] [-multi] [-q]
[-cutline datasource] [-cl layer] [-cwhere expression]
[-csql statement] [-cblend dist_in_pixels] [-crop_to_cutline]
- [-of format] [-co "NAME=VALUE"]* [-overwrite]
+ [-of format] [-co "NAME=VALUE"]* [-overwrite] [-refine_gcps tolerance minimum_gcps]
srcfile* dstfile
\endverbatim
@@ -163,6 +163,12 @@
the output format driver. Multiple -co options may be listed. See
format specific documentation for legal creation options for each format.
+
-refine_gcps tolerance minimum_gcps: refines the GCPs by automatically eliminating outliers.
+Outliers will be eliminated until minimum_gcps are left or when no outliers can be detected.
+The tolerance is passed to adjust when a GCP will be eliminated.
+Not that GCP refinement only works with polynomial interpolation.
+The tolerance is in pixel units if no projection is available, otherwise it is in SRS units.
+If minimum_gcps is not provided, the minimum GCPs according to the polynomial model is used.
-cutline datasource:Enable use of a blend cutline from the name OGR support datasource.
-cl layername:Select the named layer from the
@@ -231,7 +237,7 @@
" [-r resampling_method] [-wm memory_in_mb] [-multi] [-q]\n"
" [-cutline datasource] [-cl layer] [-cwhere expression]\n"
" [-csql statement] [-cblend dist_in_pixels] [-crop_to_cutline]\n"
- " [-of format] [-co \"NAME=VALUE\"]* [-overwrite]\n"
+ " [-of format] [-co \"NAME=VALUE\"]* [-overwrite] [-refine_gcps tolerance minimum_gcps]\n"
" srcfile* dstfile\n"
"\n"
"Available resampling methods:\n"
@@ -388,6 +394,23 @@
{
papszTO = CSLSetNameValue( papszTO, "MAX_GCP_ORDER", argv[++i] );
}
+ else if( EQUAL(argv[i],"-refine_gcps") && i < argc-2 )
+ {
+ papszTO = CSLSetNameValue( papszTO, "REFINE_TOLERANCE", argv[++i] );
+ if(atof(argv[i]) < 0)
+ {
+ printf( "The tolerance for -refine_gcps may not be negative\n");
+ Usage();
+ }
+ if(atoi(argv[i+1]) >= 0 && isdigit(argv[i+1][0]))
+ {
+ papszTO = CSLSetNameValue( papszTO, "REFINE_MINIMUM_GCPS", argv[++i] );
+ }
+ else
+ {
+ papszTO = CSLSetNameValue( papszTO, "REFINE_MINIMUM_GCPS", "-1" );
+ }
+ }
else if( EQUAL(argv[i],"-tps") )
{
papszTO = CSLSetNameValue( papszTO, "METHOD", "GCP_TPS" );
@@ -774,7 +797,7 @@
CSLDestroy( papszCreateOptions );
papszCreateOptions = NULL;
}
-
+
if( hDstDS == NULL )
exit( 1 );
@@ -799,7 +822,7 @@
/* Check that there's at least one raster band */
/* -------------------------------------------------------------------- */
if ( GDALGetRasterCount(hSrcDS) == 0 )
- {
+ {
fprintf(stderr, "Input file %s has no raster bands.\n", papszSrcFiles[iSrc] );
exit( 1 );
}
@@ -811,7 +834,7 @@
/* Warns if the file has a color table and something more */
/* complicated than nearest neighbour resampling is asked */
/* -------------------------------------------------------------------- */
-
+
if ( eResampleAlg != GRA_NearestNeighbour &&
GDALGetRasterColorTable(GDALGetRasterBand(hSrcDS, 1)) != NULL)
{
Index: alg/gdaltransformer.cpp
===================================================================
--- alg/gdaltransformer.cpp (revision 22652)
+++ alg/gdaltransformer.cpp (working copy)
@@ -1000,6 +1000,8 @@
* SRC_SRS: WKT SRS to be used as an override for hSrcDS.
* DST_SRS: WKT SRS to be used as an override for hDstDS.
* GCPS_OK: If false, GCPs will not be used, default is TRUE.
+ * REFINE_MINIMUM_GCPS: The minimum amount of GCPs that should be available after the refinement.
+ * REFINE_TOLERANCE: The tolernace that specifies when a GCP will be eliminated.
* MAX_GCP_ORDER: the maximum order to use for GCP derived polynomials if
* possible. The default is to autoselect based on the number of GCPs.
* A value of -1 triggers use of Thin Plate Spline instead of polynomials.
@@ -1030,7 +1032,8 @@
GDALRPCInfo sRPCInfo;
const char *pszMethod = CSLFetchNameValue( papszOptions, "METHOD" );
const char *pszValue;
- int nOrder = 0, bGCPUseOK = TRUE;
+ int nOrder = 0, bGCPUseOK = TRUE, nMinimumGcps = -1, bRefine = FALSE;
+ double dfTolerance = 0.0;
const char *pszSrcWKT = CSLFetchNameValue( papszOptions, "SRC_SRS" );
const char *pszDstWKT = CSLFetchNameValue( papszOptions, "DST_SRS" );
@@ -1042,6 +1045,20 @@
if( pszValue )
bGCPUseOK = CSLTestBoolean(pszValue);
+ pszValue = CSLFetchNameValue( papszOptions, "REFINE_MINIMUM_GCPS" );
+ if( pszValue )
+ {
+ if( atoi(pszValue) != -1)
+ nMinimumGcps = atoi(pszValue);
+ }
+
+ pszValue = CSLFetchNameValue( papszOptions, "REFINE_TOLERANCE" );
+ if( pszValue )
+ {
+ dfTolerance = atof(pszValue);
+ bRefine = TRUE;
+ }
+
/* -------------------------------------------------------------------- */
/* Initialize the transform info. */
/* -------------------------------------------------------------------- */
@@ -1089,10 +1106,20 @@
&& (pszMethod == NULL || EQUAL(pszMethod,"GCP_POLYNOMIAL") )
&& GDALGetGCPCount( hSrcDS ) > 0 && nOrder >= 0 )
{
- psInfo->pSrcGCPTransformArg =
- GDALCreateGCPTransformer( GDALGetGCPCount( hSrcDS ),
- GDALGetGCPs( hSrcDS ), nOrder,
- FALSE );
+ if(bRefine)
+ {
+ psInfo->pSrcGCPTransformArg =
+ GDALCreateGCPRefineTransformer( GDALGetGCPCount( hSrcDS ),
+ GDALGetGCPs( hSrcDS ), nOrder,
+ FALSE, dfTolerance, nMinimumGcps );
+ }
+ else
+ {
+ psInfo->pSrcGCPTransformArg =
+ GDALCreateGCPTransformer( GDALGetGCPCount( hSrcDS ),
+ GDALGetGCPs( hSrcDS ), nOrder,
+ FALSE );
+ }
if( psInfo->pSrcGCPTransformArg == NULL )
{
Index: alg/gdal_alg.h
===================================================================
--- alg/gdal_alg.h (revision 22652)
+++ alg/gdal_alg.h (working copy)
@@ -165,6 +165,12 @@
void CPL_DLL *
GDALCreateGCPTransformer( int nGCPCount, const GDAL_GCP *pasGCPList,
int nReqOrder, int bReversed );
+
+/* GCP based transformer with refinement of the GCPs ... forward is to georef coordinates */
+void CPL_DLL *
+GDALCreateGCPRefineTransformer( int nGCPCount, const GDAL_GCP *pasGCPList,
+ int nReqOrder, int bReversed, double tolerance, int minimumGcps);
+
void CPL_DLL GDALDestroyGCPTransformer( void *pTransformArg );
int CPL_DLL GDALGCPTransform(
void *pTransformArg, int bDstToSrc, int nPointCount,
Index: alg/gdal_crs.c
===================================================================
--- alg/gdal_crs.c (revision 22652)
+++ alg/gdal_crs.c (working copy)
@@ -26,8 +26,10 @@
Added printout of trnfile. Triggered by BDEBUG.
Last Update: 1/27/92 Brian J. Buckley
Fixed bug so that only the active control points were used.
+ Last Update: 6/29/2011 C. F. Stallmann & R. van den Dool (South African National Space Agency)
+ GCP refinement added
+
-
Copyright (c) 1992, Michigan State University
Permission is hereby granted, free of charge, to any person obtaining a
@@ -55,6 +57,7 @@
#include "cpl_minixml.h"
#include "cpl_string.h"
+
CPL_CVSID("$Id$");
#define MAXORDER 3
@@ -69,6 +72,27 @@
int *status;
};
+typedef struct
+{
+ GDALTransformerInfo sTI;
+
+ double adfToGeoX[20];
+ double adfToGeoY[20];
+
+ double adfFromGeoX[20];
+ double adfFromGeoY[20];
+
+ int nOrder;
+ int bReversed;
+
+ int nGCPCount;
+ GDAL_GCP *pasGCPList;
+ int bRefine;
+ int nMinimumGcps;
+ double dfTolerance;
+
+} GCPTransformInfo;
+
CPL_C_START
CPLXMLNode *GDALSerializeGCPTransformer( void *pTransformArg );
void *GDALDeserializeGCPTransformer( CPLXMLNode *psTree );
@@ -79,8 +103,8 @@
double [], double [], int);
static int CRS_compute_georef_equations(struct Control_Points *,
double [], double [], double [], double [], int);
+static int remove_outliers(GCPTransformInfo *);
-
static char *CRS_error_message[] = {
"Failed to compute GCP transform: Not enough points available",
"Failed to compute GCP transform: Transform is not solvable",
@@ -89,25 +113,7 @@
"Failed to compute GCP transform: Internal error"
};
-typedef struct
-{
- GDALTransformerInfo sTI;
- double adfToGeoX[20];
- double adfToGeoY[20];
-
- double adfFromGeoX[20];
- double adfFromGeoY[20];
-
- int nOrder;
- int bReversed;
-
- int nGCPCount;
- GDAL_GCP *pasGCPList;
-
-} GCPTransformInfo;
-
-
/************************************************************************/
/* GDALCreateGCPTransformer() */
/************************************************************************/
@@ -139,8 +145,8 @@
* @return the transform argument or NULL if creation fails.
*/
-void *GDALCreateGCPTransformer( int nGCPCount, const GDAL_GCP *pasGCPList,
- int nReqOrder, int bReversed )
+void *GDALCreateGCPTransformerEx( int nGCPCount, const GDAL_GCP *pasGCPList,
+ int nReqOrder, int bReversed, int bRefine, double dfTolerance, int nMinimumGcps)
{
GCPTransformInfo *psInfo;
@@ -162,6 +168,9 @@
psInfo = (GCPTransformInfo *) CPLCalloc(sizeof(GCPTransformInfo),1);
psInfo->bReversed = bReversed;
psInfo->nOrder = nReqOrder;
+ psInfo->bRefine = bRefine;
+ psInfo->dfTolerance = dfTolerance;
+ psInfo->nMinimumGcps = nMinimumGcps;
psInfo->pasGCPList = GDALDuplicateGCPs( nGCPCount, pasGCPList );
psInfo->nGCPCount = nGCPCount;
@@ -171,45 +180,51 @@
psInfo->sTI.pfnTransform = GDALGCPTransform;
psInfo->sTI.pfnCleanup = GDALDestroyGCPTransformer;
psInfo->sTI.pfnSerialize = GDALSerializeGCPTransformer;
-
+
/* -------------------------------------------------------------------- */
-/* Allocate and initialize the working points list. */
+/* Compute the forward and reverse polynomials. */
/* -------------------------------------------------------------------- */
- padfGeoX = (double *) CPLCalloc(sizeof(double),nGCPCount);
- padfGeoY = (double *) CPLCalloc(sizeof(double),nGCPCount);
- padfRasterX = (double *) CPLCalloc(sizeof(double),nGCPCount);
- padfRasterY = (double *) CPLCalloc(sizeof(double),nGCPCount);
- panStatus = (int *) CPLCalloc(sizeof(int),nGCPCount);
-
- for( iGCP = 0; iGCP < nGCPCount; iGCP++ )
+
+ if(bRefine)
{
- panStatus[iGCP] = 1;
- padfGeoX[iGCP] = pasGCPList[iGCP].dfGCPX;
- padfGeoY[iGCP] = pasGCPList[iGCP].dfGCPY;
- padfRasterX[iGCP] = pasGCPList[iGCP].dfGCPPixel;
- padfRasterY[iGCP] = pasGCPList[iGCP].dfGCPLine;
+ nCRSresult = remove_outliers(psInfo);
}
-
- sPoints.count = nGCPCount;
- sPoints.e1 = padfRasterX;
- sPoints.n1 = padfRasterY;
- sPoints.e2 = padfGeoX;
- sPoints.n2 = padfGeoY;
- sPoints.status = panStatus;
+ else
+ {
+ /* -------------------------------------------------------------------- */
+ /* Allocate and initialize the working points list. */
+ /* -------------------------------------------------------------------- */
+ padfGeoX = (double *) CPLCalloc(sizeof(double),nGCPCount);
+ padfGeoY = (double *) CPLCalloc(sizeof(double),nGCPCount);
+ padfRasterX = (double *) CPLCalloc(sizeof(double),nGCPCount);
+ padfRasterY = (double *) CPLCalloc(sizeof(double),nGCPCount);
+ panStatus = (int *) CPLCalloc(sizeof(int),nGCPCount);
-/* -------------------------------------------------------------------- */
-/* Compute the forward and reverse polynomials. */
-/* -------------------------------------------------------------------- */
- nCRSresult = CRS_compute_georef_equations( &sPoints,
- psInfo->adfToGeoX, psInfo->adfToGeoY,
- psInfo->adfFromGeoX, psInfo->adfFromGeoY,
- nReqOrder );
+ for( iGCP = 0; iGCP < nGCPCount; iGCP++ )
+ {
+ panStatus[iGCP] = 1;
+ padfGeoX[iGCP] = pasGCPList[iGCP].dfGCPX;
+ padfGeoY[iGCP] = pasGCPList[iGCP].dfGCPY;
+ padfRasterX[iGCP] = pasGCPList[iGCP].dfGCPPixel;
+ padfRasterY[iGCP] = pasGCPList[iGCP].dfGCPLine;
+ }
- CPLFree( padfGeoX );
- CPLFree( padfGeoY );
- CPLFree( padfRasterX );
- CPLFree( padfRasterY );
- CPLFree( panStatus );
+ sPoints.count = nGCPCount;
+ sPoints.e1 = padfRasterX;
+ sPoints.n1 = padfRasterY;
+ sPoints.e2 = padfGeoX;
+ sPoints.n2 = padfGeoY;
+ sPoints.status = panStatus;
+ nCRSresult = CRS_compute_georef_equations( &sPoints,
+ psInfo->adfToGeoX, psInfo->adfToGeoY,
+ psInfo->adfFromGeoX, psInfo->adfFromGeoY,
+ nReqOrder );
+ CPLFree( padfGeoX );
+ CPLFree( padfGeoY );
+ CPLFree( padfRasterX );
+ CPLFree( padfRasterY );
+ CPLFree( panStatus );
+ }
if (nCRSresult != 1)
{
@@ -223,6 +238,26 @@
}
}
+void *GDALCreateGCPTransformer( int nGCPCount, const GDAL_GCP *pasGCPList,
+ int nReqOrder, int bReversed )
+
+{
+ return GDALCreateGCPTransformerEx(nGCPCount, pasGCPList, nReqOrder, bReversed, FALSE, -1, -1);
+}
+
+void *GDALCreateGCPRefineTransformer( int nGCPCount, const GDAL_GCP *pasGCPList,
+ int nReqOrder, int bReversed, double dfTolerance, int nMinimumGcps)
+
+{
+ //If no minimumGcp parameter was passed, we use the default value according to the model
+ if(nMinimumGcps == -1)
+ {
+ nMinimumGcps = ((nReqOrder+1) * (nReqOrder+2)) / 2 + 1;
+ }
+ return GDALCreateGCPTransformerEx(nGCPCount, pasGCPList, nReqOrder, bReversed, TRUE, dfTolerance, nMinimumGcps);
+}
+
+
/************************************************************************/
/* GDALDestroyGCPTransformer() */
/************************************************************************/
@@ -337,6 +372,18 @@
CPLCreateXMLElementAndValue(
psTree, "Reversed",
CPLSPrintf( "%d", psInfo->bReversed ) );
+
+ CPLCreateXMLElementAndValue(
+ psTree, "Refine",
+ CPLSPrintf( "%d", psInfo->bRefine ) );
+
+ CPLCreateXMLElementAndValue(
+ psTree, "MinimumGcps",
+ CPLSPrintf( "%d", psInfo->nMinimumGcps ) );
+
+ CPLCreateXMLElementAndValue(
+ psTree, "Tolerance",
+ CPLSPrintf( "%f", psInfo->dfTolerance ) );
/* -------------------------------------------------------------------- */
/* Attach GCP List. */
@@ -347,6 +394,11 @@
CPLXMLNode *psGCPList = CPLCreateXMLNode( psTree, CXT_Element,
"GCPList" );
+ if(psInfo->bRefine)
+ {
+ remove_outliers(psInfo);
+ }
+
for( iGCP = 0; iGCP < psInfo->nGCPCount; iGCP++ )
{
CPLXMLNode *psXMLGCP;
@@ -392,6 +444,9 @@
void *pResult;
int nReqOrder;
int bReversed;
+ int bRefine;
+ int nMinimumGcps;
+ double dfTolerance;
/* -------------------------------------------------------------------- */
/* Check for GCPs. */
@@ -442,12 +497,23 @@
/* -------------------------------------------------------------------- */
nReqOrder = atoi(CPLGetXMLValue(psTree,"Order","3"));
bReversed = atoi(CPLGetXMLValue(psTree,"Reversed","0"));
+ bRefine = atoi(CPLGetXMLValue(psTree,"Refine","0"));
+ nMinimumGcps = atoi(CPLGetXMLValue(psTree,"MinimumGcps","6"));
+ dfTolerance = atof(CPLGetXMLValue(psTree,"Tolerance","1.0"));
/* -------------------------------------------------------------------- */
/* Generate transformation. */
/* -------------------------------------------------------------------- */
- pResult = GDALCreateGCPTransformer( nGCPCount, pasGCPList, nReqOrder,
+ if(bRefine)
+ {
+ pResult = GDALCreateGCPRefineTransformer( nGCPCount, pasGCPList, nReqOrder,
+ bReversed, dfTolerance, nMinimumGcps );
+ }
+ else
+ {
+ pResult = GDALCreateGCPTransformer( nGCPCount, pasGCPList, nReqOrder,
bReversed );
+ }
/* -------------------------------------------------------------------- */
/* Cleanup GCP copy. */
@@ -914,3 +980,203 @@
return(MSUCCESS);
}
+
+/***************************************************************************/
+/*
+ DETECTS THE WORST OUTLIER IN THE GCP LIST AND RETURNS THE INDEX OF THE
+ OUTLIER.
+
+ THE WORST OUTLIER IS CALCULATED BASED ON THE CONTROL POINTS, COEFFICIENTS
+ AND THE ALLOWED TOLERANCE:
+
+ sampleAdj = a0 + a1*sample + a2*line + a3*line*sample
+ lineAdj = b0 + b1*sample + b2*line + b3*line*sample
+
+ WHERE sampleAdj AND lineAdj ARE CORRELATED GCPS
+
+ [residualSample] = [A1][sampleCoefficients] - [b1]
+ [residualLine] = [A2][lineCoefficients] - [b2]
+
+ sampleResidual^2 = sum( [residualSample]^2 )
+ lineResidual^2 = sum( [lineSample]^2 )
+
+ residuals(i) = squareRoot( residualSample(i)^2 + residualLine(i)^2 )
+
+ THE GCP WITH THE GREATEST DISTANCE residual(i) GREATER THAN THE TOLERANCE WILL
+ CONSIDERED THE WORST OUTLIER.
+
+ IF NO OUTLIER CAN BE FOUND, -1 WILL BE RETURNED.
+*/
+/***************************************************************************/
+static int worst_outlier(struct Control_Points *cp, double E[], double N[], double dfTolerance)
+{
+ double *padfResiduals;
+ int nI, nIndex;
+ double dfThreshold, dfDifference, dfSampleResidual, dfLineResidual, dfSampleRes, dfLineRes, dfCurrentDifference;
+ double dfE1, dfN1, dfE2, dfN2, dfEn;
+
+ padfResiduals = (double *) CPLCalloc(sizeof(double),cp->count);
+ dfSampleResidual = 0.0;
+ dfLineResidual = 0.0;
+
+ for(nI = 0; nI < cp->count; nI++)
+ {
+ dfE1 = cp->e1[nI];
+ dfN1 = cp->n1[nI];
+ dfE2 = dfE1 * dfE1;
+ dfN2 = dfN1 * dfN1;
+ dfEn = dfE1 * dfN1;
+
+ dfSampleRes = E[0] + E[1] * dfE1 + E[2] * dfN1 + E[3] * dfE2 + E[4] * dfEn + E[5] * dfN2 - cp->e2[nI];
+ dfLineRes = N[0] + N[1] * dfE1 + N[2] * dfN1 + N[3] * dfE2 + N[4] * dfEn + N[5] * dfN2 - cp->n2[nI];
+
+ dfSampleResidual += dfSampleRes*dfSampleRes;
+ dfLineResidual += dfLineRes*dfLineRes;
+
+ padfResiduals[nI] = sqrt(dfSampleRes*dfSampleRes + dfLineRes*dfLineRes);
+ }
+
+ dfThreshold = dfTolerance * sqrt( (dfSampleResidual + dfLineResidual) / (double) cp->count );
+
+ nIndex = -1;
+ dfDifference = -1.0;
+ for(nI = 0; nI < cp->count; nI++)
+ {
+ dfCurrentDifference = padfResiduals[nI];
+ if(fabs(dfCurrentDifference) < 1.19209290E-07F /*FLT_EPSILON*/)
+ {
+ dfCurrentDifference = 0.0;
+ }
+ if(dfCurrentDifference > dfDifference && dfCurrentDifference >= dfThreshold)
+ {
+ dfDifference = dfCurrentDifference;
+ nIndex = nI;
+ }
+ }
+ CPLFree( padfResiduals );
+ return nIndex;
+}
+
+/***************************************************************************/
+/*
+ REMOVES THE WORST OUTLIERS ITERATIVELY UNTIL THE MINIMUM NUMBER OF GCPS
+ ARE REACHED OR NO OUTLIERS CAN BE DETECTED.
+
+ 1. WE CALCULATE THE COEFFICIENTS FOR ALL THE GCPS.
+ 2. THE GCP LIST WILL BE SCANED TO DETERMINE THE WORST OUTLIER USING
+ THE CALCULATED COEFFICIENTS.
+ 3. THE WORST OUTLIER WILL BE REMOVED FROM THE GCP LIST.
+ 4. THE COEFFICIENTS WILL BE RECALCULATED WITHOUT THE WORST OUTLIER.
+ 5. STEP 1 TO 4 ARE EXECUTED UNTIL THE MINIMUM NUMBER OF GCPS WERE REACHED
+ OR IF NO GCP IS CONSIDERED AN OUTLIER WITH THE PASSED TOLERANCE.
+*/
+/***************************************************************************/
+static int remove_outliers( GCPTransformInfo *psInfo )
+{
+ double *padfGeoX, *padfGeoY, *padfRasterX, *padfRasterY;
+ int *panStatus;
+ int nI, nCRSresult, nGCPCount, nMinimumGcps, nReqOrder;
+ double dfTolerance;
+ struct Control_Points sPoints;
+
+ nGCPCount = psInfo->nGCPCount;
+ nMinimumGcps = psInfo->nMinimumGcps;
+ nReqOrder = psInfo->nOrder;
+ dfTolerance = psInfo->dfTolerance;
+
+ padfGeoX = (double *) CPLCalloc(sizeof(double),nGCPCount);
+ padfGeoY = (double *) CPLCalloc(sizeof(double),nGCPCount);
+ padfRasterX = (double *) CPLCalloc(sizeof(double),nGCPCount);
+ padfRasterY = (double *) CPLCalloc(sizeof(double),nGCPCount);
+ panStatus = (int *) CPLCalloc(sizeof(int),nGCPCount);
+
+ for( nI = 0; nI < nGCPCount; nI++ )
+ {
+ panStatus[nI] = 1;
+ padfGeoX[nI] = psInfo->pasGCPList[nI].dfGCPX;
+ padfGeoY[nI] = psInfo->pasGCPList[nI].dfGCPY;
+ padfRasterX[nI] = psInfo->pasGCPList[nI].dfGCPPixel;
+ padfRasterY[nI] = psInfo->pasGCPList[nI].dfGCPLine;
+ }
+
+ sPoints.count = nGCPCount;
+ sPoints.e1 = padfRasterX;
+ sPoints.n1 = padfRasterY;
+ sPoints.e2 = padfGeoX;
+ sPoints.n2 = padfGeoY;
+ sPoints.status = panStatus;
+
+ nCRSresult = CRS_compute_georef_equations( &sPoints,
+ psInfo->adfToGeoX, psInfo->adfToGeoY,
+ psInfo->adfFromGeoX, psInfo->adfFromGeoY,
+ nReqOrder );
+
+ while(sPoints.count > nMinimumGcps)
+ {
+ int nGcpCounter, nIndex;
+ double *padfGeoXNew, *padfGeoYNew, *padfRasterXNew, *padfRasterYNew;
+ int *panStatusNew;
+
+ nIndex = worst_outlier(&sPoints, psInfo->adfFromGeoX, psInfo->adfFromGeoY, dfTolerance);
+
+ //If no outliers were detected, stop the GCP elimination
+ if(nIndex == -1)
+ {
+ break;
+ }
+
+ padfGeoXNew = (double *) CPLCalloc(sizeof(double),sPoints.count - 1);
+ padfGeoYNew = (double *) CPLCalloc(sizeof(double),sPoints.count - 1);
+ padfRasterXNew = (double *) CPLCalloc(sizeof(double),sPoints.count - 1);
+ padfRasterYNew = (double *) CPLCalloc(sizeof(double),sPoints.count - 1);
+ panStatusNew = (int *) CPLCalloc(sizeof(int),sPoints.count - 1);
+
+ nGcpCounter = 0;
+ for(nI = 0; nI < sPoints.count; nI++)
+ {
+ if(nI != nIndex)
+ {
+ panStatusNew[nGcpCounter] = 1;
+ padfGeoXNew[nGcpCounter] = sPoints.e2[nI];
+ padfGeoYNew[nGcpCounter] = sPoints.n2[nI];
+ padfRasterXNew[nGcpCounter] = sPoints.e1[nI];
+ padfRasterYNew[nGcpCounter] = sPoints.n1[nI];
+ nGcpCounter++;
+ }
+ }
+
+ CPLFree( sPoints.e1 );
+ CPLFree( sPoints.n1 );
+ CPLFree( sPoints.e2 );
+ CPLFree( sPoints.n2 );
+ CPLFree( sPoints.status );
+
+ sPoints.count = sPoints.count - 1;
+ sPoints.e1 = padfRasterXNew;
+ sPoints.n1 = padfRasterYNew;
+ sPoints.e2 = padfGeoXNew;
+ sPoints.n2 = padfGeoYNew;
+ sPoints.status = panStatusNew;
+
+ nCRSresult = CRS_compute_georef_equations( &sPoints,
+ psInfo->adfToGeoX, psInfo->adfToGeoY,
+ psInfo->adfFromGeoX, psInfo->adfFromGeoY,
+ nReqOrder );
+ }
+
+ for( nI = 0; nI < sPoints.count; nI++ )
+ {
+ psInfo->pasGCPList[nI].dfGCPX = sPoints.e2[nI];
+ psInfo->pasGCPList[nI].dfGCPY = sPoints.n2[nI];
+ psInfo->pasGCPList[nI].dfGCPPixel = sPoints.e1[nI];
+ psInfo->pasGCPList[nI].dfGCPLine = sPoints.n1[nI];
+ }
+ psInfo->nGCPCount = sPoints.count;
+
+ CPLFree( sPoints.e1 );
+ CPLFree( sPoints.n1 );
+ CPLFree( sPoints.e2 );
+ CPLFree( sPoints.n2 );
+ CPLFree( sPoints.status );
+ return nCRSresult;
+}