1 : /******************************************************************************
2 : * $Id: gdalsimplewarp.cpp 15436 2008-09-24 19:26:31Z rouault $
3 : *
4 : * Project: Mapinfo Image Warper
5 : * Purpose: Simple (source in memory) warp algorithm.
6 : * Author: Frank Warmerdam, warmerdam@pobox.com
7 : *
8 : ******************************************************************************
9 : * Copyright (c) 2002, i3 - information integration and imaging, Fort Collin,CO
10 : *
11 : * Permission is hereby granted, free of charge, to any person obtaining a
12 : * copy of this software and associated documentation files (the "Software"),
13 : * to deal in the Software without restriction, including without limitation
14 : * the rights to use, copy, modify, merge, publish, distribute, sublicense,
15 : * and/or sell copies of the Software, and to permit persons to whom the
16 : * Software is furnished to do so, subject to the following conditions:
17 : *
18 : * The above copyright notice and this permission notice shall be included
19 : * in all copies or substantial portions of the Software.
20 : *
21 : * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
22 : * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
23 : * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
24 : * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
25 : * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
26 : * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
27 : * DEALINGS IN THE SOFTWARE.
28 : ****************************************************************************/
29 :
30 : #include "gdal_priv.h"
31 : #include "gdal_alg.h"
32 : #include "cpl_string.h"
33 :
34 : CPL_CVSID("$Id: gdalsimplewarp.cpp 15436 2008-09-24 19:26:31Z rouault $");
35 :
36 : static void
37 : GDALSimpleWarpRemapping( int nBandCount, GByte **papabySrcData,
38 : int nSrcXSize, int nSrcYSize,
39 : char **papszWarpOptions );
40 :
41 : /************************************************************************/
42 : /* GDALSimpleImageWarp() */
43 : /************************************************************************/
44 :
45 : /**
46 : * Perform simple image warp.
47 : *
48 : * Copies an image from a source dataset to a destination dataset applying
49 : * an application defined transformation. This algorithm is called simple
50 : * because it lacks many options such as resampling kernels (other than
51 : * nearest neighbour), support for data types other than 8bit, and the
52 : * ability to warp images without holding the entire source and destination
53 : * image in memory.
54 : *
55 : * The following option(s) may be passed in papszWarpOptions.
56 : * <ul>
57 : * <li> "INIT=v[,v...]": This option indicates that the output dataset should
58 : * be initialized to the indicated value in any area valid data is not written.
59 : * Distinct values may be listed for each band separated by columns.
60 : * </ul>
61 : *
62 : * @param hSrcDS the source image dataset.
63 : * @param hDstDS the destination image dataset.
64 : * @param nBandCount the number of bands to be warped. If zero, all bands
65 : * will be processed.
66 : * @param panBandList the list of bands to translate.
67 : * @param pfnTransform the transformation function to call. See
68 : * GDALTransformerFunc().
69 : * @param pTransformArg the callback handle to pass to pfnTransform.
70 : * @param pfnProgress the function used to report progress. See
71 : * GDALProgressFunc().
72 : * @param pProgressArg the callback handle to pass to pfnProgress.
73 : * @param papszWarpOptions additional options controlling the warp.
74 : *
75 : * @return TRUE if the operation completes, or FALSE if an error occurs.
76 : */
77 :
78 : int CPL_STDCALL
79 0 : GDALSimpleImageWarp( GDALDatasetH hSrcDS, GDALDatasetH hDstDS,
80 : int nBandCount, int *panBandList,
81 : GDALTransformerFunc pfnTransform, void *pTransformArg,
82 : GDALProgressFunc pfnProgress, void *pProgressArg,
83 : char **papszWarpOptions )
84 :
85 : {
86 0 : VALIDATE_POINTER1( hSrcDS, "GDALSimpleImageWarp", 0 );
87 0 : VALIDATE_POINTER1( hDstDS, "GDALSimpleImageWarp", 0 );
88 :
89 0 : int iBand, bCancelled = FALSE;
90 :
91 : /* -------------------------------------------------------------------- */
92 : /* If no bands provided assume we should process all bands. */
93 : /* -------------------------------------------------------------------- */
94 0 : if( nBandCount == 0 )
95 : {
96 : int nResult;
97 :
98 0 : nBandCount = GDALGetRasterCount( hSrcDS );
99 0 : if (nBandCount == 0)
100 : {
101 : CPLError(CE_Failure, CPLE_AppDefined,
102 0 : "No raster band in source dataset");
103 0 : return FALSE;
104 : }
105 :
106 0 : panBandList = (int *) CPLCalloc(sizeof(int),nBandCount);
107 :
108 0 : for( iBand = 0; iBand < nBandCount; iBand++ )
109 0 : panBandList[iBand] = iBand+1;
110 :
111 : nResult = GDALSimpleImageWarp( hSrcDS, hDstDS, nBandCount, panBandList,
112 : pfnTransform, pTransformArg,
113 : pfnProgress, pProgressArg,
114 0 : papszWarpOptions );
115 0 : CPLFree( panBandList );
116 0 : return nResult;
117 : }
118 :
119 : /* -------------------------------------------------------------------- */
120 : /* Post initial progress. */
121 : /* -------------------------------------------------------------------- */
122 0 : if( pfnProgress )
123 : {
124 0 : if( !pfnProgress( 0.0, "", pProgressArg ) )
125 0 : return FALSE;
126 : }
127 :
128 : /* -------------------------------------------------------------------- */
129 : /* Load the source image band(s). */
130 : /* -------------------------------------------------------------------- */
131 0 : int nSrcXSize = GDALGetRasterXSize(hSrcDS);
132 0 : int nSrcYSize = GDALGetRasterYSize(hSrcDS);
133 : GByte **papabySrcData;
134 :
135 0 : papabySrcData = (GByte **) CPLCalloc(nBandCount,sizeof(GByte*));
136 0 : for( iBand = 0; iBand < nBandCount; iBand++ )
137 : {
138 0 : papabySrcData[iBand] = (GByte *) VSIMalloc(nSrcXSize*nSrcYSize);
139 :
140 : GDALRasterIO( GDALGetRasterBand(hSrcDS,panBandList[iBand]), GF_Read,
141 : 0, 0, nSrcXSize, nSrcYSize,
142 0 : papabySrcData[iBand], nSrcXSize, nSrcYSize, GDT_Byte,
143 0 : 0, 0 );
144 : }
145 :
146 : /* -------------------------------------------------------------------- */
147 : /* Check for remap request(s). */
148 : /* -------------------------------------------------------------------- */
149 : GDALSimpleWarpRemapping( nBandCount, papabySrcData, nSrcXSize, nSrcYSize,
150 0 : papszWarpOptions );
151 :
152 : /* -------------------------------------------------------------------- */
153 : /* Allocate scanline buffers for output image. */
154 : /* -------------------------------------------------------------------- */
155 0 : int nDstXSize = GDALGetRasterXSize( hDstDS );
156 0 : int nDstYSize = GDALGetRasterYSize( hDstDS );
157 : GByte **papabyDstLine;
158 :
159 0 : papabyDstLine = (GByte **) CPLCalloc(nBandCount,sizeof(GByte*));
160 :
161 0 : for( iBand = 0; iBand < nBandCount; iBand++ )
162 0 : papabyDstLine[iBand] = (GByte *) CPLMalloc( nDstXSize );
163 :
164 : /* -------------------------------------------------------------------- */
165 : /* Allocate x,y,z coordinate arrays for transformation ... one */
166 : /* scanlines worth of positions. */
167 : /* -------------------------------------------------------------------- */
168 : double *padfX, *padfY, *padfZ;
169 : int *pabSuccess;
170 :
171 0 : padfX = (double *) CPLMalloc(sizeof(double) * nDstXSize);
172 0 : padfY = (double *) CPLMalloc(sizeof(double) * nDstXSize);
173 0 : padfZ = (double *) CPLMalloc(sizeof(double) * nDstXSize);
174 0 : pabSuccess = (int *) CPLMalloc(sizeof(int) * nDstXSize);
175 :
176 : /* -------------------------------------------------------------------- */
177 : /* Establish the value we will use to initialize the bands. We */
178 : /* default to -1 indicating the initial value should be read */
179 : /* and preserved from the source file, but allow this to be */
180 : /* overridden by passed */
181 : /* option(s). */
182 : /* -------------------------------------------------------------------- */
183 : int *panBandInit;
184 :
185 0 : panBandInit = (int *) CPLCalloc(sizeof(int),nBandCount);
186 0 : if( CSLFetchNameValue( papszWarpOptions, "INIT" ) )
187 : {
188 : int iBand, nTokenCount;
189 : char **papszTokens =
190 : CSLTokenizeStringComplex( CSLFetchNameValue( papszWarpOptions,
191 : "INIT" ),
192 0 : " ,", FALSE, FALSE );
193 :
194 0 : nTokenCount = CSLCount(papszTokens);
195 :
196 0 : for( iBand = 0; iBand < nBandCount; iBand++ )
197 : {
198 0 : if( nTokenCount == 0 )
199 0 : panBandInit[iBand] = 0;
200 : else
201 0 : panBandInit[iBand] =
202 0 : atoi(papszTokens[MIN(iBand,nTokenCount-1)]);
203 : }
204 :
205 0 : CSLDestroy(papszTokens);
206 : }
207 :
208 : /* -------------------------------------------------------------------- */
209 : /* Loop over all the scanlines in the output image. */
210 : /* -------------------------------------------------------------------- */
211 : int iDstY;
212 :
213 0 : for( iDstY = 0; iDstY < nDstYSize; iDstY++ )
214 : {
215 : int iDstX;
216 :
217 : // Clear output buffer to "transparent" value. Shouldn't we
218 : // really be reading from the destination file to support overlay?
219 0 : for( iBand = 0; iBand < nBandCount; iBand++ )
220 : {
221 0 : if( panBandInit[iBand] == -1 )
222 : GDALRasterIO( GDALGetRasterBand(hDstDS,iBand+1), GF_Read,
223 : 0, iDstY, nDstXSize, 1,
224 0 : papabyDstLine[iBand], nDstXSize, 1, GDT_Byte,
225 0 : 0, 0 );
226 : else
227 0 : memset( papabyDstLine[iBand], panBandInit[iBand], nDstXSize );
228 : }
229 :
230 : // Set point to transform.
231 0 : for( iDstX = 0; iDstX < nDstXSize; iDstX++ )
232 : {
233 0 : padfX[iDstX] = iDstX + 0.5;
234 0 : padfY[iDstX] = iDstY + 0.5;
235 0 : padfZ[iDstX] = 0.0;
236 : }
237 :
238 : // Transform the points from destination pixel/line coordinates
239 : // to source pixel/line coordinates.
240 : pfnTransform( pTransformArg, TRUE, nDstXSize,
241 0 : padfX, padfY, padfZ, pabSuccess );
242 :
243 : // Loop over the output scanline.
244 0 : for( iDstX = 0; iDstX < nDstXSize; iDstX++ )
245 : {
246 0 : if( !pabSuccess[iDstX] )
247 0 : continue;
248 :
249 : // We test against the value before casting to avoid the
250 : // problem of asymmetric truncation effects around zero. That is
251 : // -0.5 will be 0 when cast to an int.
252 0 : if( padfX[iDstX] < 0.0 || padfY[iDstX] < 0.0 )
253 0 : continue;
254 :
255 : int iSrcX, iSrcY, iSrcOffset;
256 :
257 0 : iSrcX = (int) padfX[iDstX];
258 0 : iSrcY = (int) padfY[iDstX];
259 :
260 0 : if( iSrcX >= nSrcXSize || iSrcY >= nSrcYSize )
261 0 : continue;
262 :
263 0 : iSrcOffset = iSrcX + iSrcY * nSrcXSize;
264 :
265 0 : for( iBand = 0; iBand < nBandCount; iBand++ )
266 0 : papabyDstLine[iBand][iDstX] = papabySrcData[iBand][iSrcOffset];
267 : }
268 :
269 : // Write scanline to disk.
270 0 : for( iBand = 0; iBand < nBandCount; iBand++ )
271 : {
272 : GDALRasterIO( GDALGetRasterBand(hDstDS,iBand+1), GF_Write,
273 : 0, iDstY, nDstXSize, 1,
274 0 : papabyDstLine[iBand], nDstXSize, 1, GDT_Byte, 0, 0 );
275 : }
276 :
277 0 : if( pfnProgress != NULL )
278 : {
279 0 : if( !pfnProgress( (iDstY+1) / (double) nDstYSize,
280 : "", pProgressArg ) )
281 : {
282 0 : CPLError( CE_Failure, CPLE_UserInterrupt, "User terminated" );
283 0 : bCancelled = TRUE;
284 0 : break;
285 : }
286 : }
287 : }
288 :
289 : /* -------------------------------------------------------------------- */
290 : /* Cleanup working buffers. */
291 : /* -------------------------------------------------------------------- */
292 0 : for( iBand = 0; iBand < nBandCount; iBand++ )
293 : {
294 0 : CPLFree( papabyDstLine[iBand] );
295 0 : CPLFree( papabySrcData[iBand] );
296 : }
297 :
298 0 : CPLFree( panBandInit );
299 0 : CPLFree( papabyDstLine );
300 0 : CPLFree( papabySrcData );
301 0 : CPLFree( padfX );
302 0 : CPLFree( padfY );
303 0 : CPLFree( padfZ );
304 0 : CPLFree( pabSuccess );
305 :
306 0 : return !bCancelled;
307 : }
308 :
309 : /************************************************************************/
310 : /* GDALSimpleWarpRemapping() */
311 : /* */
312 : /* This function implements any raster remapping requested in */
313 : /* the options list. The remappings are applied to the source */
314 : /* data before warping. Two kinds are support ... REMAP */
315 : /* commands which remap selected pixel values for any band and */
316 : /* REMAP_MULTI which only remap pixels matching the input in */
317 : /* all bands at once (ie. to remap an RGB value to another). */
318 : /************************************************************************/
319 :
320 : static void
321 0 : GDALSimpleWarpRemapping( int nBandCount, GByte **papabySrcData,
322 : int nSrcXSize, int nSrcYSize,
323 : char **papszWarpOptions )
324 :
325 : {
326 :
327 : /* ==================================================================== */
328 : /* Process any and all single value REMAP commands. */
329 : /* ==================================================================== */
330 : int iRemap;
331 : char **papszRemaps = CSLFetchNameValueMultiple( papszWarpOptions,
332 0 : "REMAP" );
333 :
334 0 : for( iRemap = 0; iRemap < CSLCount(papszRemaps); iRemap++ )
335 : {
336 :
337 : /* -------------------------------------------------------------------- */
338 : /* What are the pixel values to map from and to? */
339 : /* -------------------------------------------------------------------- */
340 0 : char **papszTokens = CSLTokenizeString( papszRemaps[iRemap] );
341 : int nFromValue, nToValue;
342 :
343 0 : if( CSLCount(papszTokens) != 2 )
344 : {
345 : CPLError( CE_Warning, CPLE_AppDefined,
346 : "Ill formed REMAP `%s' ignored in GDALSimpleWarpRemapping()",
347 0 : papszRemaps[iRemap] );
348 0 : continue;
349 : }
350 :
351 0 : nFromValue = atoi(papszTokens[0]);
352 0 : nToValue = atoi(papszTokens[1]);
353 :
354 0 : CSLDestroy( papszTokens );
355 :
356 : /* -------------------------------------------------------------------- */
357 : /* Pass over each band searching for matches. */
358 : /* -------------------------------------------------------------------- */
359 0 : for( int iBand = 0; iBand < nBandCount; iBand++ )
360 : {
361 0 : GByte *pabyData = papabySrcData[iBand];
362 0 : int nPixelCount = nSrcXSize * nSrcYSize;
363 :
364 0 : while( nPixelCount != 0 )
365 : {
366 0 : if( *pabyData == nFromValue )
367 0 : *pabyData = (GByte) nToValue;
368 :
369 0 : pabyData++;
370 0 : nPixelCount--;
371 : }
372 : }
373 : }
374 :
375 0 : CSLDestroy( papszRemaps );
376 :
377 : /* ==================================================================== */
378 : /* Process any and all REMAP_MULTI commands. */
379 : /* ==================================================================== */
380 : papszRemaps = CSLFetchNameValueMultiple( papszWarpOptions,
381 0 : "REMAP_MULTI" );
382 :
383 0 : for( iRemap = 0; iRemap < CSLCount(papszRemaps); iRemap++ )
384 : {
385 : /* -------------------------------------------------------------------- */
386 : /* What are the pixel values to map from and to? */
387 : /* -------------------------------------------------------------------- */
388 0 : char **papszTokens = CSLTokenizeString( papszRemaps[iRemap] );
389 : int *panFromValue, *panToValue;
390 : int nMapBandCount, iBand;
391 :
392 0 : if( CSLCount(papszTokens) % 2 == 1
393 : || CSLCount(papszTokens) == 0
394 : || CSLCount(papszTokens) > nBandCount * 2 )
395 : {
396 : CPLError( CE_Warning, CPLE_AppDefined,
397 : "Ill formed REMAP_MULTI `%s' ignored in GDALSimpleWarpRemapping()",
398 0 : papszRemaps[iRemap] );
399 0 : continue;
400 : }
401 :
402 0 : nMapBandCount = CSLCount(papszTokens) / 2;
403 :
404 0 : panFromValue = (int *) CPLMalloc(sizeof(int) * nMapBandCount );
405 0 : panToValue = (int *) CPLMalloc(sizeof(int) * nMapBandCount );
406 :
407 0 : for( iBand = 0; iBand < nMapBandCount; iBand++ )
408 : {
409 0 : panFromValue[iBand] = atoi(papszTokens[iBand]);
410 0 : panToValue[iBand] = atoi(papszTokens[iBand+nMapBandCount]);
411 : }
412 :
413 0 : CSLDestroy( papszTokens );
414 :
415 : /* -------------------------------------------------------------------- */
416 : /* Search for matching values to replace. */
417 : /* -------------------------------------------------------------------- */
418 0 : int nPixelCount = nSrcXSize * nSrcYSize;
419 : int iPixel;
420 :
421 0 : for( iPixel = 0; iPixel < nPixelCount; iPixel++ )
422 : {
423 0 : if( papabySrcData[0][iPixel] != panFromValue[0] )
424 0 : continue;
425 :
426 0 : int bMatch = TRUE;
427 :
428 0 : for( iBand = 1; iBand < nMapBandCount; iBand++ )
429 : {
430 0 : if( papabySrcData[iBand][iPixel] != panFromValue[iBand] )
431 0 : bMatch = FALSE;
432 : }
433 :
434 0 : if( !bMatch )
435 0 : continue;
436 :
437 0 : for( iBand = 0; iBand < nMapBandCount; iBand++ )
438 0 : papabySrcData[iBand][iPixel] = (GByte) panToValue[iBand];
439 : }
440 :
441 0 : CPLFree( panFromValue );
442 0 : CPLFree( panToValue );
443 : }
444 :
445 0 : CSLDestroy( papszRemaps );
446 0 : }
|
