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