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