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 

Generated at Sat Dec 21 14:02:00 2002 for GDAL by doxygen1.2.3-20001105 written by Dimitri van Heesch, © 1997-2000