Ticket #4143: gcp_refinement.diff

File gcp_refinement.diff, 14.3 KB (added by goocreations, 5 years ago)

SVN patch for the GCP refinement

  • apps/gdalwarp.cpp

     
    7777    [-r resampling_method] [-wm memory_in_mb] [-multi] [-q]
    7878    [-cutline datasource] [-cl layer] [-cwhere expression]
    7979    [-csql statement] [-cblend dist_in_pixels] [-crop_to_cutline]
    80     [-of format] [-co "NAME=VALUE"]* [-overwrite]
     80    [-of format] [-co "NAME=VALUE"]* [-overwrite] [-refine_gcps]
    8181    srcfile* dstfile
    8282\endverbatim
    8383
     
    162162the output format driver. Multiple <b>-co</b> options may be listed. See
    163163format specific documentation for legal creation options for each format.
    164164</dd>
     165<dt> <b>-refine_gcps</b> <em>minimumGcps tolerance</em>:</dt><dd> refines the GCPs by automatically eliminating outliers.
     166Outliers will be eliminated until minimumGcps are left or when no outliers can be detected.
     167The tolerance is passed to adjust when a GCP will be eliminated.</dd>
    165168
    166169<dt> <b>-cutline</b> <em>datasource</em>:</dt><dd>Enable use of a blend cutline from the name OGR support datasource.</dd>
    167170<dt> <b>-cl</b> <em>layername</em>:</dt><dd>Select the named layer from the
     
    230233        "    [-r resampling_method] [-wm memory_in_mb] [-multi] [-q]\n"
    231234        "    [-cutline datasource] [-cl layer] [-cwhere expression]\n"
    232235        "    [-csql statement] [-cblend dist_in_pixels] [-crop_to_cutline]\n"
    233         "    [-of format] [-co \"NAME=VALUE\"]* [-overwrite]\n"
     236        "    [-of format] [-co \"NAME=VALUE\"]* [-overwrite] [-refine_gcps]\n"
    234237        "    srcfile* dstfile\n"
    235238        "\n"
    236239        "Available resampling methods:\n"
     
    387390        {
    388391            papszTO = CSLSetNameValue( papszTO, "MAX_GCP_ORDER", argv[++i] );
    389392        }
     393        else if( EQUAL(argv[i],"-refine_gcps") && i < argc-2 )
     394        {
     395            papszTO = CSLSetNameValue( papszTO, "REFINE_MINIMUM", argv[++i] );
     396            papszTO = CSLSetNameValue( papszTO, "REFINE_TOLERANCE", argv[++i] );
     397        }
    390398        else if( EQUAL(argv[i],"-tps") )
    391399        {
    392400            papszTO = CSLSetNameValue( papszTO, "METHOD", "GCP_TPS" );
     
    770778        CSLDestroy( papszCreateOptions );
    771779        papszCreateOptions = NULL;
    772780    }
    773 
     781 
    774782    if( hDstDS == NULL )
    775783        exit( 1 );
    776784
     
    795803/*      Check that there's at least one raster band                     */
    796804/* -------------------------------------------------------------------- */
    797805        if ( GDALGetRasterCount(hSrcDS) == 0 )
    798         {
     806        {     
    799807            fprintf(stderr, "Input file %s has no raster bands.\n", papszSrcFiles[iSrc] );
    800808            exit( 1 );
    801809        }
     
    807815/*      Warns if the file has a color table and something more          */
    808816/*      complicated than nearest neighbour resampling is asked          */
    809817/* -------------------------------------------------------------------- */
    810 
     818 
    811819        if ( eResampleAlg != GRA_NearestNeighbour &&
    812820             GDALGetRasterColorTable(GDALGetRasterBand(hSrcDS, 1)) != NULL)
    813821        {
  • alg/gdaltransformer.cpp

     
    10301030    GDALRPCInfo sRPCInfo;
    10311031    const char *pszMethod = CSLFetchNameValue( papszOptions, "METHOD" );
    10321032    const char *pszValue;
    1033     int nOrder = 0, bGCPUseOK = TRUE;
     1033    int nOrder = 0, bGCPUseOK = TRUE, minimumGcps = 0;
     1034    double tolerance = 0.0;
    10341035    const char *pszSrcWKT = CSLFetchNameValue( papszOptions, "SRC_SRS" );
    10351036    const char *pszDstWKT = CSLFetchNameValue( papszOptions, "DST_SRS" );
    10361037
     
    10421043    if( pszValue )
    10431044        bGCPUseOK = CSLTestBoolean(pszValue);
    10441045
     1046    pszValue = CSLFetchNameValue( papszOptions, "REFINE_MINIMUM" );
     1047    if( pszValue )
     1048        minimumGcps = atoi(pszValue);
     1049   
     1050    pszValue = CSLFetchNameValue( papszOptions, "REFINE_TOLERANCE" );
     1051    if( pszValue )
     1052        tolerance = atof(pszValue);
     1053
    10451054/* -------------------------------------------------------------------- */
    10461055/*      Initialize the transform info.                                  */
    10471056/* -------------------------------------------------------------------- */
     
    10891098             && (pszMethod == NULL || EQUAL(pszMethod,"GCP_POLYNOMIAL") )
    10901099             && GDALGetGCPCount( hSrcDS ) > 0 && nOrder >= 0 )
    10911100    {
    1092         psInfo->pSrcGCPTransformArg =
    1093             GDALCreateGCPTransformer( GDALGetGCPCount( hSrcDS ),
    1094                                       GDALGetGCPs( hSrcDS ), nOrder,
    1095                                       FALSE );
     1101        if(minimumGcps > 0)
     1102        {
     1103          psInfo->pSrcGCPTransformArg =
     1104                GDALCreateGCPRefineTransformer( GDALGetGCPCount( hSrcDS ),
     1105                                          GDALGetGCPs( hSrcDS ), nOrder,
     1106                                          FALSE, tolerance, minimumGcps );
     1107        }
     1108        else
     1109        {
     1110          psInfo->pSrcGCPTransformArg =
     1111                GDALCreateGCPTransformer( GDALGetGCPCount( hSrcDS ),
     1112                                          GDALGetGCPs( hSrcDS ), nOrder,
     1113                                          FALSE );
     1114        }
    10961115
    10971116        if( psInfo->pSrcGCPTransformArg == NULL )
    10981117        {
  • 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   
     91} GCPTransformInfo;
     92
    7293CPL_C_START
    7394CPLXMLNode *GDALSerializeGCPTransformer( void *pTransformArg );
    7495void *GDALDeserializeGCPTransformer( CPLXMLNode *psTree );
     
    79100                              double [], double [], int);
    80101static int CRS_compute_georef_equations(struct Control_Points *,
    81102    double [], double [], double [], double [], int);
     103static int remove_outliers(struct Control_Points *, GCPTransformInfo *, int, double, int);
    82104
    83 
    84105static char *CRS_error_message[] = {
    85106    "Failed to compute GCP transform: Not enough points available",
    86107    "Failed to compute GCP transform: Transform is not solvable",
     
    89110    "Failed to compute GCP transform: Internal error"
    90111};
    91112
    92 typedef struct
    93 {
    94     GDALTransformerInfo sTI;
    95113
    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 
    111114/************************************************************************/
    112115/*                      GDALCreateGCPTransformer()                      */
    113116/************************************************************************/
     
    139142 * @return the transform argument or NULL if creation fails.
    140143 */
    141144
    142 void *GDALCreateGCPTransformer( int nGCPCount, const GDAL_GCP *pasGCPList,
    143                                 int nReqOrder, int bReversed )
     145void *GDALCreateGCPTransformerEx( int nGCPCount, const GDAL_GCP *pasGCPList,
     146                                int nReqOrder, int bReversed, int refine, double tolerance, int minimumGcps)
    144147
    145148{
    146149    GCPTransformInfo *psInfo;
     
    200203/* -------------------------------------------------------------------- */
    201204/*      Compute the forward and reverse polynomials.                    */
    202205/* -------------------------------------------------------------------- */
    203     nCRSresult = CRS_compute_georef_equations( &sPoints,
    204                                       psInfo->adfToGeoX, psInfo->adfToGeoY,
    205                                       psInfo->adfFromGeoX, psInfo->adfFromGeoY,
    206                                       nReqOrder );
    207206
    208     CPLFree( padfGeoX );
    209     CPLFree( padfGeoY );
    210     CPLFree( padfRasterX );
    211     CPLFree( padfRasterY );
    212     CPLFree( panStatus );
     207    if(refine)
     208    {
     209      nCRSresult = remove_outliers(&sPoints, psInfo, minimumGcps, tolerance, nReqOrder);
     210    }
     211    else
     212    {
     213      nCRSresult = CRS_compute_georef_equations( &sPoints,
     214                                            psInfo->adfToGeoX, psInfo->adfToGeoY,
     215                                            psInfo->adfFromGeoX, psInfo->adfFromGeoY,
     216                                            nReqOrder );
     217      CPLFree( padfGeoX );
     218      CPLFree( padfGeoY );
     219      CPLFree( padfRasterX );
     220      CPLFree( padfRasterY );
     221      CPLFree( panStatus );
     222    }
    213223
    214224    if (nCRSresult != 1)
    215225    {
     
    223233    }
    224234}
    225235
     236void *GDALCreateGCPTransformer( int nGCPCount, const GDAL_GCP *pasGCPList,
     237                                int nReqOrder, int bReversed )
     238
     239{
     240  return GDALCreateGCPTransformerEx(nGCPCount, pasGCPList, nReqOrder, bReversed, FALSE, -1, -1);
     241}
     242
     243void *GDALCreateGCPRefineTransformer( int nGCPCount, const GDAL_GCP *pasGCPList,
     244                                int nReqOrder, int bReversed, double tolerance, int minimumGcps)
     245
     246{
     247  return GDALCreateGCPTransformerEx(nGCPCount, pasGCPList, nReqOrder, bReversed, TRUE, tolerance, minimumGcps);
     248}
     249
     250
    226251/************************************************************************/
    227252/*                     GDALDestroyGCPTransformer()                      */
    228253/************************************************************************/
     
    914939
    915940    return(MSUCCESS);
    916941}
     942
     943/***************************************************************************/
     944/*
     945  DETECTS THE WORST OUTLIER IN THE GCP LIST AND RETURNS THE INDEX OF THE
     946  OUTLIER.
     947*/
     948/***************************************************************************/
     949static int worst_outlier(struct Control_Points *cp, double E[], double N[], double tolerance)
     950{
     951  double *residuals;
     952  int i, index;
     953  double threshold, difference, sampleResidual, lineResidual, sampleRes, lineRes, currentDifference;
     954  double e1, n1, e2, n2, en;
     955 
     956  residuals = (double *) CPLCalloc(sizeof(double),cp->count);
     957  sampleResidual = 0.0;
     958  lineResidual = 0.0;
     959 
     960  for(i = 0; i < cp->count; i++)
     961  {
     962    e1 = cp->e1[i];
     963    n1 = cp->n1[i];
     964    e2 = e1 * e1;
     965    n2 = n1 * n1;
     966    en = e1 * n1;
     967
     968    sampleRes = E[0] + E[1] * e1 + E[2] * n1 +
     969                   E[3] * e2 + E[4] * en + E[5] * n2 - cp->e2[i];
     970    lineRes = N[0] + N[1] * e1 + N[2] * n1 +
     971                 N[3] * e2 + N[4] * en + N[5] * n2 - cp->n2[i];
     972   
     973    sampleResidual += sampleRes*sampleRes;
     974    lineResidual += lineRes*lineRes;
     975   
     976    residuals[i] = sqrt(sampleRes*sampleRes + lineRes*lineRes);
     977  }
     978 
     979  threshold = tolerance * sqrt( (sampleResidual + lineResidual) / (double) cp->count );
     980 
     981  index = -1;
     982  difference = 0.0;
     983  for(i = 0; i < cp->count; i++)
     984  {
     985    currentDifference = residuals[i] - threshold;
     986    if(currentDifference > difference)
     987    {
     988      difference = currentDifference;
     989      index = i;
     990    }
     991  }
     992  return index;
     993}
     994
     995/***************************************************************************/
     996/*
     997  REMOVES THE WORST OUTLIERS ITERATIVELY UNTIL THE MINIMUM NUMBER OF GCPS
     998  ARE REACHED OR NO OUTLIERS CAN BE DETECTED.
     999*/
     1000/***************************************************************************/
     1001static int remove_outliers(struct Control_Points *sPoints, GCPTransformInfo *psInfo, int minimumGcps, double tolerance, int nReqOrder)
     1002{
     1003  int nCRSresult;
     1004  nCRSresult = CRS_compute_georef_equations( sPoints,
     1005                                      psInfo->adfToGeoX, psInfo->adfToGeoY,
     1006                                      psInfo->adfFromGeoX, psInfo->adfFromGeoY,
     1007                                      nReqOrder );
     1008  while(sPoints->count > minimumGcps)
     1009  {
     1010    int gcpCounter, i, index;
     1011    double *padfGeoX, *padfGeoY, *padfRasterX, *padfRasterY;
     1012    int *panStatus;
     1013   
     1014    index = worst_outlier(sPoints, psInfo->adfFromGeoX, psInfo->adfFromGeoY, tolerance);
     1015
     1016    //No outliers detected
     1017    if(index == -1)
     1018    {
     1019      break;
     1020    }
     1021
     1022    padfGeoX = (double *) CPLCalloc(sizeof(double),sPoints->count - 1);
     1023    padfGeoY = (double *) CPLCalloc(sizeof(double),sPoints->count - 1);
     1024    padfRasterX = (double *) CPLCalloc(sizeof(double),sPoints->count - 1);
     1025    padfRasterY = (double *) CPLCalloc(sizeof(double),sPoints->count - 1);
     1026    panStatus = (int *) CPLCalloc(sizeof(int),sPoints->count - 1);
     1027
     1028    gcpCounter = 0;
     1029    for(i = 0; i < sPoints->count; i++)
     1030    {
     1031      if(i != index)
     1032      {
     1033        panStatus[gcpCounter] = 1;
     1034        padfGeoX[gcpCounter] = sPoints->e2[i];
     1035        padfGeoY[gcpCounter] = sPoints->n2[i];
     1036        padfRasterX[gcpCounter] = sPoints->e1[i];
     1037        padfRasterY[gcpCounter] = sPoints->n1[i];
     1038        gcpCounter++;
     1039      }
     1040    }
     1041
     1042    CPLFree( sPoints->e1 );
     1043    CPLFree( sPoints->n1 );
     1044    CPLFree( sPoints->e2 );
     1045    CPLFree( sPoints->n2 );
     1046    CPLFree( sPoints->status );
     1047
     1048    sPoints->count = sPoints->count - 1;
     1049    sPoints->e1 = padfRasterX;
     1050    sPoints->n1 = padfRasterY;
     1051    sPoints->e2 = padfGeoX;
     1052    sPoints->n2 = padfGeoY;
     1053    sPoints->status = panStatus;
     1054   
     1055    nCRSresult = CRS_compute_georef_equations( sPoints,
     1056                                      psInfo->adfToGeoX, psInfo->adfToGeoY,
     1057                                      psInfo->adfFromGeoX, psInfo->adfFromGeoY,
     1058                                      nReqOrder );
     1059  }
     1060  return nCRSresult;
     1061}