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