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 : /* -------------------------------------------------------------------- */
229 : /* If no bands provided assume we should process all bands. */
230 : /* -------------------------------------------------------------------- */
231 2 : if (nBandCount == 0)
232 : {
233 1 : nBandCount = GDALGetRasterCount(hSrcDS);
234 1 : if (nBandCount == 0)
235 : {
236 0 : CPLError(CE_Failure, CPLE_AppDefined,
237 : "No raster band in source dataset");
238 0 : return FALSE;
239 : }
240 :
241 1 : panBandList = static_cast<int *>(CPLCalloc(sizeof(int), nBandCount));
242 :
243 2 : for (int iBand = 0; iBand < nBandCount; iBand++)
244 1 : panBandList[iBand] = iBand + 1;
245 :
246 1 : const int nResult = GDALSimpleImageWarp(
247 : hSrcDS, hDstDS, nBandCount, panBandList, pfnTransform,
248 : pTransformArg, pfnProgress, pProgressArg, papszWarpOptions);
249 1 : CPLFree(panBandList);
250 1 : return nResult;
251 : }
252 :
253 : /* -------------------------------------------------------------------- */
254 : /* Post initial progress. */
255 : /* -------------------------------------------------------------------- */
256 1 : if (pfnProgress)
257 : {
258 0 : if (!pfnProgress(0.0, "", pProgressArg))
259 0 : return FALSE;
260 : }
261 :
262 : /* -------------------------------------------------------------------- */
263 : /* Load the source image band(s). */
264 : /* -------------------------------------------------------------------- */
265 1 : const int nSrcXSize = GDALGetRasterXSize(hSrcDS);
266 1 : const int nSrcYSize = GDALGetRasterYSize(hSrcDS);
267 : std::vector<std::unique_ptr<GByte, VSIFreeReleaser>> apabySrcDataUniquePtr(
268 2 : nBandCount);
269 :
270 2 : for (int iBand = 0; iBand < nBandCount; iBand++)
271 : {
272 1 : apabySrcDataUniquePtr[iBand].reset(
273 1 : static_cast<GByte *>(VSI_MALLOC2_VERBOSE(nSrcXSize, nSrcYSize)));
274 1 : if (apabySrcDataUniquePtr[iBand] == nullptr)
275 : {
276 0 : CPLError(CE_Failure, CPLE_OutOfMemory,
277 : "GDALSimpleImageWarp out of memory.");
278 0 : return FALSE;
279 : }
280 :
281 1 : if (GDALRasterIO(GDALGetRasterBand(hSrcDS, panBandList[iBand]), GF_Read,
282 : 0, 0, nSrcXSize, nSrcYSize,
283 1 : apabySrcDataUniquePtr[iBand].get(), nSrcXSize,
284 1 : nSrcYSize, GDT_Byte, 0, 0) != CE_None)
285 : {
286 0 : CPLError(CE_Failure, CPLE_FileIO,
287 : "GDALSimpleImageWarp GDALRasterIO failure %s",
288 : CPLGetLastErrorMsg());
289 0 : return FALSE;
290 : }
291 : }
292 :
293 : /* -------------------------------------------------------------------- */
294 : /* Check for remap request(s). */
295 : /* -------------------------------------------------------------------- */
296 2 : std::vector<GByte *> apabySrcData;
297 2 : for (auto &ptr : apabySrcDataUniquePtr)
298 1 : apabySrcData.push_back(ptr.get());
299 1 : GDALSimpleWarpRemapping(nBandCount, apabySrcData.data(), nSrcXSize,
300 : nSrcYSize, papszWarpOptions);
301 :
302 : /* -------------------------------------------------------------------- */
303 : /* Allocate scanline buffers for output image. */
304 : /* -------------------------------------------------------------------- */
305 1 : const int nDstXSize = GDALGetRasterXSize(hDstDS);
306 1 : const int nDstYSize = GDALGetRasterYSize(hDstDS);
307 : std::vector<std::unique_ptr<GByte, VSIFreeReleaser>> apabyDstLine(
308 2 : nBandCount);
309 :
310 2 : for (int iBand = 0; iBand < nBandCount; iBand++)
311 : {
312 1 : apabyDstLine[iBand].reset(
313 1 : static_cast<GByte *>(VSI_MALLOC_VERBOSE(nDstXSize)));
314 1 : if (!apabyDstLine[iBand])
315 0 : return FALSE;
316 : }
317 :
318 : /* -------------------------------------------------------------------- */
319 : /* Allocate x,y,z coordinate arrays for transformation ... one */
320 : /* scanlines worth of positions. */
321 : /* -------------------------------------------------------------------- */
322 : std::unique_ptr<double, VSIFreeReleaser> padfX(
323 2 : static_cast<double *>(VSI_MALLOC2_VERBOSE(sizeof(double), nDstXSize)));
324 : std::unique_ptr<double, VSIFreeReleaser> padfY(
325 2 : static_cast<double *>(VSI_MALLOC2_VERBOSE(sizeof(double), nDstXSize)));
326 : std::unique_ptr<double, VSIFreeReleaser> padfZ(
327 2 : static_cast<double *>(VSI_MALLOC2_VERBOSE(sizeof(double), nDstXSize)));
328 : std::unique_ptr<int, VSIFreeReleaser> pabSuccess(
329 2 : static_cast<int *>(VSI_MALLOC2_VERBOSE(sizeof(int), nDstXSize)));
330 1 : if (!padfX || !padfY || !padfZ || !pabSuccess)
331 0 : return FALSE;
332 :
333 : /* -------------------------------------------------------------------- */
334 : /* Establish the value we will use to initialize the bands. We */
335 : /* default to -1 indicating the initial value should be read */
336 : /* and preserved from the source file, but allow this to be */
337 : /* overridden by passed */
338 : /* option(s). */
339 : /* -------------------------------------------------------------------- */
340 2 : std::vector<int> anBandInit(nBandCount);
341 1 : if (CSLFetchNameValue(papszWarpOptions, "INIT"))
342 : {
343 : const CPLStringList aosTokens(CSLTokenizeStringComplex(
344 0 : CSLFetchNameValue(papszWarpOptions, "INIT"), " ,", FALSE, FALSE));
345 :
346 0 : const int nTokenCount = aosTokens.size();
347 :
348 0 : for (int iBand = 0; iBand < nBandCount; iBand++)
349 : {
350 0 : if (nTokenCount == 0)
351 0 : anBandInit[iBand] = 0;
352 : else
353 0 : anBandInit[iBand] =
354 0 : atoi(aosTokens[std::min(iBand, nTokenCount - 1)]);
355 : }
356 : }
357 :
358 : /* -------------------------------------------------------------------- */
359 : /* Loop over all the scanlines in the output image. */
360 : /* -------------------------------------------------------------------- */
361 21 : for (int iDstY = 0; iDstY < nDstYSize; iDstY++)
362 : {
363 : // Clear output buffer to "transparent" value. Should not we
364 : // really be reading from the destination file to support overlay?
365 40 : for (int iBand = 0; iBand < nBandCount; iBand++)
366 : {
367 20 : if (anBandInit[iBand] == -1)
368 : {
369 0 : if (GDALRasterIO(GDALGetRasterBand(hDstDS, iBand + 1), GF_Read,
370 : 0, iDstY, nDstXSize, 1,
371 0 : apabyDstLine[iBand].get(), nDstXSize, 1,
372 0 : GDT_Byte, 0, 0) != CE_None)
373 : {
374 0 : return FALSE;
375 : }
376 : }
377 : else
378 : {
379 20 : memset(apabyDstLine[iBand].get(), anBandInit[iBand], nDstXSize);
380 : }
381 : }
382 :
383 : // Set point to transform.
384 420 : for (int iDstX = 0; iDstX < nDstXSize; iDstX++)
385 : {
386 400 : (padfX.get())[iDstX] = iDstX + 0.5;
387 400 : (padfY.get())[iDstX] = iDstY + 0.5;
388 400 : (padfZ.get())[iDstX] = 0.0;
389 : }
390 :
391 : // Transform the points from destination pixel/line coordinates
392 : // to source pixel/line coordinates.
393 20 : pfnTransform(pTransformArg, TRUE, nDstXSize, padfX.get(), padfY.get(),
394 : padfZ.get(), pabSuccess.get());
395 :
396 : // Loop over the output scanline.
397 420 : for (int iDstX = 0; iDstX < nDstXSize; iDstX++)
398 : {
399 400 : if (!(pabSuccess.get())[iDstX])
400 0 : continue;
401 :
402 : // We test against the value before casting to avoid the
403 : // problem of asymmetric truncation effects around zero. That is
404 : // -0.5 will be 0 when cast to an int.
405 400 : if ((padfX.get())[iDstX] < 0.0 || (padfY.get())[iDstX] < 0.0)
406 0 : continue;
407 :
408 400 : const int iSrcX = static_cast<int>((padfX.get())[iDstX]);
409 400 : const int iSrcY = static_cast<int>((padfY.get())[iDstX]);
410 :
411 400 : if (iSrcX >= nSrcXSize || iSrcY >= nSrcYSize)
412 0 : continue;
413 :
414 400 : const int iSrcOffset = iSrcX + iSrcY * nSrcXSize;
415 :
416 800 : for (int iBand = 0; iBand < nBandCount; iBand++)
417 400 : (apabyDstLine[iBand].get())[iDstX] =
418 400 : apabySrcData[iBand][iSrcOffset];
419 : }
420 :
421 : // Write scanline to disk.
422 40 : for (int iBand = 0; iBand < nBandCount; iBand++)
423 : {
424 20 : if (GDALRasterIO(GDALGetRasterBand(hDstDS, iBand + 1), GF_Write, 0,
425 20 : iDstY, nDstXSize, 1, apabyDstLine[iBand].get(),
426 20 : nDstXSize, 1, GDT_Byte, 0, 0) != CE_None)
427 : {
428 0 : return FALSE;
429 : }
430 : }
431 :
432 20 : if (pfnProgress != nullptr)
433 : {
434 0 : if (!pfnProgress((iDstY + 1) / static_cast<double>(nDstYSize), "",
435 : pProgressArg))
436 : {
437 0 : CPLError(CE_Failure, CPLE_UserInterrupt, "User terminated");
438 0 : return FALSE;
439 : }
440 : }
441 : }
442 :
443 1 : return TRUE;
444 : }
|
