Line data Source code
1 : /*
2 : * Copyright 2014-2021 Esri
3 : * Licensed under the Apache License, Version 2.0 (the "License");
4 : * you may not use this file except in compliance with the License.
5 : * You may obtain a copy of the License at
6 : *
7 : * http://www.apache.org/licenses/LICENSE-2.0
8 : *
9 : * Unless required by applicable law or agreed to in writing, software
10 : * distributed under the License is distributed on an "AS IS" BASIS,
11 : * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
12 : * See the License for the specific language governing permissions and
13 : * limitations under the License.
14 : */
15 :
16 : /******************************************************************************
17 : *
18 : * Project: Meta Raster File Format Driver Implementation, overlay support
19 : * Purpose: Implementation overlay support for MRF
20 : *
21 : * Author: Lucian Plesea, Lucian.Plesea jpl.nasa.gov, lplesea esri.com
22 : *
23 : ******************************************************************************
24 : * This source file contains the non GDAL standard part of the MRF overview
25 : *building The PatchOverview method only handles powers of 2 overviews!!
26 : ****************************************************************************/
27 :
28 : #include "marfa.h"
29 : #include <vector>
30 :
31 : NAMESPACE_MRF_START
32 :
33 : //
34 : // Scales by 2x2 a buffer in place, using Nearest resampling
35 : // Always pick the top-left corner
36 : //
37 19 : template <typename T> static void NearByFour(T *buff, int xsz, int ysz)
38 : {
39 19 : T *obuff = buff;
40 115 : for (int line = 0; line < ysz; line++)
41 : {
42 : // Copy every other pixel
43 640 : for (int col = 0; col < xsz; col++, buff++)
44 : {
45 544 : *obuff++ = *buff++;
46 : }
47 : // Skip every other line
48 96 : buff += xsz * 2;
49 : }
50 19 : }
51 :
52 : //
53 : // If the NoData value exists, pick a valid pixel if possible
54 : //
55 7 : template <typename T> static void NearByFour(T *buff, int xsz, int ysz, T ndv)
56 : {
57 7 : T *obuff = buff;
58 7 : T *evenline = buff;
59 :
60 77 : for (int line = 0; line < ysz; line++)
61 : {
62 70 : T *oddline = evenline + xsz * 2;
63 770 : for (int col = 0; col < xsz; col++)
64 : {
65 :
66 700 : if (evenline[0] != ndv)
67 623 : *obuff++ = evenline[0];
68 77 : else if (evenline[1] != ndv)
69 77 : *obuff++ = evenline[1];
70 0 : else if (oddline[0] != ndv)
71 0 : *obuff++ = oddline[0];
72 : else
73 0 : *obuff++ = oddline[1];
74 :
75 700 : evenline += 2;
76 700 : oddline += 2;
77 : }
78 70 : evenline += xsz * 2; // Skips the other input line
79 : }
80 7 : }
81 :
82 : // Scales by 2x2 using averaging
83 : // There are lots of these AverageByFour templates, because some types have to
84 : // be treated slightly different than others. Some could be folded by using
85 : // is_integral(), but support is not universal There are two categories,
86 : // depending on NoData presence
87 : //
88 :
89 : // Integer data types shorter than 32 bit use integer math safely
90 3 : template <typename T> static void AverageByFour(T *buff, int xsz, int ysz)
91 : {
92 3 : T *obuff = buff;
93 3 : T *evenline = buff;
94 :
95 33 : for (int line = 0; line < ysz; line++)
96 : {
97 30 : T *oddline = evenline + xsz * 2;
98 330 : for (int col = 0; col < xsz; col++)
99 : {
100 300 : *obuff++ =
101 300 : (2 + evenline[0] + evenline[1] + oddline[0] + oddline[1]) / 4;
102 300 : evenline += 2;
103 300 : oddline += 2;
104 : }
105 30 : evenline += xsz * 2; // Skips the other line
106 : }
107 3 : }
108 :
109 : // 32bit int specialization, avoiding overflow by using 64bit int math
110 1 : template <> void AverageByFour<GInt32>(GInt32 *buff, int xsz, int ysz)
111 : {
112 1 : GInt32 *obuff = buff;
113 1 : GInt32 *evenline = buff;
114 :
115 11 : for (int line = 0; line < ysz; line++)
116 : {
117 10 : GInt32 *oddline = evenline + xsz * 2;
118 110 : for (int col = 0; col < xsz; col++)
119 : {
120 100 : *obuff++ = (GIntBig(2) + evenline[0] + evenline[1] + oddline[0] +
121 100 : oddline[1]) /
122 : 4;
123 100 : evenline += 2;
124 100 : oddline += 2;
125 : }
126 10 : evenline += xsz * 2; // Skips the other line
127 : }
128 1 : }
129 :
130 : // Same for 32bit unsigned int specialization
131 1 : template <> void AverageByFour<GUInt32>(GUInt32 *buff, int xsz, int ysz)
132 : {
133 1 : GUInt32 *obuff = buff;
134 1 : GUInt32 *evenline = buff;
135 :
136 11 : for (int line = 0; line < ysz; line++)
137 : {
138 10 : GUInt32 *oddline = evenline + xsz * 2;
139 110 : for (int col = 0; col < xsz; col++)
140 : {
141 100 : *obuff++ = (GIntBig(2) + evenline[0] + evenline[1] + oddline[0] +
142 100 : oddline[1]) /
143 : 4;
144 100 : evenline += 2;
145 100 : oddline += 2;
146 : }
147 10 : evenline += xsz * 2; // Skips the other line
148 : }
149 1 : }
150 :
151 : // float specialization
152 1 : template <> void AverageByFour<float>(float *buff, int xsz, int ysz)
153 : {
154 1 : float *obuff = buff;
155 1 : float *evenline = buff;
156 :
157 11 : for (int line = 0; line < ysz; line++)
158 : {
159 10 : float *oddline = evenline + xsz * 2;
160 110 : for (int col = 0; col < xsz; col++)
161 : {
162 100 : *obuff++ =
163 100 : (evenline[0] + evenline[1] + oddline[0] + oddline[1]) * 0.25f;
164 100 : evenline += 2;
165 100 : oddline += 2;
166 : }
167 10 : evenline += xsz * 2; // Skips the other line
168 : }
169 1 : }
170 :
171 : // double specialization
172 1 : template <> void AverageByFour<double>(double *buff, int xsz, int ysz)
173 : {
174 1 : double *obuff = buff;
175 1 : double *evenline = buff;
176 :
177 11 : for (int line = 0; line < ysz; line++)
178 : {
179 10 : double *oddline = evenline + xsz * 2;
180 110 : for (int col = 0; col < xsz; col++)
181 : {
182 100 : *obuff++ =
183 100 : (evenline[0] + evenline[1] + oddline[0] + oddline[1]) * 0.25;
184 100 : evenline += 2;
185 100 : oddline += 2;
186 : }
187 10 : evenline += xsz * 2; // Skips the other line
188 : }
189 1 : }
190 :
191 : //
192 : // Integer type specialization, with roundup and integer math, avoids overflow
193 : // using GIntBig accumulator Speedup by specialization for smaller byte count
194 : // int types is probably not worth much since there are so many conditions here
195 : //
196 : template <typename T>
197 5 : static void AverageByFour(T *buff, int xsz, int ysz, T ndv)
198 : {
199 5 : T *obuff = buff;
200 5 : T *evenline = buff;
201 :
202 55 : for (int line = 0; line < ysz; line++)
203 : {
204 50 : T *oddline = evenline + xsz * 2;
205 550 : for (int col = 0; col < xsz; col++)
206 : {
207 500 : GIntBig acc = 0;
208 500 : int count = 0;
209 :
210 : // Temporary macro to accumulate the sum, uses the value, increments the pointer
211 : // Careful with this one, it has side effects
212 : #define use(valp) \
213 : if (*valp != ndv) \
214 : { \
215 : acc += *valp; \
216 : count++; \
217 : }; \
218 : valp++;
219 500 : use(evenline);
220 500 : use(evenline);
221 500 : use(oddline);
222 500 : use(oddline);
223 : #undef use
224 : // The count/2 is the bias to obtain correct rounding
225 500 : *obuff++ = T((count != 0) ? ((acc + count / 2) / count) : ndv);
226 : }
227 50 : evenline += xsz * 2; // Skips every other line
228 : }
229 5 : }
230 :
231 : // float specialization
232 1 : template <> void AverageByFour<float>(float *buff, int xsz, int ysz, float ndv)
233 : {
234 1 : float *obuff = buff;
235 1 : float *evenline = buff;
236 :
237 11 : for (int line = 0; line < ysz; line++)
238 : {
239 10 : float *oddline = evenline + xsz * 2;
240 110 : for (int col = 0; col < xsz; col++)
241 : {
242 100 : double acc = 0;
243 100 : double count = 0;
244 :
245 : // Temporary macro to accumulate the sum, uses the value, increments the pointer
246 : // Careful with this one, it has side effects
247 : #define use(valp) \
248 : if (*valp != ndv) \
249 : { \
250 : acc += *valp; \
251 : count += 1.0; \
252 : }; \
253 : valp++;
254 100 : use(evenline);
255 100 : use(evenline);
256 100 : use(oddline);
257 100 : use(oddline);
258 : #undef use
259 : // Output value is eiher accumulator divided by count or the
260 : // NoDataValue
261 100 : *obuff++ = float((count != 0.0) ? acc / count : ndv);
262 : }
263 10 : evenline += xsz * 2; // Skips every other line
264 : }
265 1 : }
266 :
267 : // double specialization, same as above
268 : template <>
269 1 : void AverageByFour<double>(double *buff, int xsz, int ysz, double ndv)
270 : {
271 1 : double *obuff = buff;
272 1 : double *evenline = buff;
273 :
274 11 : for (int line = 0; line < ysz; line++)
275 : {
276 10 : double *oddline = evenline + xsz * 2;
277 110 : for (int col = 0; col < xsz; col++)
278 : {
279 100 : double acc = 0;
280 100 : double count = 0;
281 :
282 : // Temporary macro to accumulate the sum, uses the value, increments the pointer
283 : // Careful with this one, it has side effects
284 : #define use(valp) \
285 : if (*valp != ndv) \
286 : { \
287 : acc += *valp; \
288 : count += 1.0; \
289 : }; \
290 : valp++;
291 100 : use(evenline);
292 100 : use(evenline);
293 100 : use(oddline);
294 100 : use(oddline);
295 : #undef use
296 : // Output value is eiher accumulator divided by count or the
297 : // NoDataValue
298 100 : *obuff++ = ((count != 0.0) ? acc / count : ndv);
299 : }
300 10 : evenline += xsz * 2; // Skips every other line
301 : }
302 1 : }
303 :
304 : /*
305 : *\brief Patches an overview for the selected area
306 : * arguments are in blocks in the source level, if toTheTop is false it only
307 : *does the next level It will read adjacent blocks if they are needed, so actual
308 : *area read might be padded by one block in either side
309 : */
310 :
311 26 : CPLErr MRFDataset::PatchOverview(int BlockX, int BlockY, int Width, int Height,
312 : int srcLevel, int recursive, int sampling_mode)
313 : {
314 26 : CPLErr status = CE_None;
315 26 : GDALRasterBand *b0 = GetRasterBand(1);
316 26 : if (b0->GetOverviewCount() <= srcLevel)
317 0 : return CE_None;
318 :
319 26 : int BlockXOut = BlockX / 2; // Round down
320 26 : Width += BlockX & 1; // Increment width if rounding down
321 26 : int BlockYOut = BlockY / 2; // Round down
322 26 : Height += BlockY & 1; // Increment height if rounding down
323 :
324 26 : int WidthOut = Width / 2 + (Width & 1); // Round up
325 26 : int HeightOut = Height / 2 + (Height & 1); // Round up
326 :
327 26 : int bands = GetRasterCount();
328 : int tsz_x, tsz_y;
329 26 : b0->GetBlockSize(&tsz_x, &tsz_y);
330 26 : GDALDataType eDataType = b0->GetRasterDataType();
331 :
332 : int pixel_size =
333 26 : GDALGetDataTypeSizeBytes(eDataType); // Bytes per pixel per band
334 26 : int line_size = tsz_x * pixel_size; // A line has this many bytes
335 26 : int buffer_size = line_size * tsz_y; // A block size in bytes
336 :
337 : // Build a vector of input and output bands
338 52 : std::vector<GDALRasterBand *> src_b;
339 52 : std::vector<GDALRasterBand *> dst_b;
340 :
341 52 : for (int band = 1; band <= bands; band++)
342 : {
343 26 : if (srcLevel == 0)
344 23 : src_b.push_back(GetRasterBand(band));
345 : else
346 3 : src_b.push_back(GetRasterBand(band)->GetOverview(srcLevel - 1));
347 26 : dst_b.push_back(GetRasterBand(band)->GetOverview(srcLevel));
348 : }
349 :
350 : // Allocate input space for four blocks
351 52 : std::vector<GByte> buffer(buffer_size * 4);
352 :
353 : // If the page is interleaved, we only need to check the page exists
354 : // otherwise we need to check each band block
355 26 : int check_bands = (bands == current.pagesize.c) ? 1 : bands;
356 :
357 : //
358 : // The inner loop is the band, so it is efficient for interleaved data.
359 : // There is no penalty for separate bands either.
360 : //
361 56 : for (int y = 0; y < HeightOut && CE_None == status; y++)
362 : {
363 30 : int dst_offset_y = BlockYOut + y;
364 30 : int src_offset_y = dst_offset_y * 2;
365 70 : for (int x = 0; x < WidthOut && CE_None == status; x++)
366 : {
367 40 : int dst_offset_x = BlockXOut + x;
368 40 : int src_offset_x = dst_offset_x * 2;
369 :
370 : // If none of the source blocks exists, there is no need to
371 : // read/write the blocks themselves
372 40 : bool has_data = false;
373 80 : for (int band = 0; band < check_bands; band++)
374 : {
375 : MRFRasterBand *bsrc =
376 40 : reinterpret_cast<MRFRasterBand *>(src_b[band]);
377 40 : has_data |= bsrc->TestBlock(src_offset_x, src_offset_y);
378 40 : has_data |= bsrc->TestBlock(src_offset_x + 1, src_offset_y);
379 40 : has_data |= bsrc->TestBlock(src_offset_x, src_offset_y + 1);
380 40 : has_data |= bsrc->TestBlock(src_offset_x + 1, src_offset_y + 1);
381 : }
382 :
383 : // No data in any of the bands for this output block
384 40 : if (!has_data)
385 : {
386 : // check that the output is already empty, otherwise force write
387 : // an empty block
388 0 : for (int band = 0; band < check_bands; band++)
389 : {
390 : MRFRasterBand *bdst =
391 0 : reinterpret_cast<MRFRasterBand *>(dst_b[band]);
392 0 : if (bdst->TestBlock(dst_offset_x, dst_offset_y))
393 : {
394 : // Output block exists, but it should not, force it
395 : ILSize req(dst_offset_x, dst_offset_y, 0, band,
396 0 : bdst->m_l);
397 0 : WriteTile(nullptr, IdxOffset(req, bdst->img));
398 : }
399 : }
400 : // No blocks in -> No block out
401 0 : continue;
402 : }
403 :
404 : // Do it band at a time so we can work in grayscale
405 80 : for (int band = 0; band < bands; band++)
406 : { // Counting from zero in a vector
407 :
408 40 : int sz_x = 2 * tsz_x, sz_y = 2 * tsz_y;
409 40 : MRFRasterBand *bsrc = static_cast<MRFRasterBand *>(src_b[band]);
410 40 : MRFRasterBand *bdst = static_cast<MRFRasterBand *>(dst_b[band]);
411 :
412 : //
413 : // Clip to the size to the input image
414 : // This is one of the worst features of GDAL, it doesn't
415 : // tolerate any padding
416 : //
417 40 : bool adjusted = false;
418 40 : if (bsrc->GetXSize() < (src_offset_x + 2) * tsz_x)
419 : {
420 9 : sz_x = bsrc->GetXSize() - src_offset_x * tsz_x;
421 9 : adjusted = true;
422 : }
423 40 : if (bsrc->GetYSize() < (src_offset_y + 2) * tsz_y)
424 : {
425 9 : sz_y = bsrc->GetYSize() - src_offset_y * tsz_y;
426 9 : adjusted = true;
427 : }
428 :
429 40 : if (adjusted)
430 : { // Fill with no data for partial buffer, instead of padding
431 : // afterwards
432 13 : size_t bsb = bsrc->blockSizeBytes();
433 13 : auto b = buffer.data();
434 13 : bsrc->FillBlock(b);
435 13 : bsrc->FillBlock(b + bsb);
436 13 : bsrc->FillBlock(b + 2 * bsb);
437 13 : bsrc->FillBlock(b + 3 * bsb);
438 : }
439 :
440 40 : int hasNoData = 0;
441 40 : double ndv = bsrc->GetNoDataValue(&hasNoData);
442 :
443 40 : status = bsrc->RasterIO(
444 : GF_Read, src_offset_x * tsz_x,
445 : src_offset_y * tsz_y, // offset in input image
446 : sz_x, sz_y, // Size in output image
447 40 : buffer.data(), sz_x, sz_y, // Buffer and size in buffer
448 40 : eDataType, pixel_size, 2 * line_size, nullptr);
449 :
450 40 : if (CE_None != status)
451 : {
452 0 : CPLError(CE_Failure, CPLE_AppDefined,
453 : "MRF: Patch - RasterIO() read failed");
454 0 : break; // Get out now
455 : }
456 :
457 : // Count the NoData values
458 40 : int count = 0; // Assume all points are data
459 40 : if (sampling_mode == SAMPLING_Avg)
460 : {
461 :
462 : // Dispatch based on data type
463 : // Use a temporary macro to make it look easy
464 : // Runs the optimized version if the page is full with data
465 : #define resample(T) \
466 : if (hasNoData) \
467 : { \
468 : count = MatchCount((T *)buffer.data(), 4 * tsz_x * tsz_y, T(ndv)); \
469 : if (4 * tsz_x * tsz_y == count) \
470 : bdst->FillBlock(buffer.data()); \
471 : else if (0 != count) \
472 : AverageByFour((T *)buffer.data(), tsz_x, tsz_y, T(ndv)); \
473 : } \
474 : if (0 == count) \
475 : AverageByFour((T *)buffer.data(), tsz_x, tsz_y); \
476 : break;
477 :
478 14 : switch (eDataType)
479 : {
480 2 : case GDT_Byte:
481 2 : resample(GByte);
482 2 : case GDT_UInt16:
483 2 : resample(GUInt16);
484 2 : case GDT_Int16:
485 2 : resample(GInt16);
486 2 : case GDT_UInt32:
487 2 : resample(GUInt32);
488 2 : case GDT_Int32:
489 2 : resample(GInt32);
490 2 : case GDT_Float32:
491 2 : resample(float);
492 2 : case GDT_Float64:
493 2 : resample(double);
494 0 : default:
495 0 : CPLAssert(false);
496 : break;
497 : }
498 : #undef resample
499 : }
500 26 : else if (sampling_mode == SAMPLING_Near)
501 : {
502 :
503 : #define resample(T) \
504 : if (hasNoData) \
505 : { \
506 : count = MatchCount((T *)buffer.data(), 4 * tsz_x * tsz_y, T(ndv)); \
507 : if (4 * tsz_x * tsz_y == count) \
508 : bdst->FillBlock(buffer.data()); \
509 : else if (0 != count) \
510 : NearByFour((T *)buffer.data(), tsz_x, tsz_y, T(ndv)); \
511 : } \
512 : if (0 == count) \
513 : NearByFour((T *)buffer.data(), tsz_x, tsz_y); \
514 : break;
515 26 : switch (eDataType)
516 : {
517 20 : case GDT_Byte:
518 20 : resample(GByte);
519 1 : case GDT_UInt16:
520 1 : resample(GUInt16);
521 1 : case GDT_Int16:
522 1 : resample(GInt16);
523 1 : case GDT_UInt32:
524 1 : resample(GUInt32);
525 1 : case GDT_Int32:
526 1 : resample(GInt32);
527 1 : case GDT_Float32:
528 1 : resample(float);
529 1 : case GDT_Float64:
530 1 : resample(double);
531 0 : default:
532 0 : CPLAssert(false);
533 : break;
534 : }
535 : #undef resample
536 : }
537 :
538 : // Done filling the buffer
539 : // Argh, still need to clip the output to the band size on the
540 : // right and bottom The offset should be fine, just the size
541 : // might need adjustments
542 40 : sz_x = tsz_x;
543 40 : sz_y = tsz_y;
544 :
545 40 : if (bdst->GetXSize() < dst_offset_x * sz_x + sz_x)
546 9 : sz_x = bdst->GetXSize() - dst_offset_x * sz_x;
547 40 : if (bdst->GetYSize() < dst_offset_y * sz_y + sz_y)
548 9 : sz_y = bdst->GetYSize() - dst_offset_y * sz_y;
549 :
550 40 : status = bdst->RasterIO(
551 : GF_Write, dst_offset_x * tsz_x,
552 : dst_offset_y * tsz_y, // offset in output image
553 : sz_x, sz_y, // Size in output image
554 40 : buffer.data(), sz_x, sz_y, // Buffer and size in buffer
555 : eDataType, pixel_size, line_size, nullptr);
556 :
557 40 : if (CE_None != status)
558 : {
559 0 : CPLError(CE_Failure, CPLE_AppDefined,
560 : "MRF: Patch - RasterIO() write failed");
561 0 : break;
562 : }
563 : } // Band loop
564 :
565 : // Mark input data as no longer needed, saves RAM
566 80 : for (int band = 0; band < bands; band++)
567 40 : src_b[band]->FlushCache(false);
568 : }
569 : }
570 :
571 26 : if (CE_None != status)
572 0 : return status; // Report problems
573 :
574 52 : for (int band = 0; band < bands; band++)
575 26 : dst_b[band]->FlushCache(
576 26 : false); // Commit destination to disk after each overview
577 :
578 26 : if (!recursive)
579 26 : return CE_None;
580 0 : return PatchOverview(BlockXOut, BlockYOut, WidthOut, HeightOut,
581 0 : srcLevel + 1, true);
582 : }
583 :
584 : NAMESPACE_MRF_END
|
