Index: apps/gdalwarp.cpp
===================================================================
--- apps/gdalwarp.cpp (revision 22612)
+++ apps/gdalwarp.cpp (working copy)
@@ -77,7 +77,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]
srcfile* dstfile
\endverbatim
@@ -162,6 +162,9 @@
the output format driver. Multiple -co options may be listed. See
format specific documentation for legal creation options for each format.
+
-refine_gcps minimumGcps tolerance: refines the GCPs by automatically eliminating outliers.
+Outliers will be eliminated until minimumGcps are left or when no outliers can be detected.
+The tolerance is passed to adjust when a GCP will be eliminated.
-cutline datasource:Enable use of a blend cutline from the name OGR support datasource.
-cl layername:Select the named layer from the
@@ -230,7 +233,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]\n"
" srcfile* dstfile\n"
"\n"
"Available resampling methods:\n"
@@ -387,6 +390,11 @@
{
papszTO = CSLSetNameValue( papszTO, "MAX_GCP_ORDER", argv[++i] );
}
+ else if( EQUAL(argv[i],"-refine_gcps") && i < argc-2 )
+ {
+ papszTO = CSLSetNameValue( papszTO, "REFINE_MINIMUM", argv[++i] );
+ papszTO = CSLSetNameValue( papszTO, "REFINE_TOLERANCE", argv[++i] );
+ }
else if( EQUAL(argv[i],"-tps") )
{
papszTO = CSLSetNameValue( papszTO, "METHOD", "GCP_TPS" );
@@ -770,7 +778,7 @@
CSLDestroy( papszCreateOptions );
papszCreateOptions = NULL;
}
-
+
if( hDstDS == NULL )
exit( 1 );
@@ -795,7 +803,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 );
}
@@ -807,7 +815,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 22612)
+++ alg/gdaltransformer.cpp (working copy)
@@ -1030,7 +1030,8 @@
GDALRPCInfo sRPCInfo;
const char *pszMethod = CSLFetchNameValue( papszOptions, "METHOD" );
const char *pszValue;
- int nOrder = 0, bGCPUseOK = TRUE;
+ int nOrder = 0, bGCPUseOK = TRUE, minimumGcps = 0;
+ double tolerance = 0.0;
const char *pszSrcWKT = CSLFetchNameValue( papszOptions, "SRC_SRS" );
const char *pszDstWKT = CSLFetchNameValue( papszOptions, "DST_SRS" );
@@ -1042,6 +1043,14 @@
if( pszValue )
bGCPUseOK = CSLTestBoolean(pszValue);
+ pszValue = CSLFetchNameValue( papszOptions, "REFINE_MINIMUM" );
+ if( pszValue )
+ minimumGcps = atoi(pszValue);
+
+ pszValue = CSLFetchNameValue( papszOptions, "REFINE_TOLERANCE" );
+ if( pszValue )
+ tolerance = atof(pszValue);
+
/* -------------------------------------------------------------------- */
/* Initialize the transform info. */
/* -------------------------------------------------------------------- */
@@ -1089,10 +1098,20 @@
&& (pszMethod == NULL || EQUAL(pszMethod,"GCP_POLYNOMIAL") )
&& GDALGetGCPCount( hSrcDS ) > 0 && nOrder >= 0 )
{
- psInfo->pSrcGCPTransformArg =
- GDALCreateGCPTransformer( GDALGetGCPCount( hSrcDS ),
- GDALGetGCPs( hSrcDS ), nOrder,
- FALSE );
+ if(minimumGcps > 0)
+ {
+ psInfo->pSrcGCPTransformArg =
+ GDALCreateGCPRefineTransformer( GDALGetGCPCount( hSrcDS ),
+ GDALGetGCPs( hSrcDS ), nOrder,
+ FALSE, tolerance, minimumGcps );
+ }
+ else
+ {
+ psInfo->pSrcGCPTransformArg =
+ GDALCreateGCPTransformer( GDALGetGCPCount( hSrcDS ),
+ GDALGetGCPs( hSrcDS ), nOrder,
+ FALSE );
+ }
if( psInfo->pSrcGCPTransformArg == NULL )
{
Index: alg/gdal_alg.h
===================================================================
--- alg/gdal_alg.h (revision 22612)
+++ 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 22612)
+++ 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,24 @@
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;
+
+} GCPTransformInfo;
+
CPL_C_START
CPLXMLNode *GDALSerializeGCPTransformer( void *pTransformArg );
void *GDALDeserializeGCPTransformer( CPLXMLNode *psTree );
@@ -79,8 +100,8 @@
double [], double [], int);
static int CRS_compute_georef_equations(struct Control_Points *,
double [], double [], double [], double [], int);
+static int remove_outliers(struct Control_Points *, GCPTransformInfo *, int, double, int);
-
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 +110,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 +142,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 refine, double tolerance, int minimumGcps)
{
GCPTransformInfo *psInfo;
@@ -200,16 +203,23 @@
/* -------------------------------------------------------------------- */
/* Compute the forward and reverse polynomials. */
/* -------------------------------------------------------------------- */
- 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(refine)
+ {
+ nCRSresult = remove_outliers(&sPoints, psInfo, minimumGcps, tolerance, nReqOrder);
+ }
+ else
+ {
+ 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 +233,21 @@
}
}
+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 tolerance, int minimumGcps)
+
+{
+ return GDALCreateGCPTransformerEx(nGCPCount, pasGCPList, nReqOrder, bReversed, TRUE, tolerance, minimumGcps);
+}
+
+
/************************************************************************/
/* GDALDestroyGCPTransformer() */
/************************************************************************/
@@ -914,3 +939,123 @@
return(MSUCCESS);
}
+
+/***************************************************************************/
+/*
+ DETECTS THE WORST OUTLIER IN THE GCP LIST AND RETURNS THE INDEX OF THE
+ OUTLIER.
+*/
+/***************************************************************************/
+static int worst_outlier(struct Control_Points *cp, double E[], double N[], double tolerance)
+{
+ double *residuals;
+ int i, index;
+ double threshold, difference, sampleResidual, lineResidual, sampleRes, lineRes, currentDifference;
+ double e1, n1, e2, n2, en;
+
+ residuals = (double *) CPLCalloc(sizeof(double),cp->count);
+ sampleResidual = 0.0;
+ lineResidual = 0.0;
+
+ for(i = 0; i < cp->count; i++)
+ {
+ e1 = cp->e1[i];
+ n1 = cp->n1[i];
+ e2 = e1 * e1;
+ n2 = n1 * n1;
+ en = e1 * n1;
+
+ sampleRes = E[0] + E[1] * e1 + E[2] * n1 +
+ E[3] * e2 + E[4] * en + E[5] * n2 - cp->e2[i];
+ lineRes = N[0] + N[1] * e1 + N[2] * n1 +
+ N[3] * e2 + N[4] * en + N[5] * n2 - cp->n2[i];
+
+ sampleResidual += sampleRes*sampleRes;
+ lineResidual += lineRes*lineRes;
+
+ residuals[i] = sqrt(sampleRes*sampleRes + lineRes*lineRes);
+ }
+
+ threshold = tolerance * sqrt( (sampleResidual + lineResidual) / (double) cp->count );
+
+ index = -1;
+ difference = 0.0;
+ for(i = 0; i < cp->count; i++)
+ {
+ currentDifference = residuals[i] - threshold;
+ if(currentDifference > difference)
+ {
+ difference = currentDifference;
+ index = i;
+ }
+ }
+ return index;
+}
+
+/***************************************************************************/
+/*
+ REMOVES THE WORST OUTLIERS ITERATIVELY UNTIL THE MINIMUM NUMBER OF GCPS
+ ARE REACHED OR NO OUTLIERS CAN BE DETECTED.
+*/
+/***************************************************************************/
+static int remove_outliers(struct Control_Points *sPoints, GCPTransformInfo *psInfo, int minimumGcps, double tolerance, int nReqOrder)
+{
+ int nCRSresult;
+ nCRSresult = CRS_compute_georef_equations( sPoints,
+ psInfo->adfToGeoX, psInfo->adfToGeoY,
+ psInfo->adfFromGeoX, psInfo->adfFromGeoY,
+ nReqOrder );
+ while(sPoints->count > minimumGcps)
+ {
+ int gcpCounter, i, index;
+ double *padfGeoX, *padfGeoY, *padfRasterX, *padfRasterY;
+ int *panStatus;
+
+ index = worst_outlier(sPoints, psInfo->adfFromGeoX, psInfo->adfFromGeoY, tolerance);
+
+ //No outliers detected
+ if(index == -1)
+ {
+ break;
+ }
+
+ padfGeoX = (double *) CPLCalloc(sizeof(double),sPoints->count - 1);
+ padfGeoY = (double *) CPLCalloc(sizeof(double),sPoints->count - 1);
+ padfRasterX = (double *) CPLCalloc(sizeof(double),sPoints->count - 1);
+ padfRasterY = (double *) CPLCalloc(sizeof(double),sPoints->count - 1);
+ panStatus = (int *) CPLCalloc(sizeof(int),sPoints->count - 1);
+
+ gcpCounter = 0;
+ for(i = 0; i < sPoints->count; i++)
+ {
+ if(i != index)
+ {
+ panStatus[gcpCounter] = 1;
+ padfGeoX[gcpCounter] = sPoints->e2[i];
+ padfGeoY[gcpCounter] = sPoints->n2[i];
+ padfRasterX[gcpCounter] = sPoints->e1[i];
+ padfRasterY[gcpCounter] = sPoints->n1[i];
+ gcpCounter++;
+ }
+ }
+
+ CPLFree( sPoints->e1 );
+ CPLFree( sPoints->n1 );
+ CPLFree( sPoints->e2 );
+ CPLFree( sPoints->n2 );
+ CPLFree( sPoints->status );
+
+ sPoints->count = sPoints->count - 1;
+ 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 );
+ }
+ return nCRSresult;
+}