Main Page Class Hierarchy Compound List File List Compound Members File Members Related Pages

# gdaltransformer.cpp

    00001 /******************************************************************************
    00002  * $Id: gdaltransformer_cpp-source.html,v 1.2 2002/12/21 19:13:13 warmerda Exp $
    00003  *
    00004  * Project:  Mapinfo Image Warper
    00005  * Purpose:  Implementation of one or more GDALTrasformerFunc types, including
    00006  *           the GenImgProj (general image reprojector) transformer.
    00007  * Author:   Frank Warmerdam, warmerdam@pobox.com
    00008  *
    00009  ******************************************************************************
    00010  * Copyright (c) 2002, i3 - information integration and imaging
    00011  *                          Fort Collin, CO
    00012  *
    00013  * Permission is hereby granted, free of charge, to any person obtaining a
    00014  * copy of this software and associated documentation files (the "Software"),
    00015  * to deal in the Software without restriction, including without limitation
    00016  * the rights to use, copy, modify, merge, publish, distribute, sublicense,
    00017  * and/or sell copies of the Software, and to permit persons to whom the
    00018  * Software is furnished to do so, subject to the following conditions:
    00019  *
    00020  * The above copyright notice and this permission notice shall be included
    00021  * in all copies or substantial portions of the Software.
    00022  *
    00023  * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
    00024  * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
    00025  * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
    00026  * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
    00027  * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
    00028  * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
    00029  * DEALINGS IN THE SOFTWARE.
    00030  ******************************************************************************
    00031  *
    00032  * $Log: gdaltransformer_cpp-source.html,v $
    00032  * Revision 1.2  2002/12/21 19:13:13  warmerda
    00032  * updated
    00032  *
    00033  * Revision 1.8  2002/12/13 15:55:27  warmerda
    00034  * fix suggested output to preserve orig size exact on null transform
    00035  *
    00036  * Revision 1.7  2002/12/13 04:57:49  warmerda
    00037  * remove approx transformer debug output
    00038  *
    00039  * Revision 1.6  2002/12/09 16:08:32  warmerda
    00040  * added approximating transformer
    00041  *
    00042  * Revision 1.5  2002/12/07 17:09:38  warmerda
    00043  * added order flag to GenImgProjTransformer
    00044  *
    00045  * Revision 1.4  2002/12/06 21:43:56  warmerda
    00046  * added GCP support to general transformer
    00047  *
    00048  * Revision 1.3  2002/12/05 21:45:01  warmerda
    00049  * fixed a few GenImgProj bugs
    00050  *
    00051  * Revision 1.2  2002/12/05 16:37:46  warmerda
    00052  * fixed method name
    00053  *
    00054  * Revision 1.1  2002/12/05 05:43:23  warmerda
    00055  * New
    00056  *
    00057  */
    00058
    00059 #include "gdal_priv.h"
    00060 #include "gdal_alg.h"
    00061 #include "ogr_spatialref.h"
    00062
    00063 CPL_CVSID("$Id: gdaltransformer_cpp-source.html,v 1.2 2002/12/21 19:13:13 warmerda Exp $");
    00064
    00065 /************************************************************************/
    00066 /*                          InvGeoTransform()                           */
    00067 /*                                                                      */
    00068 /*      Invert a standard 3x2 "GeoTransform" style matrix with an       */
    00069 /*      implicit [1 0 0] final row.                                     */
    00070 /************************************************************************/
    00071
    00072 static int InvGeoTransform( double *gt_in, double *gt_out )
    00073
    00074 {
    00075     double      det, inv_det;
    00076
    00077     /* we assume a 3rd row that is [1 0 0] */
    00078
    00079     /* Compute determinate */
    00080
    00081     det = gt_in[1] * gt_in[5] - gt_in[2] * gt_in[4];
    00082
    00083     if( fabs(det) < 0.000000000000001 )
    00084         return 0;
    00085
    00086     inv_det = 1.0 / det;
    00087
    00088     /* compute adjoint, and devide by determinate */
    00089
    00090     gt_out[1] =  gt_in[5] * inv_det;
    00091     gt_out[4] = -gt_in[4] * inv_det;
    00092
    00093     gt_out[2] = -gt_in[2] * inv_det;
    00094     gt_out[5] =  gt_in[1] * inv_det;
    00095
    00096     gt_out[0] = ( gt_in[2] * gt_in[3] - gt_in[0] * gt_in[5]) * inv_det;
    00097     gt_out[3] = (-gt_in[1] * gt_in[3] + gt_in[0] * gt_in[4]) * inv_det;
    00098
    00099     return 1;
    00100 }
    00101
    00102 /************************************************************************/
    00103 /*                      GDALSuggestedWarpOutput()                       */
    00104 /************************************************************************/
    00105
    00106 CPLErr GDALSuggestedWarpOutput( GDALDatasetH hSrcDS,
    00107                                 GDALTransformerFunc pfnTransformer,
    00108                                 void *pTransformArg,
    00109                                 double *padfGeoTransformOut,
    00110                                 int *pnPixels, int *pnLines )
    00111
    00112 {
    00113 /* -------------------------------------------------------------------- */
    00114 /*      Setup sample points all around the edge of the input raster.    */
    00115 /* -------------------------------------------------------------------- */
    00116     int    nSamplePoints=0, abSuccess[84];
    00117     double adfX[84], adfY[84], adfZ[84], dfRatio;
    00118     int    nInXSize = GDALGetRasterXSize( hSrcDS );
    00119     int    nInYSize = GDALGetRasterYSize( hSrcDS );
    00120
    00121     // Take 20 steps
    00122     for( dfRatio = 0.0; dfRatio <= 1.01; dfRatio += 0.05 )
    00123     {
    00124
    00125         // Ensure we end exactly at the end.
    00126         if( dfRatio > 0.99 )
    00127             dfRatio = 1.0;
    00128
    00129         // Along top
    00130         adfX[nSamplePoints]   = dfRatio * nInXSize;
    00131         adfY[nSamplePoints]   = 0.0;
    00132         adfZ[nSamplePoints++] = 0.0;
    00133
    00134         // Along bottom
    00135         adfX[nSamplePoints]   = dfRatio * nInXSize;
    00136         adfY[nSamplePoints]   = nInYSize;
    00137         adfZ[nSamplePoints++] = 0.0;
    00138
    00139         // Along left
    00140         adfX[nSamplePoints]   = 0.0;
    00141         adfY[nSamplePoints] = dfRatio * nInYSize;
    00142         adfZ[nSamplePoints++] = 0.0;
    00143
    00144         // Along right
    00145         adfX[nSamplePoints]   = nInXSize;
    00146         adfY[nSamplePoints] = dfRatio * nInYSize;
    00147         adfZ[nSamplePoints++] = 0.0;
    00148     }
    00149
    00150     CPLAssert( nSamplePoints == 84 );
    00151
    00152 /* -------------------------------------------------------------------- */
    00153 /*      Transform them to the output coordinate system.                 */
    00154 /* -------------------------------------------------------------------- */
    00155     if( !pfnTransformer( pTransformArg, FALSE, nSamplePoints,
    00156                          adfX, adfY, adfZ, abSuccess ) )
    00157     {
    00158         CPLError( CE_Failure, CPLE_AppDefined,
    00159                   "GDALSuggestedWarpOutput() failed because the passed\n"
    00160                   "transformer failed." );
    00161         return CE_Failure;
    00162     }
    00163
    00164 /* -------------------------------------------------------------------- */
    00165 /*      Collect the bounds, ignoring any failed points.                 */
    00166 /* -------------------------------------------------------------------- */
    00167     double dfMinXOut, dfMinYOut, dfMaxXOut, dfMaxYOut;
    00168     int    bGotInitialPoint = FALSE;
    00169     int    nFailedCount = 0, i;
    00170
    00171     for( i = 0; i < nSamplePoints; i++ )
    00172     {
    00173         if( !abSuccess[i] )
    00174         {
    00175             nFailedCount++;
    00176             continue;
    00177         }
    00178
    00179         if( !bGotInitialPoint )
    00180         {
    00181             bGotInitialPoint = TRUE;
    00182             dfMinXOut = dfMaxXOut = adfX[i];
    00183             dfMinYOut = dfMaxYOut = adfY[i];
    00184         }
    00185         else
    00186         {
    00187             dfMinXOut = MIN(dfMinXOut,adfX[i]);
    00188             dfMinYOut = MIN(dfMinYOut,adfY[i]);
    00189             dfMaxXOut = MAX(dfMaxXOut,adfX[i]);
    00190             dfMaxYOut = MAX(dfMaxYOut,adfY[i]);
    00191         }
    00192     }
    00193
    00194     if( nFailedCount > nSamplePoints - 10 )
    00195     {
    00196         CPLError( CE_Failure, CPLE_AppDefined,
    00197                   "Too many points (%d out of %d) failed to transform,\n"
    00198                   "unable to compute output bounds.",
    00199                   nFailedCount, nSamplePoints );
    00200         return CE_Failure;
    00201     }
    00202
    00203     if( nFailedCount > 0 )
    00204         CPLDebug( "GDAL",
    00205                   "GDALSuggestedWarpOutput(): %d out of %d points failed to transform.",
    00206                   nFailedCount, nSamplePoints );
    00207
    00208 /* -------------------------------------------------------------------- */
    00209 /*      Compute the distance in "georeferenced" units from the top      */
    00210 /*      corner of the transformed input image to the bottom left        */
    00211 /*      corner of the transformed input.  Use this distance to          */
    00212 /*      compute an approximate pixel size in the output                 */
    00213 /*      georeferenced coordinates.                                      */
    00214 /* -------------------------------------------------------------------- */
    00215     double dfDiagonalDist, dfDeltaX, dfDeltaY;
    00216
    00217     dfDeltaX = adfX[nSamplePoints-1] - adfX[0];
    00218     dfDeltaY = adfY[nSamplePoints-1] - adfY[0];
    00219
    00220     dfDiagonalDist = sqrt( dfDeltaX * dfDeltaX + dfDeltaY * dfDeltaY );
    00221
    00222 /* -------------------------------------------------------------------- */
    00223 /*      Compute a pixel size from this.                                 */
    00224 /* -------------------------------------------------------------------- */
    00225     double dfPixelSize;
    00226
    00227     dfPixelSize = dfDiagonalDist
    00228         / sqrt(((double)nInXSize)*nInXSize + ((double)nInYSize)*nInYSize);
    00229
    00230     *pnPixels = (int) ((dfMaxXOut - dfMinXOut) / dfPixelSize + 0.5);
    00231     *pnLines = (int) ((dfMaxYOut - dfMinYOut) / dfPixelSize + 0.5);
    00232
    00233 /* -------------------------------------------------------------------- */
    00234 /*      Set the output geotransform.                                    */
    00235 /* -------------------------------------------------------------------- */
    00236     padfGeoTransformOut[0] = dfMinXOut;
    00237     padfGeoTransformOut[1] = dfPixelSize;
    00238     padfGeoTransformOut[2] = 0.0;
    00239     padfGeoTransformOut[3] = dfMaxYOut;
    00240     padfGeoTransformOut[4] = 0.0;
    00241     padfGeoTransformOut[5] = - dfPixelSize;
    00242
    00243     return CE_None;
    00244 }
    00245
    00246 /************************************************************************/
    00247 /* ==================================================================== */
    00248 /*                       GDALGenImgProjTransformer                      */
    00249 /* ==================================================================== */
    00250 /************************************************************************/
    00251
    00252 typedef struct {
    00253
    00254     double   adfSrcGeoTransform[6];
    00255     double   adfSrcInvGeoTransform[6];
    00256
    00257     void     *pSrcGCPTransformArg;
    00258
    00259     void     *pReprojectArg;
    00260
    00261     double   adfDstGeoTransform[6];
    00262     double   adfDstInvGeoTransform[6];
    00263
    00264     void     *pDstGCPTransformArg;
    00265
    00266 } GDALGenImgProjTransformInfo;
    00267
    00268 /************************************************************************/
    00269 /*                  GDALCreateGenImgProjTransformer()                   */
    00270 /************************************************************************/
    00271
    00272 void *
    00273 GDALCreateGenImgProjTransformer( GDALDatasetH hSrcDS, const char *pszSrcWKT,
    00274                                  GDALDatasetH hDstDS, const char *pszDstWKT,
    00275                                  int bGCPUseOK, double dfGCPErrorThreshold,
    00276                                  int nOrder )
    00277
    00278 {
    00279     GDALGenImgProjTransformInfo *psInfo;
    00280
    00281 /* -------------------------------------------------------------------- */
    00282 /*      Initialize the transform info.                                  */
    00283 /* -------------------------------------------------------------------- */
    00284     psInfo = (GDALGenImgProjTransformInfo *)
    00285         CPLCalloc(sizeof(GDALGenImgProjTransformInfo),1);
    00286
    00287 /* -------------------------------------------------------------------- */
    00288 /*      Get forward and inverse geotransform for the source image.      */
    00289 /* -------------------------------------------------------------------- */
    00290     if( GDALGetGeoTransform( hSrcDS, psInfo->adfSrcGeoTransform ) == CE_None
    00291         && (psInfo->adfSrcGeoTransform[0] != 0.0
    00292             || psInfo->adfSrcGeoTransform[1] != 1.0
    00293             || psInfo->adfSrcGeoTransform[2] != 0.0
    00294             || psInfo->adfSrcGeoTransform[3] != 0.0
    00295             || psInfo->adfSrcGeoTransform[4] != 0.0
    00296             || ABS(psInfo->adfSrcGeoTransform[5]) != 1.0) )
    00297     {
    00298         InvGeoTransform( psInfo->adfSrcGeoTransform,
    00299                          psInfo->adfSrcInvGeoTransform );
    00300     }
    00301     else if( bGCPUseOK && GDALGetGCPCount( hSrcDS ) > 0 )
    00302     {
    00303         psInfo->pSrcGCPTransformArg =
    00304             GDALCreateGCPTransformer( GDALGetGCPCount( hSrcDS ),
    00305                                       GDALGetGCPs( hSrcDS ), nOrder, FALSE );
    00306         if( psInfo->pSrcGCPTransformArg == NULL )
    00307         {
    00308             GDALDestroyGenImgProjTransformer( psInfo );
    00309             return NULL;
    00310         }
    00311     }
    00312     else
    00313     {
    00314         CPLError( CE_Failure, CPLE_AppDefined,
    00315                   "Unable to compute a transformation between pixel/line\n"
    00316                   "and georeferenced coordinates for %s.\n"
    00317                   "There is no affine transformation and no GCPs.",
    00318                   GDALGetDescription( hSrcDS ) );
    00319
    00320         GDALDestroyGenImgProjTransformer( psInfo );
    00321         return NULL;
    00322     }
    00323
    00324 /* -------------------------------------------------------------------- */
    00325 /*      Setup reprojection.                                             */
    00326 /* -------------------------------------------------------------------- */
    00327     if( pszSrcWKT != NULL && pszDstWKT != NULL
    00328         && !EQUAL(pszSrcWKT,pszDstWKT) )
    00329     {
    00330         psInfo->pReprojectArg =
    00331             GDALCreateReprojectionTransformer( pszSrcWKT, pszDstWKT );
    00332     }
    00333
    00334 /* -------------------------------------------------------------------- */
    00335 /*      Get forward and inverse geotransform for destination image.     */
    00336 /*      If we have no destination use a unit transform.                 */
    00337 /* -------------------------------------------------------------------- */
    00338     if( hDstDS )
    00339     {
    00340         GDALGetGeoTransform( hDstDS, psInfo->adfDstGeoTransform );
    00341         InvGeoTransform( psInfo->adfDstGeoTransform,
    00342                          psInfo->adfDstInvGeoTransform );
    00343     }
    00344     else
    00345     {
    00346         psInfo->adfDstGeoTransform[0] = 0.0;
    00347         psInfo->adfDstGeoTransform[1] = 1.0;
    00348         psInfo->adfDstGeoTransform[2] = 0.0;
    00349         psInfo->adfDstGeoTransform[3] = 0.0;
    00350         psInfo->adfDstGeoTransform[4] = 0.0;
    00351         psInfo->adfDstGeoTransform[5] = 1.0;
    00352         memcpy( psInfo->adfDstInvGeoTransform, psInfo->adfDstGeoTransform,
    00353                 sizeof(double) * 6 );
    00354     }
    00355
    00356     return psInfo;
    00357 }
    00358
    00359
    00360 /************************************************************************/
    00361 /*                  GDALDestroyGenImgProjTransformer()                  */
    00362 /************************************************************************/
    00363
    00364 void GDALDestroyGenImgProjTransformer( void *hTransformArg )
    00365
    00366 {
    00367     GDALGenImgProjTransformInfo *psInfo =
    00368         (GDALGenImgProjTransformInfo *) hTransformArg;
    00369
    00370     if( psInfo->pSrcGCPTransformArg != NULL )
    00371         GDALDestroyGCPTransformer( psInfo->pSrcGCPTransformArg );
    00372
    00373     if( psInfo->pDstGCPTransformArg != NULL )
    00374         GDALDestroyGCPTransformer( psInfo->pDstGCPTransformArg );
    00375
    00376     if( psInfo->pReprojectArg != NULL )
    00377         GDALDestroyReprojectionTransformer( psInfo->pReprojectArg );
    00378
    00379     CPLFree( psInfo );
    00380 }
    00381
    00382 /************************************************************************/
    00383 /*                      GDALGenImgProjTransform()                       */
    00384 /************************************************************************/
    00385
    00386 int GDALGenImgProjTransform( void *pTransformArg, int bDstToSrc,
    00387                              int nPointCount,
    00388                              double *padfX, double *padfY, double *padfZ,
    00389                              int *panSuccess )
    00390 {
    00391     GDALGenImgProjTransformInfo *psInfo =
    00392         (GDALGenImgProjTransformInfo *) pTransformArg;
    00393     int   i;
    00394     double *padfGeoTransform;
    00395     void *pGCPTransformArg;
    00396
    00397 /* -------------------------------------------------------------------- */
    00398 /*      Convert from src (dst) pixel/line to src (dst)                  */
    00399 /*      georeferenced coordinates.                                      */
    00400 /* -------------------------------------------------------------------- */
    00401     if( bDstToSrc )
    00402     {
    00403         padfGeoTransform = psInfo->adfDstGeoTransform;
    00404         pGCPTransformArg = psInfo->pDstGCPTransformArg;
    00405     }
    00406     else
    00407     {
    00408         padfGeoTransform = psInfo->adfSrcGeoTransform;
    00409         pGCPTransformArg = psInfo->pSrcGCPTransformArg;
    00410     }
    00411
    00412     if( pGCPTransformArg == NULL )
    00413     {
    00414         for( i = 0; i < nPointCount; i++ )
    00415         {
    00416             double dfNewX, dfNewY;
    00417
    00418             dfNewX = padfGeoTransform[0]
    00419                 + padfX[i] * padfGeoTransform[1]
    00420                 + padfY[i] * padfGeoTransform[2];
    00421             dfNewY = padfGeoTransform[3]
    00422                 + padfX[i] * padfGeoTransform[4]
    00423                 + padfY[i] * padfGeoTransform[5];
    00424
    00425             padfX[i] = dfNewX;
    00426             padfY[i] = dfNewY;
    00427         }
    00428     }
    00429     else
    00430     {
    00431         if( !GDALGCPTransform( pGCPTransformArg, FALSE,
    00432                                nPointCount, padfX, padfY, padfZ,
    00433                                panSuccess ) )
    00434             return FALSE;
    00435     }
    00436
    00437 /* -------------------------------------------------------------------- */
    00438 /*      Reproject if needed.                                            */
    00439 /* -------------------------------------------------------------------- */
    00440     if( psInfo->pReprojectArg )
    00441     {
    00442         if( !GDALReprojectionTransform( psInfo->pReprojectArg, bDstToSrc,
    00443                                         nPointCount, padfX, padfY, padfZ,
    00444                                         panSuccess ) )
    00445             return FALSE;
    00446     }
    00447     else
    00448     {
    00449         for( i = 0; i < nPointCount; i++ )
    00450             panSuccess[i] = 1;
    00451     }
    00452
    00453 /* -------------------------------------------------------------------- */
    00454 /*      Convert dst (src) georef coordinates back to pixel/line.        */
    00455 /* -------------------------------------------------------------------- */
    00456     if( bDstToSrc )
    00457     {
    00458         padfGeoTransform = psInfo->adfSrcInvGeoTransform;
    00459         pGCPTransformArg = psInfo->pSrcGCPTransformArg;
    00460     }
    00461     else
    00462     {
    00463         padfGeoTransform = psInfo->adfDstInvGeoTransform;
    00464         pGCPTransformArg = psInfo->pDstGCPTransformArg;
    00465     }
    00466
    00467     if( pGCPTransformArg == NULL )
    00468     {
    00469         for( i = 0; i < nPointCount; i++ )
    00470         {
    00471             double dfNewX, dfNewY;
    00472
    00473             dfNewX = padfGeoTransform[0]
    00474                 + padfX[i] * padfGeoTransform[1]
    00475                 + padfY[i] * padfGeoTransform[2];
    00476             dfNewY = padfGeoTransform[3]
    00477                 + padfX[i] * padfGeoTransform[4]
    00478                 + padfY[i] * padfGeoTransform[5];
    00479
    00480             padfX[i] = dfNewX;
    00481             padfY[i] = dfNewY;
    00482         }
    00483     }
    00484     else
    00485     {
    00486         if( !GDALGCPTransform( pGCPTransformArg, TRUE,
    00487                                nPointCount, padfX, padfY, padfZ,
    00488                                panSuccess ) )
    00489             return FALSE;
    00490     }
    00491
    00492     return TRUE;
    00493 }
    00494
    00495
    00496 /************************************************************************/
    00497 /* ==================================================================== */
    00498 /*                       GDALReprojectionTransformer                    */
    00499 /* ==================================================================== */
    00500 /************************************************************************/
    00501
    00502 typedef struct {
    00503     OGRCoordinateTransformation *poForwardTransform;
    00504     OGRCoordinateTransformation *poReverseTransform;
    00505 } GDALReprojectionTransformInfo;
    00506
    00507 /************************************************************************/
    00508 /*                 GDALCreateReprojectionTransformer()                  */
    00509 /************************************************************************/
    00510
    00511 void *GDALCreateReprojectionTransformer( const char *pszSrcWKT,
    00512                                          const char *pszDstWKT )
    00513
    00514 {
    00515     OGRSpatialReference oSrcSRS, oDstSRS;
    00516     OGRCoordinateTransformation *poForwardTransform;
    00517
    00518 /* -------------------------------------------------------------------- */
    00519 /*      Ingest the SRS definitions.                                     */
    00520 /* -------------------------------------------------------------------- */
    00521     if( oSrcSRS.importFromWkt( (char **) &pszSrcWKT ) != OGRERR_NONE )
    00522     {
    00523         CPLError( CE_Failure, CPLE_AppDefined,
    00524                   "Failed to import coordinate system `%s'.",
    00525                   pszSrcWKT );
    00526         return NULL;
    00527     }
    00528     if( oDstSRS.importFromWkt( (char **) &pszDstWKT ) != OGRERR_NONE )
    00529     {
    00530         CPLError( CE_Failure, CPLE_AppDefined,
    00531                   "Failed to import coordinate system `%s'.",
    00532                   pszSrcWKT );
    00533         return NULL;
    00534     }
    00535
    00536 /* -------------------------------------------------------------------- */
    00537 /*      Build the forward coordinate transformation.                    */
    00538 /* -------------------------------------------------------------------- */
    00539     poForwardTransform = OGRCreateCoordinateTransformation(&oSrcSRS,&oDstSRS);
    00540
    00541     if( poForwardTransform == NULL )
    00542         // OGRCreateCoordinateTransformation() will report errors on its own.
    00543         return NULL;
    00544
    00545 /* -------------------------------------------------------------------- */
    00546 /*      Create a structure to hold the transform info, and also         */
    00547 /*      build reverse transform.  We assume that if the forward         */
    00548 /*      transform can be created, then so can the reverse one.          */
    00549 /* -------------------------------------------------------------------- */
    00550     GDALReprojectionTransformInfo *psInfo;
    00551
    00552     psInfo = (GDALReprojectionTransformInfo *)
    00553         CPLCalloc(sizeof(GDALReprojectionTransformInfo),1);
    00554
    00555     psInfo->poForwardTransform = poForwardTransform;
    00556     psInfo->poReverseTransform =
    00557         OGRCreateCoordinateTransformation(&oDstSRS,&oSrcSRS);
    00558
    00559     return psInfo;
    00560 }
    00561
    00562 /************************************************************************/
    00563 /*                  GDALDestroyReprojectionTransform()                  */
    00564 /************************************************************************/
    00565
    00566 void GDALDestroyReprojectionTransformer( void *pTransformAlg )
    00567
    00568 {
    00569     GDALReprojectionTransformInfo *psInfo =
    00570         (GDALReprojectionTransformInfo *) pTransformAlg;
    00571
    00572     if( psInfo->poForwardTransform )
    00573         delete psInfo->poForwardTransform;
    00574
    00575     if( psInfo->poReverseTransform )
    00576     delete psInfo->poReverseTransform;
    00577
    00578     CPLFree( psInfo );
    00579 }
    00580
    00581 /************************************************************************/
    00582 /*                     GDALReprojectionTransform()                      */
    00583 /************************************************************************/
    00584
    00585 int GDALReprojectionTransform( void *pTransformArg, int bDstToSrc,
    00586                                 int nPointCount,
    00587                                 double *padfX, double *padfY, double *padfZ,
    00588                                 int *panSuccess )
    00589
    00590 {
    00591     GDALReprojectionTransformInfo *psInfo =
    00592         (GDALReprojectionTransformInfo *) pTransformArg;
    00593     int bSuccess;
    00594
    00595     if( bDstToSrc )
    00596         bSuccess = psInfo->poReverseTransform->Transform(
    00597             nPointCount, padfX, padfY, padfZ );
    00598     else
    00599         bSuccess = psInfo->poForwardTransform->Transform(
    00600             nPointCount, padfX, padfY, padfZ );
    00601
    00602     if( bSuccess )
    00603         memset( panSuccess, 1, sizeof(int) * nPointCount );
    00604     else
    00605         memset( panSuccess, 0, sizeof(int) * nPointCount );
    00606
    00607     return bSuccess;
    00608 }
    00609
    00610
    00611 /************************************************************************/
    00612 /* ==================================================================== */
    00613 /*      Approximate transformer.                                        */
    00614 /* ==================================================================== */
    00615 /************************************************************************/
    00616
    00617 typedef struct
    00618 {
    00619     GDALTransformerFunc pfnBaseTransformer;
    00620     void             *pBaseCBData;
    00621     double            dfMaxError;
    00622 } ApproxTransformInfo;
    00623
    00624 /************************************************************************/
    00625 /*                    GDALCreateApproxTransformer()                     */
    00626 /************************************************************************/
    00627
    00628 void *GDALCreateApproxTransformer( GDALTransformerFunc pfnBaseTransformer,
    00629                                    void *pBaseTransformArg, double dfMaxError)
    00630
    00631 {
    00632     ApproxTransformInfo *psATInfo;
    00633
    00634     psATInfo = (ApproxTransformInfo*) CPLMalloc(sizeof(ApproxTransformInfo));
    00635     psATInfo->pfnBaseTransformer = pfnBaseTransformer;
    00636     psATInfo->pBaseCBData = pBaseTransformArg;
    00637     psATInfo->dfMaxError = dfMaxError;
    00638
    00639     return psATInfo;
    00640 }
    00641
    00642 /************************************************************************/
    00643 /*                    GDALDestroyApproxTransformer()                    */
    00644 /************************************************************************/
    00645
    00646 void GDALDestroyApproxTransformer( void * pCBData )
    00647
    00648 {
    00649     CPLFree( pCBData );
    00650 }
    00651
    00652 /************************************************************************/
    00653 /*                        GDALApproxTransform()                         */
    00654 /************************************************************************/
    00655
    00656 int GDALApproxTransform( void *pCBData, int bDstToSrc, int nPoints,
    00657                          double *x, double *y, double *z, int *panSuccess )
    00658
    00659 {
    00660     ApproxTransformInfo *psATInfo = (ApproxTransformInfo *) pCBData;
    00661     double x2[3], y2[3], z2[3], dfDeltaX, dfDeltaY, dfError, dfDist, dfDeltaZ;
    00662     int nMiddle, anSuccess2[3], i, bSuccess;
    00663
    00664     nMiddle = (nPoints-1)/2;
    00665
    00666 /* -------------------------------------------------------------------- */
    00667 /*      Bail if our preconditions are not met, or if error is not       */
    00668 /*      acceptable.                                                     */
    00669 /* -------------------------------------------------------------------- */
    00670     if( y[0] != y[nPoints-1] || y[0] != y[nMiddle]
    00671         || x[0] == x[nPoints-1] || x[0] == x[nMiddle]
    00672         || psATInfo->dfMaxError == 0.0 || nPoints <= 5 )
    00673     {
    00674         return psATInfo->pfnBaseTransformer( psATInfo->pBaseCBData, bDstToSrc,
    00675                                              nPoints, x, y, z, panSuccess );
    00676     }
    00677
    00678 /* -------------------------------------------------------------------- */
    00679 /*      Transform first, last and middle point.                         */
    00680 /* -------------------------------------------------------------------- */
    00681     x2[0] = x[0];
    00682     y2[0] = y[0];
    00683     z2[0] = z[0];
    00684     x2[1] = x[nMiddle];
    00685     y2[1] = y[nMiddle];
    00686     z2[1] = z[nMiddle];
    00687     x2[2] = x[nPoints-1];
    00688     y2[2] = y[nPoints-1];
    00689     z2[2] = z[nPoints-1];
    00690
    00691     bSuccess =
    00692         psATInfo->pfnBaseTransformer( psATInfo->pBaseCBData, bDstToSrc, 3,
    00693                                       x2, y2, z2, anSuccess2 );
    00694     if( !bSuccess )
    00695         return psATInfo->pfnBaseTransformer( psATInfo->pBaseCBData, bDstToSrc,
    00696                                              nPoints, x, y, z, panSuccess );
    00697
    00698 /* -------------------------------------------------------------------- */
    00699 /*      Is the error at the middle acceptable relative to an            */
    00700 /*      interpolation of the middle position?                           */
    00701 /* -------------------------------------------------------------------- */
    00702     dfDeltaX = (x2[2] - x2[0]) / (x[nPoints-1] - x[0]);
    00703     dfDeltaY = (y2[2] - y2[0]) / (x[nPoints-1] - x[0]);
    00704     dfDeltaZ = (z2[2] - z2[0]) / (x[nPoints-1] - x[0]);
    00705
    00706     dfError = fabs((x2[0] + dfDeltaX * (x[nMiddle] - x[0])) - x2[1])
    00707         + fabs((y2[0] + dfDeltaY * (x[nMiddle] - x[0])) - y2[1]);
    00708
    00709     if( dfError > psATInfo->dfMaxError )
    00710     {
    00711 #ifdef notdef
    00712         CPLDebug( "GDAL", "ApproxTransformer - "
    00713                   "error %g over threshold %g, subdivide %d points.",
    00714                   dfError, psATInfo->dfMaxError, nPoints );
    00715 #endif
    00716
    00717         bSuccess =
    00718             GDALApproxTransform( psATInfo, bDstToSrc, nMiddle,
    00719                                  x, y, z, panSuccess );
    00720
    00721         if( !bSuccess )
    00722             return FALSE;
    00723
    00724         bSuccess =
    00725             GDALApproxTransform( psATInfo, bDstToSrc, nPoints - nMiddle,
    00726                                  x+nMiddle, y+nMiddle, z+nMiddle,
    00727                                  panSuccess+nMiddle );
    00728
    00729         if( !bSuccess )
    00730             return FALSE;
    00731
    00732         return TRUE;
    00733     }
    00734
    00735 /* -------------------------------------------------------------------- */
    00736 /*      Error is OK, linearly interpolate all points along line.        */
    00737 /* -------------------------------------------------------------------- */
    00738     for( i = nPoints-1; i >= 0; i-- )
    00739     {
    00740         dfDist = (x[i] - x[0]);
    00741         y[i] = y2[0] + dfDeltaY * dfDist;
    00742         x[i] = x2[0] + dfDeltaX * dfDist;
    00743         z[i] = z2[0] + dfDeltaZ * dfDist;
    00744         panSuccess[i] = TRUE;
    00745     }
    00746
    00747     return TRUE;
    00748 }
    00749