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