Line data Source code
1 : /******************************************************************************
2 : *
3 : * Project: NITF Read/Write Library
4 : * Purpose: Creates RPFHDR, RPFIMG, RPFDES TREs and RPF image data
5 : * Author: Even Rouault, even dot rouault at spatialys dot com
6 : *
7 : **********************************************************************
8 : * Copyright (c) 2026, T-Kartor
9 : *
10 : * SPDX-License-Identifier: MIT
11 : ****************************************************************************/
12 :
13 : #include "rpfframewriter.h"
14 :
15 : #include "cpl_enumerate.h"
16 : #include "cpl_string.h"
17 : #include "cpl_time.h"
18 : #include "cpl_worker_thread_pool.h"
19 : #include "gdal_alg.h"
20 : #include "gdal_alg_priv.h"
21 : #include "gdal_colortable.h"
22 : #include "gdal_dataset.h"
23 : #include "gdal_geotransform.h"
24 : #include "gdal_rasterband.h"
25 : #include "gdal_thread_pool.h"
26 : #include "gdal_utils.h"
27 : #include "ogr_spatialref.h"
28 : #include "offsetpatcher.h"
29 : #include "nitfdataset.h"
30 : #include "nitflib.h"
31 : #include "kdtree_vq_cadrg.h"
32 : #include "rpftocwriter.h"
33 :
34 : #include <algorithm>
35 : #include <array>
36 : #include <map>
37 : #include <mutex>
38 : #include <utility>
39 : #include <vector>
40 :
41 : // Refer to following documents to (hopefully) understand this file:
42 : // MIL-C-89038, CADRG: http://everyspec.com/MIL-PRF/MIL-PRF-080000-99999/MIL-PRF-89038_25371/
43 : // MIL-STD-2411, RPF: http://everyspec.com/MIL-STD/MIL-STD-2000-2999/MIL-STD-2411_6903/
44 : // MIL-STD-2411-1, RPF details: http://everyspec.com/MIL-STD/MIL-STD-2000-2999/MIL-STD-2411-1_6909/
45 : // MIL-STD-2411-2, RPF-in-NITF: http://everyspec.com/MIL-STD/MIL-STD-2000-2999/MIL-STD-2411-2_6908/
46 : // MIL-A-89007, ADRG: http://everyspec.com/MIL-SPECS/MIL-SPECS-MIL-A/MIL-A-89007_51725/
47 : // MIL-STD-188/199, VQ NITF: http://everyspec.com/MIL-STD/MIL-STD-0100-0299/MIL_STD_188_199_1730/
48 :
49 : constexpr int SUBSAMPLING = 4;
50 : constexpr int BLOCK_SIZE = 256;
51 : constexpr int CODEBOOK_MAX_SIZE = 4096;
52 : constexpr int TRANSPARENT_CODEBOOK_CODE = CODEBOOK_MAX_SIZE - 1;
53 : constexpr int TRANSPARENT_COLOR_TABLE_ENTRY = CADRG_MAX_COLOR_ENTRY_COUNT;
54 :
55 : constexpr int CADRG_SECOND_CT_COUNT = 32;
56 : constexpr int CADRG_THIRD_CT_COUNT = 16;
57 :
58 : constexpr double CADRG_PITCH_IN_CM = 0.0150; // 150 micrometers
59 :
60 : constexpr int DEFAULT_DENSIFY_PTS = 21;
61 :
62 : /************************************************************************/
63 : /* RPFCADRGIsValidZone() */
64 : /************************************************************************/
65 :
66 1549 : static bool RPFCADRGIsValidZone(int nZone)
67 : {
68 1549 : return nZone >= MIN_ZONE && nZone <= MAX_ZONE;
69 : }
70 :
71 : /************************************************************************/
72 : /* RPFCADRGZoneNumToChar() */
73 : /************************************************************************/
74 :
75 109 : char RPFCADRGZoneNumToChar(int nZone)
76 : {
77 109 : CPLAssert(RPFCADRGIsValidZone(nZone));
78 109 : if (nZone <= MAX_ZONE_NORTHERN_HEMISPHERE)
79 102 : return static_cast<char>('0' + nZone);
80 7 : else if (nZone < MAX_ZONE)
81 4 : return static_cast<char>('A' + (nZone - MIN_ZONE_SOUTHERN_HEMISPHERE));
82 : else
83 3 : return 'J';
84 : }
85 :
86 : /************************************************************************/
87 : /* RPFCADRGZoneCharToNum() */
88 : /************************************************************************/
89 :
90 70 : int RPFCADRGZoneCharToNum(char chZone)
91 : {
92 70 : if (chZone >= '1' && chZone <= '9')
93 59 : return chZone - '0';
94 11 : else if (chZone >= 'A' && chZone <= 'H')
95 5 : return (chZone - 'A') + MIN_ZONE_SOUTHERN_HEMISPHERE;
96 6 : else if (chZone >= 'a' && chZone <= 'h')
97 3 : return (chZone - 'a') + MIN_ZONE_SOUTHERN_HEMISPHERE;
98 3 : else if (chZone == 'J' || chZone == 'j')
99 3 : return MAX_ZONE;
100 : else
101 0 : return 0;
102 : }
103 :
104 : /************************************************************************/
105 : /* RPFCADRGGetScaleFromDataSeriesCode() */
106 : /************************************************************************/
107 :
108 70 : int RPFCADRGGetScaleFromDataSeriesCode(const char *pszCode)
109 : {
110 70 : int nVal = 0;
111 70 : const auto *psEntry = NITFGetRPFSeriesInfoFromCode(pszCode);
112 70 : if (psEntry)
113 : {
114 70 : nVal = NITFGetScaleFromScaleResolution(psEntry->scaleResolution);
115 : }
116 70 : return nVal;
117 : }
118 :
119 : /************************************************************************/
120 : /* RPFCADRGIsKnownDataSeriesCode() */
121 : /************************************************************************/
122 :
123 66 : bool RPFCADRGIsKnownDataSeriesCode(const char *pszCode)
124 : {
125 66 : return NITFIsKnownRPFDataSeriesCode(pszCode, "CADRG");
126 : }
127 :
128 : /************************************************************************/
129 : /* CADRGInformation::Private */
130 : /************************************************************************/
131 :
132 : class CADRGInformation::Private
133 : {
134 : public:
135 : std::vector<BucketItem<ColorTableBased4x4Pixels>> codebook{};
136 : std::vector<short> VQImage{};
137 : bool bHasTransparentPixels = false;
138 : };
139 :
140 : /************************************************************************/
141 : /* CADRGInformation::CADRGInformation() */
142 : /************************************************************************/
143 :
144 33 : CADRGInformation::CADRGInformation(std::unique_ptr<Private> priv)
145 33 : : m_private(std::move(priv))
146 : {
147 33 : }
148 :
149 : /************************************************************************/
150 : /* CADRGInformation::~CADRGInformation() */
151 : /************************************************************************/
152 :
153 : CADRGInformation::~CADRGInformation() = default;
154 :
155 : /************************************************************************/
156 : /* CADRGInformation::HasTransparentPixels() */
157 : /************************************************************************/
158 :
159 65 : bool CADRGInformation::HasTransparentPixels() const
160 : {
161 65 : return m_private->bHasTransparentPixels;
162 : }
163 :
164 : /************************************************************************/
165 : /* StrPadTruncate() */
166 : /************************************************************************/
167 :
168 : #ifndef StrPadTruncate_defined
169 : #define StrPadTruncate_defined
170 :
171 289 : static std::string StrPadTruncate(const std::string &osIn, size_t nSize)
172 : {
173 289 : std::string osOut(osIn);
174 289 : osOut.resize(nSize, ' ');
175 289 : return osOut;
176 : }
177 : #endif
178 :
179 : /************************************************************************/
180 : /* Create_CADRG_RPFHDR() */
181 : /************************************************************************/
182 :
183 53 : void Create_CADRG_RPFHDR(GDALOffsetPatcher::OffsetPatcher *offsetPatcher,
184 : const std::string &osFilename,
185 : CPLStringList &aosOptions)
186 : {
187 53 : auto poRPFHDR = offsetPatcher->CreateBuffer(
188 : "RPFHDR", /* bEndiannessIsLittle = */ false);
189 53 : CPLAssert(poRPFHDR);
190 : #ifdef INCLUDE_HEADER_AND_LOCATION
191 : poRPFHDR->DeclareOffsetAtCurrentPosition("HEADER_COMPONENT_LOCATION");
192 : #endif
193 53 : poRPFHDR->AppendByte(0); // big endian order
194 53 : poRPFHDR->AppendUInt16RefForSizeOfBuffer("RPFHDR");
195 53 : poRPFHDR->AppendString(
196 106 : StrPadTruncate(CPLGetFilename(osFilename.c_str()), 12));
197 53 : poRPFHDR->AppendByte(0); // update indicator: initial release
198 53 : poRPFHDR->AppendString("MIL-C-89038 "); // GOVERNING_STANDARD_NUMBER
199 53 : poRPFHDR->AppendString("19941006"); // GOVERNING_STANDARD_DATE
200 : // SECURITY_CLASSIFICATION
201 53 : poRPFHDR->AppendString(
202 106 : StrPadTruncate(aosOptions.FetchNameValueDef("FCLASS", "U"), 1));
203 53 : poRPFHDR->AppendString(StrPadTruncate(
204 : aosOptions.FetchNameValueDef("SECURITY_COUNTRY_CODE", " "), 2));
205 53 : poRPFHDR->AppendString(" "); // SECURITY_RELEASE_MARKING
206 53 : poRPFHDR->AppendUInt32RefForOffset("LOCATION_COMPONENT_LOCATION");
207 :
208 53 : char *pszEscaped = CPLEscapeString(
209 53 : reinterpret_cast<const char *>(poRPFHDR->GetBuffer().data()),
210 53 : static_cast<int>(poRPFHDR->GetBuffer().size()),
211 : CPLES_BackslashQuotable);
212 : aosOptions.AddString(
213 53 : std::string("FILE_TRE=RPFHDR=").append(pszEscaped).c_str());
214 53 : CPLFree(pszEscaped);
215 53 : }
216 :
217 : /************************************************************************/
218 : /* Create_CADRG_LocationComponent() */
219 : /************************************************************************/
220 :
221 : static void
222 33 : Create_CADRG_LocationComponent(GDALOffsetPatcher::OffsetPatcher *offsetPatcher)
223 : {
224 33 : auto poBuffer = offsetPatcher->CreateBuffer(
225 : "LocationComponent", /* bEndiannessIsLittle = */ false);
226 33 : CPLAssert(poBuffer);
227 33 : poBuffer->DeclareOffsetAtCurrentPosition("LOCATION_COMPONENT_LOCATION");
228 :
229 : static const struct
230 : {
231 : uint16_t locationId;
232 : const char *locationBufferName;
233 : const char *locationOffsetName;
234 : } asLocations[] = {
235 : #ifdef INCLUDE_HEADER_AND_LOCATION
236 : // While it shouldn't hurt, it doesn't seem idiomatical to include
237 : // those locations in the location table.
238 : {LID_HeaderComponent /* 128 */, "RPFHDR", "HEADER_COMPONENT_LOCATION"},
239 : {LID_LocationComponent /* 129 */, "LocationComponent",
240 : "LOCATION_COMPONENT_LOCATION"},
241 : #endif
242 : {LID_CoverageSectionSubheader /* 130 */, "CoverageSectionSubheader",
243 : "COVERAGE_SECTION_LOCATION"},
244 : {LID_CompressionSectionSubsection /* 131 */, "CompressionSection",
245 : "COMPRESSION_SECTION_LOCATION"},
246 : {LID_CompressionLookupSubsection /* 132 */,
247 : "CompressionLookupSubsection", "COMPRESSION_LOOKUP_LOCATION"},
248 : /* no LID_CompressionParameterSubsection = 133 in CADRG */
249 : {LID_ColorGrayscaleSectionSubheader /* 134 */,
250 : "ColorGrayscaleSectionSubheader", "COLOR_GRAYSCALE_LOCATION"},
251 : {LID_ColormapSubsection /* 135 */, "ColormapSubsection",
252 : "COLORMAP_LOCATION"},
253 : {LID_ImageDescriptionSubheader /* 136 */, "ImageDescriptionSubheader",
254 : "IMAGE_DESCRIPTION_SECTION_LOCATION"},
255 : {LID_ImageDisplayParametersSubheader /* 137 */,
256 : "ImageDisplayParametersSubheader",
257 : "IMAGE_DISPLAY_PARAMETERS_SECTION_LOCATION"},
258 : {LID_MaskSubsection /* 138 */, "MaskSubsection",
259 : "MASK_SUBSECTION_LOCATION"},
260 : {LID_ColorConverterSubsection /* 139 */, "ColorConverterSubsection",
261 : "COLOR_CONVERTER_SUBSECTION"},
262 : {LID_SpatialDataSubsection /* 140 */, "SpatialDataSubsection",
263 : "SPATIAL_DATA_SUBSECTION_LOCATION"},
264 : {LID_AttributeSectionSubheader /* 141 */, "AttributeSectionSubheader",
265 : "ATTRIBUTE_SECTION_SUBHEADER_LOCATION"},
266 : {LID_AttributeSubsection /* 142 */, "AttributeSubsection",
267 : "ATTRIBUTE_SUBSECTION_LOCATION"},
268 : };
269 :
270 66 : std::string sumOfSizes;
271 33 : uint16_t nComponents = 0;
272 429 : for (const auto &sLocation : asLocations)
273 : {
274 396 : ++nComponents;
275 396 : if (!sumOfSizes.empty())
276 363 : sumOfSizes += '+';
277 396 : sumOfSizes += sLocation.locationBufferName;
278 : }
279 :
280 33 : constexpr uint16_t COMPONENT_LOCATION_OFFSET = 14;
281 33 : constexpr uint16_t COMPONENT_LOCATION_RECORD_LENGTH = 10;
282 33 : poBuffer->AppendUInt16RefForSizeOfBuffer("LocationComponent");
283 33 : poBuffer->AppendUInt32(COMPONENT_LOCATION_OFFSET);
284 33 : poBuffer->AppendUInt16(nComponents);
285 33 : poBuffer->AppendUInt16(COMPONENT_LOCATION_RECORD_LENGTH);
286 : // COMPONENT_AGGREGATE_LENGTH
287 33 : poBuffer->AppendUInt32RefForSizeOfBuffer(sumOfSizes);
288 :
289 429 : for (const auto &sLocation : asLocations)
290 : {
291 396 : poBuffer->AppendUInt16(sLocation.locationId);
292 396 : poBuffer->AppendUInt32RefForSizeOfBuffer(sLocation.locationBufferName);
293 396 : poBuffer->AppendUInt32RefForOffset(sLocation.locationOffsetName);
294 : }
295 33 : }
296 :
297 : /************************************************************************/
298 : /* asARCZoneDefinitions */
299 : /************************************************************************/
300 :
301 : constexpr double ARC_B = 400384;
302 :
303 : // Content of MIL-A-89007 (ADRG specification), appendix 70, table III
304 : static constexpr struct
305 : {
306 : int nZone; // zone number (for northern hemisphere. Add 9 for southern hemisphere)
307 : int minLat; // minimum latitude of the zone
308 : int maxLat; // maximum latitude of the zone
309 : double A; // longitudinal pixel spacing constant at 1:1M
310 : double B; // latitudinal pixel spacing constant at 1:1M
311 : double latRes; // in microns
312 : double lonRes; // in microns
313 : } asARCZoneDefinitions[] = {
314 : {1, 0, 32, 369664, ARC_B, 99.9, 99.9},
315 : {2, 32, 48, 302592, ARC_B, 99.9, 99.9},
316 : {3, 48, 56, 245760, ARC_B, 100.0, 99.9},
317 : {4, 56, 64, 199168, ARC_B, 99.9, 99.9},
318 : {5, 64, 68, 163328, ARC_B, 99.7, 99.9},
319 : {6, 68, 72, 137216, ARC_B, 99.7, 99.9},
320 : {7, 72, 76, 110080, ARC_B, 99.8, 99.9},
321 : {8, 76, 80, 82432, ARC_B, 100.0, 99.9},
322 : {9, 80, 90, ARC_B, ARC_B, 99.9, 99.9},
323 : };
324 :
325 : constexpr double RATIO_PITCH_CADRG_OVER_ADRG = 150.0 / 100.0;
326 : constexpr double REF_SCALE = 1e6;
327 : constexpr int ADRG_BLOCK_SIZE = 512;
328 :
329 : /************************************************************************/
330 : /* GetARCZoneFromLat() */
331 : /************************************************************************/
332 :
333 0 : static int GetARCZoneFromLat(double dfLat)
334 : {
335 0 : for (const auto &sZoneDef : asARCZoneDefinitions)
336 : {
337 0 : if (std::fabs(dfLat) >= sZoneDef.minLat &&
338 0 : std::fabs(dfLat) <= sZoneDef.maxLat)
339 : {
340 0 : return dfLat >= 0 ? sZoneDef.nZone
341 0 : : sZoneDef.nZone + MAX_ZONE_NORTHERN_HEMISPHERE;
342 : }
343 : }
344 0 : return 0;
345 : }
346 :
347 : /************************************************************************/
348 : /* GetPolarConstant() */
349 : /************************************************************************/
350 :
351 70 : static double GetPolarConstant(int nReciprocalScale)
352 : {
353 : // Cf MIL-A-89007 (ADRG specification), appendix 70, table III
354 70 : const double N = REF_SCALE / nReciprocalScale;
355 :
356 : // Cf MIL-C-89038 (CADRG specification), para 60.4
357 70 : const double B_s = ARC_B * N;
358 70 : const double latCst_ADRG =
359 70 : std::ceil(B_s / ADRG_BLOCK_SIZE) * ADRG_BLOCK_SIZE;
360 70 : return std::round(latCst_ADRG * 20.0 / 360.0 / RATIO_PITCH_CADRG_OVER_ADRG /
361 70 : ADRG_BLOCK_SIZE) *
362 70 : ADRG_BLOCK_SIZE * 360 / 20;
363 : }
364 :
365 : /************************************************************************/
366 : /* GetYPixelSize() */
367 : /************************************************************************/
368 :
369 : /** Return the size of a pixel (in degree for non-polar zones, in meters for
370 : * polar zones), along the latitude/Y axis,
371 : * at specified scale and zone */
372 586 : static double GetYPixelSize(int nZone, int nReciprocalScale)
373 : {
374 586 : CPLAssert(RPFCADRGIsValidZone(nZone));
375 586 : const int nZoneIdx = (nZone - 1) % MAX_ZONE_NORTHERN_HEMISPHERE;
376 :
377 586 : if (nZoneIdx + 1 == MAX_ZONE_NORTHERN_HEMISPHERE)
378 : {
379 18 : const double polarCst = GetPolarConstant(nReciprocalScale);
380 18 : return SRS_WGS84_SEMIMAJOR * 2 * M_PI / polarCst;
381 : }
382 :
383 568 : const auto &sZoneDef = asARCZoneDefinitions[nZoneIdx];
384 :
385 : // Cf MIL-A-89007 (ADRG specification), appendix 70, table III
386 568 : const double N = REF_SCALE / nReciprocalScale;
387 :
388 : // Cf MIL-C-89038 (CADRG specification), para 60.1.1 and following
389 568 : const double B_s = sZoneDef.B * N;
390 568 : const double latCst_ADRG =
391 568 : std::ceil(B_s / ADRG_BLOCK_SIZE) * ADRG_BLOCK_SIZE;
392 568 : const double latCst_CADRG =
393 568 : std::round(latCst_ADRG / RATIO_PITCH_CADRG_OVER_ADRG / 4 / BLOCK_SIZE) *
394 : BLOCK_SIZE;
395 568 : const double latInterval = 90.0 / latCst_CADRG;
396 :
397 568 : return latInterval;
398 : }
399 :
400 : /************************************************************************/
401 : /* GetXPixelSize() */
402 : /************************************************************************/
403 :
404 : /** Return the size of a pixel (in degree for non-polar zones, in meters for
405 : * polar zones), along the longitude/X axis,
406 : * at specified scale and zone */
407 299 : static double GetXPixelSize(int nZone, int nReciprocalScale)
408 : {
409 299 : CPLAssert(RPFCADRGIsValidZone(nZone));
410 299 : const int nZoneIdx = (nZone - MIN_ZONE) % MAX_ZONE_NORTHERN_HEMISPHERE;
411 :
412 299 : if (nZoneIdx + 1 == MAX_ZONE_NORTHERN_HEMISPHERE)
413 : {
414 18 : const double polarCst = GetPolarConstant(nReciprocalScale);
415 18 : return SRS_WGS84_SEMIMAJOR * 2 * M_PI / polarCst;
416 : }
417 :
418 281 : const auto &sZoneDef = asARCZoneDefinitions[nZoneIdx];
419 :
420 : // Cf MIL-A-89007 (ADRG specification), appendix 70, table III
421 281 : const double N = REF_SCALE / nReciprocalScale;
422 :
423 : // Cf MIL-C-89038 (CADRG specification), para 60.1.1 and following
424 281 : const double A_s = sZoneDef.A * N;
425 281 : const double lonCst_ADRG =
426 281 : std::ceil(A_s / ADRG_BLOCK_SIZE) * ADRG_BLOCK_SIZE;
427 281 : const double lonCst_CADRG =
428 281 : std::round(lonCst_ADRG / RATIO_PITCH_CADRG_OVER_ADRG / BLOCK_SIZE) *
429 : BLOCK_SIZE;
430 281 : const double lonInterval = 360.0 / lonCst_CADRG;
431 :
432 281 : return lonInterval;
433 : }
434 :
435 : /************************************************************************/
436 : /* RPFGetCADRGResolutionAndInterval() */
437 : /************************************************************************/
438 :
439 56 : void RPFGetCADRGResolutionAndInterval(int nZone, int nReciprocalScale,
440 : double &latResolution,
441 : double &lonResolution,
442 : double &latInterval, double &lonInterval)
443 : {
444 56 : CPLAssert(nReciprocalScale > 0);
445 56 : CPLAssert(RPFCADRGIsValidZone(nZone));
446 56 : const int nZoneIdx = (nZone - MIN_ZONE) % MAX_ZONE_NORTHERN_HEMISPHERE;
447 56 : const auto &sZoneDef = asARCZoneDefinitions[nZoneIdx];
448 :
449 : // Cf MIL-A-89007 (ADRG specification), appendix 70, table III
450 56 : const double N = REF_SCALE / std::max(1, nReciprocalScale);
451 :
452 : // Cf MIL-C-89038 (CADRG specification), para 60.1.1 and following
453 56 : const double B_s = sZoneDef.B * N;
454 56 : const double latCst_ADRG =
455 56 : std::ceil(B_s / ADRG_BLOCK_SIZE) * ADRG_BLOCK_SIZE;
456 56 : const double latCst_CADRG =
457 56 : std::round(latCst_ADRG / RATIO_PITCH_CADRG_OVER_ADRG / 4 / BLOCK_SIZE) *
458 : BLOCK_SIZE;
459 56 : double latResolutionLocal =
460 56 : sZoneDef.latRes / N * latCst_ADRG / (4 * latCst_CADRG);
461 :
462 56 : const double A_s = sZoneDef.A * N;
463 56 : const double lonCst_ADRG =
464 56 : std::ceil(A_s / ADRG_BLOCK_SIZE) * ADRG_BLOCK_SIZE;
465 56 : const double lonCst_CADRG =
466 56 : std::round(lonCst_ADRG / RATIO_PITCH_CADRG_OVER_ADRG / BLOCK_SIZE) *
467 : BLOCK_SIZE;
468 56 : double lonResolutionLocal =
469 56 : sZoneDef.lonRes / N * lonCst_ADRG / lonCst_CADRG;
470 :
471 56 : double latIntervalLocal = 90.0 / latCst_CADRG;
472 56 : double lonIntervalLocal = 360.0 / lonCst_CADRG;
473 :
474 56 : if (nZoneIdx + MIN_ZONE == MAX_ZONE_NORTHERN_HEMISPHERE)
475 : {
476 6 : lonResolutionLocal = latResolutionLocal;
477 6 : lonIntervalLocal = latIntervalLocal;
478 : }
479 :
480 56 : latResolution = latResolutionLocal;
481 56 : lonResolution = lonResolutionLocal;
482 56 : latInterval = latIntervalLocal;
483 56 : lonInterval = lonIntervalLocal;
484 56 : }
485 :
486 : /************************************************************************/
487 : /* GetMinMaxLatWithOverlap() */
488 : /************************************************************************/
489 :
490 : /** Return the actual minimum and maximum latitude of a zone, for a given
491 : * reciprocal scale, taking into potential overlap between zones.
492 : */
493 287 : static std::pair<double, double> GetMinMaxLatWithOverlap(int nZone,
494 : int nReciprocalScale)
495 : {
496 287 : if (nZone == MAX_ZONE_NORTHERN_HEMISPHERE)
497 : {
498 0 : return std::pair(80.0, 90.0);
499 : }
500 287 : else if (nZone == MAX_ZONE)
501 : {
502 0 : return std::pair(-90.0, -80.0);
503 : }
504 287 : CPLAssert(RPFCADRGIsValidZone(nZone));
505 287 : const int nZoneIdx = (nZone - MIN_ZONE) % MAX_ZONE_NORTHERN_HEMISPHERE;
506 287 : const auto &sZoneDef = asARCZoneDefinitions[nZoneIdx];
507 :
508 287 : const double latInterval = GetYPixelSize(nZone, nReciprocalScale);
509 287 : const double deltaLatFrame = latInterval * CADRG_FRAME_PIXEL_COUNT;
510 :
511 287 : const double dfMinLat =
512 287 : std::floor(sZoneDef.minLat / deltaLatFrame) * deltaLatFrame;
513 287 : const double dfMaxLat =
514 287 : std::ceil(sZoneDef.maxLat / deltaLatFrame) * deltaLatFrame;
515 : return nZone >= MIN_ZONE_SOUTHERN_HEMISPHERE
516 287 : ? std::pair(-dfMaxLat, -dfMinLat)
517 287 : : std::pair(dfMinLat, dfMaxLat);
518 : }
519 :
520 : /************************************************************************/
521 : /* GetPolarFrameCount() */
522 : /************************************************************************/
523 :
524 : constexpr double EPSILON_1Em3 = 1e-3;
525 :
526 34 : static int GetPolarFrameCount(int nReciprocalScale)
527 : {
528 34 : const double numberSubFrames = std::round(
529 34 : GetPolarConstant(nReciprocalScale) * 20.0 / 360.0 / BLOCK_SIZE);
530 34 : constexpr int SUBFRAMES_PER_FRAME = CADRG_FRAME_PIXEL_COUNT / BLOCK_SIZE;
531 34 : int numberFrames = static_cast<int>(
532 34 : std::ceil(numberSubFrames / SUBFRAMES_PER_FRAME - EPSILON_1Em3));
533 34 : if ((numberFrames % 2) == 0)
534 0 : ++numberFrames;
535 34 : return numberFrames;
536 : }
537 :
538 : /************************************************************************/
539 : /* GetFrameCountAlongX() */
540 : /************************************************************************/
541 :
542 : constexpr double dfMinLonZone = -180.0;
543 : constexpr double dfMaxLonZone = 180.0;
544 :
545 96 : static int GetFrameCountAlongX(int nZone, int nReciprocalScale)
546 : {
547 96 : if (nZone == MAX_ZONE_NORTHERN_HEMISPHERE || nZone == MAX_ZONE)
548 8 : return GetPolarFrameCount(nReciprocalScale);
549 88 : const double lonInterval = GetXPixelSize(nZone, nReciprocalScale);
550 88 : const double eastWestPixelCst = 360.0 / lonInterval;
551 88 : CPLDebugOnly("CADRG", "eastWestPixelCst=%f, count=%f", eastWestPixelCst,
552 : eastWestPixelCst / CADRG_FRAME_PIXEL_COUNT);
553 88 : const int nFrameCount = static_cast<int>(
554 88 : std::ceil(eastWestPixelCst / CADRG_FRAME_PIXEL_COUNT - EPSILON_1Em3));
555 88 : return nFrameCount;
556 : }
557 :
558 : /************************************************************************/
559 : /* GetFrameCountAlongY() */
560 : /************************************************************************/
561 :
562 96 : static int GetFrameCountAlongY(int nZone, int nReciprocalScale)
563 : {
564 96 : if (nZone == MAX_ZONE_NORTHERN_HEMISPHERE || nZone == MAX_ZONE)
565 8 : return GetPolarFrameCount(nReciprocalScale);
566 88 : const double latInterval = GetYPixelSize(nZone, nReciprocalScale);
567 88 : const double deltaLatFrame = latInterval * CADRG_FRAME_PIXEL_COUNT;
568 88 : const auto [dfMinLatZone, dfMaxLatZone] =
569 88 : GetMinMaxLatWithOverlap(nZone, nReciprocalScale);
570 176 : return std::max(1, static_cast<int>(std::ceil(dfMaxLatZone - dfMinLatZone) /
571 88 : deltaLatFrame -
572 88 : EPSILON_1Em3));
573 : }
574 :
575 : /************************************************************************/
576 : /* RPFGetCADRGFramesForEnvelope() */
577 : /************************************************************************/
578 :
579 : enum class HemiphereType
580 : {
581 : BOTH,
582 : NORTH,
583 : SOUTH,
584 : };
585 :
586 : static std::vector<RPFFrameDef>
587 56 : RPFGetCADRGFramesForEnvelope(int nZoneIn, int nReciprocalScale,
588 : GDALDataset *poSrcDS, HemiphereType hemisphere)
589 : {
590 56 : CPLAssert(nZoneIn == 0 || RPFCADRGIsValidZone(nZoneIn));
591 :
592 56 : OGREnvelope sExtentNativeCRS;
593 56 : OGREnvelope sExtentWGS84;
594 112 : if (poSrcDS->GetExtent(&sExtentNativeCRS) != CE_None ||
595 56 : poSrcDS->GetExtentWGS84LongLat(&sExtentWGS84) != CE_None)
596 : {
597 0 : CPLError(CE_Failure, CPLE_AppDefined, "Cannot get dataset extent");
598 0 : return {};
599 : }
600 :
601 56 : if (hemisphere == HemiphereType::BOTH)
602 : {
603 55 : if (sExtentWGS84.MinY < 0 && sExtentWGS84.MaxY > 0)
604 : {
605 2 : std::vector<RPFFrameDef> res1, res2;
606 1 : if (nZoneIn == 0 || (nZoneIn >= MIN_ZONE_SOUTHERN_HEMISPHERE))
607 : {
608 2 : res1 = RPFGetCADRGFramesForEnvelope(
609 1 : nZoneIn, nReciprocalScale, poSrcDS, HemiphereType::SOUTH);
610 : }
611 1 : if (nZoneIn == 0 || (nZoneIn >= MIN_ZONE &&
612 : nZoneIn <= MAX_ZONE_NORTHERN_HEMISPHERE))
613 : {
614 0 : res2 = RPFGetCADRGFramesForEnvelope(
615 0 : nZoneIn, nReciprocalScale, poSrcDS, HemiphereType::NORTH);
616 : }
617 1 : res1.insert(res1.end(), res2.begin(), res2.end());
618 1 : return res1;
619 : }
620 54 : else if (sExtentWGS84.MinY >= 0)
621 50 : hemisphere = HemiphereType::NORTH;
622 : else
623 4 : hemisphere = HemiphereType::SOUTH;
624 : }
625 :
626 55 : CPLAssert(hemisphere == HemiphereType::NORTH ||
627 : hemisphere == HemiphereType::SOUTH);
628 :
629 55 : if (!(sExtentWGS84.MinX >= dfMinLonZone &&
630 53 : sExtentWGS84.MinX < dfMaxLonZone))
631 : {
632 2 : CPLError(CE_Failure, CPLE_AppDefined,
633 : "Minimum longitude of extent not in [%f,%f[ range",
634 : dfMinLonZone, dfMaxLonZone);
635 2 : return {};
636 : }
637 :
638 53 : double dfLastMaxLatZone = 0;
639 106 : std::vector<RPFFrameDef> res;
640 : const double dfYMinNorth = (hemisphere == HemiphereType::NORTH)
641 58 : ? std::max(0.0, sExtentWGS84.MinY)
642 5 : : std::max(0.0, -sExtentWGS84.MaxY);
643 53 : const double dfYMaxNorth = (hemisphere == HemiphereType::NORTH)
644 53 : ? sExtentWGS84.MaxY
645 5 : : -sExtentWGS84.MinY;
646 53 : const int nZoneOffset =
647 53 : (hemisphere == HemiphereType::NORTH) ? 0 : MAX_ZONE_NORTHERN_HEMISPHERE;
648 53 : CPLDebugOnly("CADRG", "Source minLat=%f, maxLat=%f", sExtentWGS84.MinY,
649 : sExtentWGS84.MaxY);
650 :
651 205 : for (int nZoneNorth = MIN_ZONE;
652 186 : nZoneNorth < MAX_ZONE_NORTHERN_HEMISPHERE &&
653 388 : nZoneIn != MAX_ZONE_NORTHERN_HEMISPHERE && nZoneIn != MAX_ZONE;
654 : ++nZoneNorth)
655 : {
656 181 : const int nZone = nZoneIn ? nZoneIn : nZoneNorth + nZoneOffset;
657 181 : const int nZoneIdx = (nZone - MIN_ZONE) % MAX_ZONE_NORTHERN_HEMISPHERE;
658 181 : const auto &sZoneDef = asARCZoneDefinitions[nZoneIdx];
659 181 : if (dfYMinNorth < sZoneDef.maxLat && dfYMaxNorth > sZoneDef.minLat)
660 : {
661 56 : const auto [dfMinLatZone, dfMaxLatZone] =
662 56 : GetMinMaxLatWithOverlap(nZone, nReciprocalScale);
663 56 : CPLDebugOnly("CADRG",
664 : "Zone %d: minLat_nominal=%d, maxLat_nominal=%d, "
665 : "minLat_overlap=%f, maxLat_overlap=%f",
666 : nZone,
667 : (hemisphere == HemiphereType::NORTH) ? sZoneDef.minLat
668 : : sZoneDef.maxLat,
669 : (hemisphere == HemiphereType::NORTH)
670 : ? sZoneDef.maxLat
671 : : -sZoneDef.minLat,
672 : dfMinLatZone, dfMaxLatZone);
673 : const double dfMaxLatZoneNorth =
674 56 : std::max(std::fabs(dfMinLatZone), std::fabs(dfMaxLatZone));
675 56 : if (dfMaxLatZoneNorth == dfLastMaxLatZone)
676 : {
677 : // Skip zone if fully within a previous one
678 : // This is the case of zones 5, 6 and 8 at scale 5 M.
679 : // See MIL-C-89038 page 70
680 6 : CPLDebugOnly(
681 : "CADRG",
682 : "Skipping zone %d as fully contained in previous one",
683 : nZone);
684 6 : continue;
685 : }
686 50 : dfLastMaxLatZone = dfMaxLatZoneNorth;
687 50 : const double latInterval = GetYPixelSize(nZone, nReciprocalScale);
688 50 : const double deltaLatFrame = latInterval * CADRG_FRAME_PIXEL_COUNT;
689 : const double dfFrameMinY =
690 50 : (std::max(sExtentWGS84.MinY, dfMinLatZone) - dfMinLatZone) /
691 50 : deltaLatFrame;
692 : const double dfFrameMaxY =
693 50 : (std::min(sExtentWGS84.MaxY, dfMaxLatZone) - dfMinLatZone) /
694 50 : deltaLatFrame;
695 50 : CPLDebugOnly("CADRG", "dfFrameMinY = %f, dfFrameMaxY=%f",
696 : dfFrameMinY, dfFrameMaxY);
697 50 : const int nFrameMinY = static_cast<int>(dfFrameMinY + EPSILON_1Em3);
698 50 : const int nFrameMaxY = static_cast<int>(dfFrameMaxY - EPSILON_1Em3);
699 :
700 50 : const double lonInterval = GetXPixelSize(nZone, nReciprocalScale);
701 50 : const double deltaLonFrame = lonInterval * CADRG_FRAME_PIXEL_COUNT;
702 : const double dfFrameMinX =
703 50 : (std::max(sExtentWGS84.MinX, dfMinLonZone) - dfMinLonZone) /
704 50 : deltaLonFrame;
705 : const double dfFrameMaxX =
706 50 : (std::min(sExtentWGS84.MaxX, dfMaxLonZone) - dfMinLonZone) /
707 50 : deltaLonFrame;
708 50 : CPLDebugOnly("CADRG", "dfFrameMinX = %f, dfFrameMaxX=%f",
709 : dfFrameMinX, dfFrameMaxX);
710 50 : const int nFrameMinX = static_cast<int>(dfFrameMinX + EPSILON_1Em3);
711 50 : const int nFrameMaxX = static_cast<int>(dfFrameMaxX - EPSILON_1Em3);
712 :
713 : RPFFrameDef sDef;
714 50 : sDef.nZone = nZone;
715 50 : sDef.nReciprocalScale = nReciprocalScale;
716 50 : sDef.nFrameMinX = nFrameMinX;
717 50 : sDef.nFrameMinY = nFrameMinY;
718 50 : sDef.nFrameMaxX = nFrameMaxX;
719 50 : sDef.nFrameMaxY = nFrameMaxY;
720 50 : sDef.dfResX = lonInterval;
721 50 : sDef.dfResY = latInterval;
722 50 : CPLDebug("CADRG", "Zone %d: frame (x,y)=(%d,%d) to (%d,%d)", nZone,
723 : nFrameMinX, nFrameMinY, nFrameMaxX, nFrameMaxY);
724 50 : res.push_back(sDef);
725 : }
726 175 : if (nZoneIn)
727 29 : break;
728 : }
729 :
730 : // Polar zone
731 53 : constexpr double MIN_POLAR_LAT_MANUAL = 70;
732 53 : constexpr double MIN_POLAR_LAT_AUTO = 80;
733 53 : if ((nZoneIn == 0 && dfYMaxNorth > MIN_POLAR_LAT_AUTO) ||
734 52 : (nZoneIn == MAX_ZONE_NORTHERN_HEMISPHERE + nZoneOffset &&
735 : dfYMaxNorth > MIN_POLAR_LAT_MANUAL))
736 : {
737 6 : const int nZone = MAX_ZONE_NORTHERN_HEMISPHERE + nZoneOffset;
738 12 : OGRSpatialReference oPolarSRS;
739 6 : oPolarSRS.importFromWkt(hemisphere == HemiphereType::NORTH
740 : ? pszNorthPolarProjection
741 : : pszSouthPolarProjection);
742 6 : const auto poSrcSRS = poSrcDS->GetSpatialRef();
743 : auto poCT = std::unique_ptr<OGRCoordinateTransformation>(
744 12 : OGRCreateCoordinateTransformation(poSrcSRS, &oPolarSRS));
745 6 : double dfXMin = 0;
746 6 : double dfYMin = 0;
747 6 : double dfXMax = 0;
748 6 : double dfYMax = 0;
749 :
750 6 : if (poSrcSRS->IsGeographic())
751 : {
752 1 : if (hemisphere == HemiphereType::NORTH)
753 0 : sExtentNativeCRS.MinY =
754 0 : std::max(MIN_POLAR_LAT_MANUAL, sExtentNativeCRS.MinY);
755 : else
756 1 : sExtentNativeCRS.MaxY =
757 1 : std::min(-MIN_POLAR_LAT_MANUAL, sExtentNativeCRS.MaxY);
758 : }
759 :
760 6 : if (poCT->TransformBounds(sExtentNativeCRS.MinX, sExtentNativeCRS.MinY,
761 : sExtentNativeCRS.MaxX, sExtentNativeCRS.MaxY,
762 : &dfXMin, &dfYMin, &dfXMax, &dfYMax,
763 6 : DEFAULT_DENSIFY_PTS))
764 : {
765 6 : const int numberFrames = GetPolarFrameCount(nReciprocalScale);
766 :
767 : // Cf MIL-C-89038 (CADRG specification), para 30.5
768 6 : const double R = numberFrames / 2.0 * CADRG_FRAME_PIXEL_COUNT;
769 :
770 : RPFFrameDef sDef;
771 6 : sDef.nZone = nZone;
772 6 : sDef.nReciprocalScale = nReciprocalScale;
773 6 : sDef.dfResX = GetXPixelSize(nZone, nReciprocalScale);
774 : // will lead to same value as dfResX
775 6 : sDef.dfResY = GetYPixelSize(nZone, nReciprocalScale);
776 6 : sDef.nFrameMinX =
777 12 : std::clamp(static_cast<int>((dfXMin / sDef.dfResX + R) /
778 6 : CADRG_FRAME_PIXEL_COUNT +
779 : EPSILON_1Em3),
780 6 : 0, numberFrames - 1);
781 6 : sDef.nFrameMinY =
782 12 : std::clamp(static_cast<int>((dfYMin / sDef.dfResY + R) /
783 6 : CADRG_FRAME_PIXEL_COUNT +
784 : EPSILON_1Em3),
785 6 : 0, numberFrames - 1);
786 6 : sDef.nFrameMaxX =
787 12 : std::clamp(static_cast<int>((dfXMax / sDef.dfResX + R) /
788 6 : CADRG_FRAME_PIXEL_COUNT -
789 : EPSILON_1Em3),
790 6 : 0, numberFrames - 1);
791 6 : sDef.nFrameMaxY =
792 12 : std::clamp(static_cast<int>((dfYMax / sDef.dfResY + R) /
793 6 : CADRG_FRAME_PIXEL_COUNT -
794 : EPSILON_1Em3),
795 6 : 0, numberFrames - 1);
796 6 : CPLDebug("CADRG", "Zone %d: frame (x,y)=(%d,%d) to (%d,%d)",
797 : sDef.nZone, sDef.nFrameMinX, sDef.nFrameMinY,
798 : sDef.nFrameMaxX, sDef.nFrameMaxY);
799 6 : res.push_back(sDef);
800 : }
801 0 : else if (nZoneIn == 0)
802 : {
803 0 : CPLError(CE_Failure, CPLE_AppDefined,
804 : "Cannot reproject source dataset extent to polar "
805 : "Azimuthal Equidistant");
806 : }
807 : }
808 :
809 53 : return res;
810 : }
811 :
812 : /** Given a dataset and a reciprocal scale (e.g. 1,000,000), returns the min/max
813 : * frame coordinate indices in all zones that intersect that area of interest
814 : * (when nZoneIn is 0), or in the specified zone (when nZoneIn is not 0)
815 : */
816 55 : std::vector<RPFFrameDef> RPFGetCADRGFramesForEnvelope(int nZoneIn,
817 : int nReciprocalScale,
818 : GDALDataset *poSrcDS)
819 : {
820 : return RPFGetCADRGFramesForEnvelope(nZoneIn, nReciprocalScale, poSrcDS,
821 55 : HemiphereType::BOTH);
822 : }
823 :
824 : /************************************************************************/
825 : /* RPFGetCADRGFrameNumberAsString() */
826 : /************************************************************************/
827 :
828 : /** Returns the 5 first character of the filename corresponding to the
829 : * frame specified by the provided parameters.
830 : */
831 96 : std::string RPFGetCADRGFrameNumberAsString(int nZone, int nReciprocalScale,
832 : int nFrameX, int nFrameY)
833 : {
834 : // Cf MIL-C-89038, page 60, 30.6 "Frame naming convention"
835 :
836 96 : const int nFrameCountAlongX = GetFrameCountAlongX(nZone, nReciprocalScale);
837 96 : const int nFrameCountAlongY = GetFrameCountAlongY(nZone, nReciprocalScale);
838 96 : CPLDebugOnly("CADRG",
839 : "Zone %d -> nFrameCountAlongX = %d, nFrameCountAlongY = %d",
840 : nZone, nFrameCountAlongX, nFrameCountAlongY);
841 96 : CPL_IGNORE_RET_VAL(nFrameCountAlongY);
842 96 : CPLAssert(nFrameX >= 0 && nFrameX < nFrameCountAlongX);
843 96 : CPLAssert(nFrameY >= 0 && nFrameY < nFrameCountAlongY);
844 96 : const int nFrameIdx = nFrameX + nFrameY * nFrameCountAlongX;
845 96 : CPLDebugOnly("CADRG", "Frame number (%d, %d) -> %d", nFrameX, nFrameY,
846 : nFrameIdx);
847 :
848 96 : std::string osRes;
849 :
850 96 : constexpr int BASE_34 = 34;
851 : // clang-format off
852 : // letters 'I' and 'O' are omitted to avoid confusiong with one and zero
853 96 : constexpr char ALPHABET_BASE_34[] = {
854 : '0', '1', '2', '3', '4', '5', '6', '7', '8', '9',
855 : 'A', 'B', 'C', 'D', 'E', 'F', 'G', 'H', 'J',
856 : 'K', 'L', 'M', 'N', 'P', 'Q', 'R', 'S', 'T',
857 : 'U', 'V', 'W', 'X', 'Y', 'Z'
858 : };
859 : // clang-format on
860 : static_assert(sizeof(ALPHABET_BASE_34) == BASE_34);
861 :
862 96 : int nCur = nFrameIdx;
863 36 : do
864 : {
865 132 : osRes += ALPHABET_BASE_34[nCur % BASE_34];
866 132 : nCur /= BASE_34;
867 132 : } while (nCur > 0);
868 :
869 96 : std::reverse(osRes.begin(), osRes.end());
870 : // Pad to 5 characters with leading zeroes.
871 444 : while (osRes.size() < 5)
872 348 : osRes = '0' + osRes;
873 192 : return osRes;
874 : }
875 :
876 : /************************************************************************/
877 : /* RPFGetCADRGFrameExtent() */
878 : /************************************************************************/
879 :
880 155 : void RPFGetCADRGFrameExtent(int nZone, int nReciprocalScale, int nFrameX,
881 : int nFrameY, double &dfXMin, double &dfYMin,
882 : double &dfXMax, double &dfYMax)
883 : {
884 155 : CPLAssert(RPFCADRGIsValidZone(nZone));
885 155 : const double dfXRes = GetXPixelSize(nZone, nReciprocalScale);
886 155 : const double dfYRes = GetYPixelSize(nZone, nReciprocalScale);
887 :
888 : double dfXMinLocal, dfYMinLocal;
889 155 : if (nZone == MAX_ZONE_NORTHERN_HEMISPHERE || nZone == MAX_ZONE)
890 : {
891 12 : const int numberFrames = GetPolarFrameCount(nReciprocalScale);
892 12 : const double dfXYOrigin =
893 12 : -numberFrames / 2.0 * dfXRes * CADRG_FRAME_PIXEL_COUNT;
894 12 : dfXMinLocal = dfXYOrigin + nFrameX * dfXRes * CADRG_FRAME_PIXEL_COUNT;
895 12 : dfYMinLocal = dfXYOrigin + nFrameY * dfYRes * CADRG_FRAME_PIXEL_COUNT;
896 : }
897 : else
898 : {
899 143 : const auto [dfMinLatZone, ignored] =
900 143 : GetMinMaxLatWithOverlap(nZone, nReciprocalScale);
901 143 : dfXMinLocal = dfMinLonZone + nFrameX * dfXRes * CADRG_FRAME_PIXEL_COUNT;
902 143 : dfYMinLocal = dfMinLatZone + nFrameY * dfYRes * CADRG_FRAME_PIXEL_COUNT;
903 : }
904 :
905 155 : dfXMin = dfXMinLocal;
906 155 : dfYMin = dfYMinLocal;
907 155 : dfXMax = dfXMinLocal + dfXRes * CADRG_FRAME_PIXEL_COUNT;
908 155 : dfYMax = dfYMinLocal + dfYRes * CADRG_FRAME_PIXEL_COUNT;
909 155 : }
910 :
911 : /************************************************************************/
912 : /* RPFGetCADRGClosestReciprocalScale() */
913 : /************************************************************************/
914 :
915 0 : int RPFGetCADRGClosestReciprocalScale(GDALDataset *poSrcDS,
916 : double dfDPIOverride, bool &bGotDPI)
917 : {
918 0 : bGotDPI = false;
919 :
920 0 : constexpr double INCH_TO_CM = 2.54;
921 0 : double dfDPI = dfDPIOverride;
922 0 : if (dfDPI <= 0)
923 : {
924 : const char *pszUnit =
925 0 : poSrcDS->GetMetadataItem("TIFFTAG_RESOLUTIONUNIT");
926 0 : const char *pszYRes = poSrcDS->GetMetadataItem("TIFFTAG_YRESOLUTION");
927 0 : if (pszUnit && pszYRes)
928 : {
929 0 : if (EQUAL(pszUnit, "2")) // Inch
930 0 : dfDPI = CPLAtof(pszYRes);
931 0 : else if (EQUAL(pszUnit, "3")) // Centimeter
932 : {
933 0 : dfDPI = CPLAtof(pszYRes) * INCH_TO_CM;
934 : }
935 : }
936 : }
937 0 : if (dfDPI <= 0)
938 : {
939 0 : CPLError(CE_Failure, CPLE_AppDefined,
940 : "Failed to estimate a reciprocal scale due to lack of DPI "
941 : "information. Please specify the DPI creation option.\n"
942 : "For reference, ADRG DPI is 254 and CADRG DPI is 169.333");
943 0 : return 0;
944 : }
945 :
946 0 : bGotDPI = true;
947 :
948 0 : constexpr double CADRG_DPI = INCH_TO_CM / CADRG_PITCH_IN_CM;
949 :
950 0 : GDALGeoTransform gt;
951 0 : if (poSrcDS->GetGeoTransform(gt) != CE_None)
952 : {
953 :
954 0 : CPLError(CE_Failure, CPLE_AppDefined, "Cannot get geotransform");
955 0 : return 0;
956 : }
957 :
958 0 : const double dfYRes = std::fabs(gt[5]);
959 0 : const double dfYResAtNominalCADRGDPI = dfYRes * dfDPI / CADRG_DPI;
960 :
961 0 : std::set<int> anSetReciprocalScales;
962 0 : for (int nIdx = 0;; ++nIdx)
963 : {
964 0 : const auto *psEntry = NITFGetRPFSeriesInfoFromIndex(nIdx);
965 0 : if (!psEntry)
966 0 : break;
967 0 : if (EQUAL(psEntry->rpfDataType, "CADRG"))
968 : {
969 : const int nVal =
970 0 : NITFGetScaleFromScaleResolution(psEntry->scaleResolution);
971 0 : if (nVal)
972 0 : anSetReciprocalScales.insert(nVal);
973 : }
974 0 : }
975 :
976 0 : int nCandidateReciprocalScale = 0;
977 0 : double dfBestProximityRatio = 0;
978 :
979 : // We tolerate up to a factor of 2 between the CADRG resolution at a
980 : //given scale and the dataset resolution
981 0 : constexpr double MAX_PROXIMITY_RATIO = 2;
982 :
983 0 : for (int nReciprocalScale : anSetReciprocalScales)
984 : {
985 : // This is actually zone independent
986 : const double dfLatInterval =
987 0 : GetYPixelSize(/* nZone = */ 1, nReciprocalScale);
988 0 : if (nCandidateReciprocalScale == 0 &&
989 0 : dfLatInterval / dfYResAtNominalCADRGDPI > MAX_PROXIMITY_RATIO)
990 : {
991 0 : break;
992 : }
993 0 : if (dfYResAtNominalCADRGDPI <= dfLatInterval)
994 : {
995 0 : const double dfThisProximityRatio =
996 : dfLatInterval / dfYResAtNominalCADRGDPI;
997 0 : if (nCandidateReciprocalScale == 0 ||
998 : dfThisProximityRatio < dfBestProximityRatio)
999 : {
1000 0 : nCandidateReciprocalScale = nReciprocalScale;
1001 : }
1002 0 : break;
1003 : }
1004 : else
1005 : {
1006 0 : const double dfThisProximityRatio =
1007 : dfYResAtNominalCADRGDPI / dfLatInterval;
1008 0 : if (dfThisProximityRatio < MAX_PROXIMITY_RATIO &&
1009 0 : (nCandidateReciprocalScale == 0 ||
1010 : dfThisProximityRatio < dfBestProximityRatio))
1011 : {
1012 0 : nCandidateReciprocalScale = nReciprocalScale;
1013 0 : dfBestProximityRatio = dfThisProximityRatio;
1014 : }
1015 : }
1016 : }
1017 :
1018 0 : if (nCandidateReciprocalScale == 0)
1019 : {
1020 0 : CPLError(CE_Failure, CPLE_AppDefined,
1021 : "Cannot find a pre-established scale matching source dataset "
1022 : "pixel size and scan resolution");
1023 : }
1024 0 : return nCandidateReciprocalScale;
1025 : }
1026 :
1027 : /************************************************************************/
1028 : /* Create_CADRG_CoverageSection() */
1029 : /************************************************************************/
1030 :
1031 33 : static bool Create_CADRG_CoverageSection(
1032 : GDALOffsetPatcher::OffsetPatcher *offsetPatcher, GDALDataset *poSrcDS,
1033 : const std::string &osFilename, const CPLStringList &aosOptions,
1034 : int nReciprocalScale)
1035 : {
1036 33 : auto poBuffer = offsetPatcher->CreateBuffer(
1037 : "CoverageSectionSubheader", /* bEndiannessIsLittle = */ false);
1038 33 : CPLAssert(poBuffer);
1039 33 : poBuffer->DeclareOffsetAtCurrentPosition("COVERAGE_SECTION_LOCATION");
1040 :
1041 33 : OGREnvelope sExtent;
1042 33 : const auto poSrcSRS = poSrcDS->GetSpatialRef();
1043 33 : if (!poSrcSRS || poSrcDS->GetExtent(&sExtent) != CE_None)
1044 : {
1045 0 : CPLAssert(false); // already checked in calling sites
1046 : return false;
1047 : }
1048 :
1049 33 : double dfULX = sExtent.MinX;
1050 33 : double dfULY = sExtent.MaxY;
1051 33 : double dfLLX = sExtent.MinX;
1052 33 : double dfLLY = sExtent.MinY;
1053 33 : double dfURX = sExtent.MaxX;
1054 33 : double dfURY = sExtent.MaxY;
1055 33 : double dfLRX = sExtent.MaxX;
1056 33 : double dfLRY = sExtent.MinY;
1057 :
1058 66 : OGRSpatialReference oSRS_WGS84;
1059 33 : oSRS_WGS84.SetWellKnownGeogCS("WGS84");
1060 33 : oSRS_WGS84.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
1061 33 : if (!poSrcSRS->IsSame(&oSRS_WGS84))
1062 : {
1063 : auto poCT = std::unique_ptr<OGRCoordinateTransformation>(
1064 4 : OGRCreateCoordinateTransformation(poSrcSRS, &oSRS_WGS84));
1065 8 : if (!(poCT && poCT->Transform(1, &dfULX, &dfULY) &&
1066 4 : poCT->Transform(1, &dfLLX, &dfLLY) &&
1067 4 : poCT->Transform(1, &dfURX, &dfURY) &&
1068 4 : poCT->Transform(1, &dfLRX, &dfLRY)))
1069 : {
1070 0 : CPLError(CE_Failure, CPLE_AppDefined,
1071 : "Create_CADRG_CoverageSection(): cannot compute corner "
1072 : "lat/lon");
1073 0 : return false;
1074 : }
1075 : }
1076 :
1077 264 : const auto RoundIfCloseToInt = [](double dfX)
1078 : {
1079 264 : double dfRounded = std::round(dfX);
1080 264 : if (std::abs(dfX - dfRounded) < 1e-12)
1081 108 : return dfRounded;
1082 156 : return dfX;
1083 : };
1084 :
1085 132 : const auto NormalizeToMinusPlus180 = [](double dfLon)
1086 : {
1087 132 : constexpr double EPSILON_SMALL = 1e-9;
1088 132 : if (dfLon < -180 - EPSILON_SMALL)
1089 0 : dfLon += 360;
1090 132 : else if (dfLon < -180)
1091 0 : dfLon = -180;
1092 132 : else if (dfLon > 180 + EPSILON_SMALL)
1093 0 : dfLon -= 360;
1094 132 : else if (dfLon > 180)
1095 4 : dfLon = 180;
1096 132 : return dfLon;
1097 : };
1098 :
1099 33 : int nZone = atoi(aosOptions.FetchNameValueDef("ZONE", "0"));
1100 33 : if (nZone < MIN_ZONE || nZone > MAX_ZONE)
1101 : {
1102 3 : const std::string osExt = CPLGetExtensionSafe(osFilename.c_str());
1103 3 : if (osExt.size() == 3)
1104 3 : nZone = RPFCADRGZoneCharToNum(osExt.back());
1105 3 : if (nZone < MIN_ZONE || nZone > MAX_ZONE)
1106 : {
1107 0 : const double dfMeanLat = (dfULY + dfLLY) / 2;
1108 0 : nZone = GetARCZoneFromLat(dfMeanLat);
1109 0 : if (nZone == 0)
1110 0 : return false;
1111 : }
1112 : }
1113 :
1114 : // Upper left corner lat, lon
1115 33 : poBuffer->AppendFloat64(RoundIfCloseToInt(dfULY));
1116 33 : poBuffer->AppendFloat64(RoundIfCloseToInt(NormalizeToMinusPlus180(dfULX)));
1117 : // Lower left corner lat, lon
1118 33 : poBuffer->AppendFloat64(RoundIfCloseToInt(dfLLY));
1119 33 : poBuffer->AppendFloat64(RoundIfCloseToInt(NormalizeToMinusPlus180(dfLLX)));
1120 : // Upper right corner lat, lon
1121 33 : poBuffer->AppendFloat64(RoundIfCloseToInt(dfURY));
1122 33 : poBuffer->AppendFloat64(RoundIfCloseToInt(NormalizeToMinusPlus180(dfURX)));
1123 : // Lower right corner lat, lon
1124 33 : poBuffer->AppendFloat64(RoundIfCloseToInt(dfLRY));
1125 33 : poBuffer->AppendFloat64(RoundIfCloseToInt(NormalizeToMinusPlus180(dfLRX)));
1126 :
1127 33 : double latResolution = 0;
1128 33 : double lonResolution = 0;
1129 33 : double latInterval = 0;
1130 33 : double lonInterval = 0;
1131 33 : RPFGetCADRGResolutionAndInterval(nZone, nReciprocalScale, latResolution,
1132 : lonResolution, latInterval, lonInterval);
1133 :
1134 33 : poBuffer->AppendFloat64(latResolution);
1135 33 : poBuffer->AppendFloat64(lonResolution);
1136 33 : if (!poSrcSRS->IsSame(&oSRS_WGS84))
1137 : {
1138 4 : poBuffer->AppendFloat64(latInterval);
1139 4 : poBuffer->AppendFloat64(lonInterval);
1140 : }
1141 : else
1142 : {
1143 29 : GDALGeoTransform gt;
1144 29 : CPL_IGNORE_RET_VAL(poSrcDS->GetGeoTransform(gt));
1145 :
1146 : // Theoretical value: latInterval = 90.0 / latCst_CADRG;
1147 29 : poBuffer->AppendFloat64(std::fabs(gt[5]));
1148 : // Theoretical value: lonInterval = 360.0 / lonCst_CADRG;
1149 29 : poBuffer->AppendFloat64(gt[1]);
1150 : }
1151 33 : return true;
1152 : }
1153 :
1154 : /************************************************************************/
1155 : /* Create_CADRG_ColorGrayscaleSection() */
1156 : /************************************************************************/
1157 :
1158 : constexpr GByte NUM_COLOR_TABLES = 3;
1159 : constexpr GByte NUM_COLOR_CONVERTERS = 2;
1160 :
1161 33 : static void Create_CADRG_ColorGrayscaleSection(
1162 : GDALOffsetPatcher::OffsetPatcher *offsetPatcher)
1163 : {
1164 33 : auto poBuffer = offsetPatcher->CreateBuffer(
1165 : "ColorGrayscaleSectionSubheader", /* bEndiannessIsLittle = */ false);
1166 33 : CPLAssert(poBuffer);
1167 33 : poBuffer->DeclareOffsetAtCurrentPosition("COLOR_GRAYSCALE_LOCATION");
1168 33 : poBuffer->AppendByte(NUM_COLOR_TABLES);
1169 33 : poBuffer->AppendByte(NUM_COLOR_CONVERTERS);
1170 : // EXTERNAL_COLOR_GRAYSCALE_FILENAME
1171 33 : poBuffer->AppendString(std::string(12, ' '));
1172 33 : }
1173 :
1174 : /************************************************************************/
1175 : /* Create_CADRG_ImageDescriptionSection() */
1176 : /************************************************************************/
1177 :
1178 33 : static void Create_CADRG_ImageDescriptionSection(
1179 : GDALOffsetPatcher::OffsetPatcher *offsetPatcher, GDALDataset *poSrcDS,
1180 : bool bHasTransparentPixels)
1181 : {
1182 33 : CPLAssert((poSrcDS->GetRasterXSize() % BLOCK_SIZE) == 0);
1183 33 : CPLAssert((poSrcDS->GetRasterYSize() % BLOCK_SIZE) == 0);
1184 33 : CPLAssert(poSrcDS->GetRasterXSize() <= UINT16_MAX * BLOCK_SIZE);
1185 33 : CPLAssert(poSrcDS->GetRasterYSize() <= UINT16_MAX * BLOCK_SIZE);
1186 : const uint16_t nSubFramesPerRow =
1187 33 : static_cast<uint16_t>(poSrcDS->GetRasterXSize() / BLOCK_SIZE);
1188 : const uint16_t nSubFramesPerCol =
1189 33 : static_cast<uint16_t>(poSrcDS->GetRasterYSize() / BLOCK_SIZE);
1190 33 : CPLAssert(nSubFramesPerRow * nSubFramesPerCol < UINT16_MAX);
1191 :
1192 33 : auto poBuffer = offsetPatcher->CreateBuffer(
1193 : "ImageDescriptionSubheader", /* bEndiannessIsLittle = */ false);
1194 33 : CPLAssert(poBuffer);
1195 33 : poBuffer->DeclareOffsetAtCurrentPosition(
1196 : "IMAGE_DESCRIPTION_SECTION_LOCATION");
1197 33 : poBuffer->AppendUInt16(1); // NUMBER_OF_SPECTRAL_GROUPS
1198 33 : poBuffer->AppendUInt16(static_cast<uint16_t>(
1199 33 : nSubFramesPerRow * nSubFramesPerCol)); // NUMBER_OF_SUBFRAME_TABLES
1200 33 : poBuffer->AppendUInt16(1); // NUMBER_OF_SPECTRAL_BAND_TABLES
1201 33 : poBuffer->AppendUInt16(1); // NUMBER_OF_SPECTRAL_BAND_LINES_PER_IMAGE_ROW
1202 33 : poBuffer->AppendUInt16(
1203 : nSubFramesPerRow); // NUMBER_OF_SUBFRAME_IN_EAST_WEST_DIRECTION
1204 33 : poBuffer->AppendUInt16(
1205 : nSubFramesPerCol); // NUMBER_OF_SUBFRAME_IN_NORTH_SOUTH_DIRECTION
1206 33 : poBuffer->AppendUInt32(
1207 : BLOCK_SIZE); // NUMBER_OF_OUTPUT_COLUMNS_PER_SUBFRAME
1208 33 : poBuffer->AppendUInt32(BLOCK_SIZE); // NUMBER_OF_OUTPUT_ROWS_PER_SUBFRAME
1209 33 : poBuffer->AppendUInt32(UINT32_MAX); // SUBFRAME_MASK_TABLE_OFFSET
1210 33 : if (!bHasTransparentPixels)
1211 : {
1212 4 : poBuffer->AppendUInt32(UINT32_MAX); // TRANSPARENCY_MASK_TABLE_OFFSET
1213 : }
1214 : else
1215 : {
1216 : // Offset in bytes from the beginning of the mask subsection to the
1217 : // beginning of the transparency mask table.
1218 29 : poBuffer->AppendUInt32(7);
1219 : }
1220 33 : }
1221 :
1222 : /************************************************************************/
1223 : /* Create_CADRG_ColormapSection() */
1224 : /************************************************************************/
1225 :
1226 : constexpr uint16_t SZ_UINT16 = 2;
1227 : constexpr uint16_t SZ_UINT32 = 4;
1228 : constexpr uint32_t COLORMAP_OFFSET_TABLE_OFFSET = SZ_UINT32 + SZ_UINT16;
1229 : constexpr uint16_t COLOR_GRAYSCALE_OFFSET_RECORD_LENGTH = 17;
1230 :
1231 33 : static void Create_CADRG_ColormapSection(
1232 : GDALOffsetPatcher::OffsetPatcher *offsetPatcher, GDALDataset *poSrcDS,
1233 : bool bHasTransparentPixels,
1234 : const std::vector<BucketItem<ColorTableBased4x4Pixels>> &codebook,
1235 : const CADRGCreateCopyContext ©Context)
1236 : {
1237 33 : auto poBuffer = offsetPatcher->CreateBuffer(
1238 : "ColormapSubsection", /* bEndiannessIsLittle = */ false);
1239 33 : CPLAssert(poBuffer);
1240 33 : poBuffer->DeclareOffsetAtCurrentPosition("COLORMAP_LOCATION");
1241 33 : poBuffer->AppendUInt32(COLORMAP_OFFSET_TABLE_OFFSET);
1242 33 : poBuffer->AppendUInt16(COLOR_GRAYSCALE_OFFSET_RECORD_LENGTH);
1243 33 : CPLAssert(poBuffer->GetBuffer().size() == COLORMAP_OFFSET_TABLE_OFFSET);
1244 :
1245 33 : uint32_t nColorTableOffset =
1246 : COLORMAP_OFFSET_TABLE_OFFSET +
1247 : NUM_COLOR_TABLES * COLOR_GRAYSCALE_OFFSET_RECORD_LENGTH;
1248 : // 4=R,G,B,M
1249 33 : constexpr GByte COLOR_TABLE_ENTRY_SIZE = 4;
1250 33 : uint32_t nHistogramTableOffset =
1251 : COLORMAP_OFFSET_TABLE_OFFSET +
1252 : NUM_COLOR_TABLES * COLOR_GRAYSCALE_OFFSET_RECORD_LENGTH +
1253 : (CADRG_MAX_COLOR_ENTRY_COUNT + CADRG_SECOND_CT_COUNT +
1254 : CADRG_THIRD_CT_COUNT) *
1255 : COLOR_TABLE_ENTRY_SIZE;
1256 132 : for (int i = 0; i < NUM_COLOR_TABLES; ++i)
1257 : {
1258 99 : poBuffer->AppendUInt16(2); // color/grayscale table id
1259 : // number of colors
1260 99 : const uint32_t nColorCount = (i == 0) ? CADRG_MAX_COLOR_ENTRY_COUNT
1261 : : (i == 1) ? CADRG_SECOND_CT_COUNT
1262 : : CADRG_THIRD_CT_COUNT;
1263 99 : poBuffer->AppendUInt32(nColorCount);
1264 99 : poBuffer->AppendByte(
1265 : COLOR_TABLE_ENTRY_SIZE); // color/grayscale element length
1266 99 : poBuffer->AppendUInt16(static_cast<uint16_t>(
1267 : sizeof(uint32_t))); // histogram record length
1268 : // color/grayscale table offset
1269 99 : poBuffer->AppendUInt32(nColorTableOffset);
1270 99 : nColorTableOffset += nColorCount * COLOR_TABLE_ENTRY_SIZE;
1271 : // histogram table offset
1272 99 : poBuffer->AppendUInt32(nHistogramTableOffset);
1273 99 : nHistogramTableOffset +=
1274 99 : nColorCount * static_cast<uint32_t>(sizeof(uint32_t));
1275 : }
1276 33 : CPLAssert(poBuffer->GetBuffer().size() ==
1277 : COLORMAP_OFFSET_TABLE_OFFSET +
1278 : NUM_COLOR_TABLES * COLOR_GRAYSCALE_OFFSET_RECORD_LENGTH);
1279 :
1280 : // Write color tables
1281 132 : for (int iCT = 0; iCT < NUM_COLOR_TABLES; ++iCT)
1282 : {
1283 : const auto poCT = (iCT == 0)
1284 165 : ? poSrcDS->GetRasterBand(1)->GetColorTable()
1285 66 : : (iCT == 1) ? &(copyContext.oCT2)
1286 99 : : &(copyContext.oCT3);
1287 99 : const int nMaxCTEntries =
1288 : (iCT == 0)
1289 165 : ? CADRG_MAX_COLOR_ENTRY_COUNT + (bHasTransparentPixels ? 1 : 0)
1290 66 : : (iCT == 1) ? CADRG_SECOND_CT_COUNT
1291 : : CADRG_THIRD_CT_COUNT;
1292 8840 : for (int i = 0; i < nMaxCTEntries; ++i)
1293 : {
1294 8741 : if (i < poCT->GetColorEntryCount())
1295 : {
1296 8074 : const auto psEntry = poCT->GetColorEntry(i);
1297 8074 : poBuffer->AppendByte(static_cast<GByte>(psEntry->c1));
1298 8074 : poBuffer->AppendByte(static_cast<GByte>(psEntry->c2));
1299 8074 : poBuffer->AppendByte(static_cast<GByte>(psEntry->c3));
1300 : // Standard formula to convert R,G,B to gray scale level
1301 8074 : const int M = (psEntry->c1 * 299 + psEntry->c2 * 587 +
1302 8074 : psEntry->c3 * 114) /
1303 : 1000;
1304 8074 : poBuffer->AppendByte(static_cast<GByte>(M));
1305 : }
1306 : else
1307 : {
1308 667 : poBuffer->AppendUInt32(0);
1309 : }
1310 : }
1311 : }
1312 :
1313 : // Compute the number of pixels in the output image per colormap entry
1314 : // (exclude the entry for transparent pixels)
1315 66 : std::vector<uint32_t> anHistogram(CADRG_MAX_COLOR_ENTRY_COUNT);
1316 66 : std::vector<uint32_t> anHistogram2(CADRG_SECOND_CT_COUNT);
1317 66 : std::vector<uint32_t> anHistogram3(CADRG_THIRD_CT_COUNT);
1318 33 : size_t nTotalCount = 0;
1319 135201 : for (const auto &[i, item] : cpl::enumerate(codebook))
1320 : {
1321 135168 : if (bHasTransparentPixels && i == TRANSPARENT_CODEBOOK_CODE)
1322 : {
1323 29 : nTotalCount += item.m_count;
1324 : }
1325 : else
1326 : {
1327 2297360 : for (GByte byVal : item.m_vec.vals())
1328 : {
1329 2162220 : anHistogram[byVal] += item.m_count;
1330 2162220 : nTotalCount += item.m_count;
1331 2162220 : anHistogram2[copyContext.anMapCT1ToCT2[byVal]] += item.m_count;
1332 2162220 : anHistogram3[copyContext.anMapCT1ToCT3[byVal]] += item.m_count;
1333 : }
1334 : }
1335 : }
1336 33 : CPLAssert(nTotalCount == static_cast<size_t>(poSrcDS->GetRasterXSize()) *
1337 : poSrcDS->GetRasterYSize());
1338 33 : CPL_IGNORE_RET_VAL(nTotalCount);
1339 :
1340 : // Write histograms
1341 7161 : for (auto nCount : anHistogram)
1342 : {
1343 7128 : poBuffer->AppendUInt32(nCount);
1344 : }
1345 1089 : for (auto nCount : anHistogram2)
1346 : {
1347 1056 : poBuffer->AppendUInt32(nCount);
1348 : }
1349 561 : for (auto nCount : anHistogram3)
1350 : {
1351 528 : poBuffer->AppendUInt32(nCount);
1352 : }
1353 33 : }
1354 :
1355 : /************************************************************************/
1356 : /* Create_CADRG_ColorConverterSubsection() */
1357 : /************************************************************************/
1358 :
1359 33 : static void Create_CADRG_ColorConverterSubsection(
1360 : GDALOffsetPatcher::OffsetPatcher *offsetPatcher, bool bHasTransparentPixels,
1361 : const CADRGCreateCopyContext ©Context)
1362 : {
1363 33 : auto poBuffer = offsetPatcher->CreateBuffer(
1364 : "ColorConverterSubsection", /* bEndiannessIsLittle = */ false);
1365 33 : CPLAssert(poBuffer);
1366 33 : poBuffer->DeclareOffsetAtCurrentPosition("COLOR_CONVERTER_SUBSECTION");
1367 :
1368 33 : constexpr uint32_t COLOR_CONVERTER_OFFSET_TABLE_OFFSET =
1369 : SZ_UINT32 + SZ_UINT16 + SZ_UINT16;
1370 33 : poBuffer->AppendUInt32(COLOR_CONVERTER_OFFSET_TABLE_OFFSET);
1371 33 : constexpr uint16_t COLOR_CONVERTER_OFFSET_RECORD_LENGTH =
1372 : SZ_UINT16 + SZ_UINT32 + SZ_UINT32 + SZ_UINT32;
1373 33 : poBuffer->AppendUInt16(COLOR_CONVERTER_OFFSET_RECORD_LENGTH);
1374 33 : constexpr uint16_t COLOR_CONVERTER_RECORD_LENGTH = SZ_UINT32;
1375 33 : poBuffer->AppendUInt16(COLOR_CONVERTER_RECORD_LENGTH);
1376 33 : const int numberColorConverterRecords =
1377 33 : CADRG_MAX_COLOR_ENTRY_COUNT + (bHasTransparentPixels ? 1 : 0);
1378 :
1379 99 : for (int i = 0; i < NUM_COLOR_CONVERTERS; ++i)
1380 : {
1381 66 : constexpr uint16_t COLOR_CONVERTER_TABLE_ID = 5;
1382 66 : poBuffer->AppendUInt16(COLOR_CONVERTER_TABLE_ID);
1383 66 : poBuffer->AppendUInt32(numberColorConverterRecords);
1384 66 : uint32_t colorConverterTableOffset =
1385 : COLOR_CONVERTER_OFFSET_TABLE_OFFSET +
1386 : 2 * COLOR_CONVERTER_OFFSET_RECORD_LENGTH +
1387 66 : i * numberColorConverterRecords * COLOR_CONVERTER_RECORD_LENGTH;
1388 66 : poBuffer->AppendUInt32(colorConverterTableOffset);
1389 66 : constexpr uint32_t SOURCE_OFFSET_TABLE_OFFSET =
1390 : COLORMAP_OFFSET_TABLE_OFFSET;
1391 66 : poBuffer->AppendUInt32(SOURCE_OFFSET_TABLE_OFFSET);
1392 66 : const uint32_t targetOffsetTableOffset =
1393 : COLORMAP_OFFSET_TABLE_OFFSET +
1394 66 : i * COLOR_GRAYSCALE_OFFSET_RECORD_LENGTH;
1395 66 : poBuffer->AppendUInt32(targetOffsetTableOffset);
1396 : }
1397 :
1398 99 : for (int iCvt = 0; iCvt < NUM_COLOR_CONVERTERS; ++iCvt)
1399 : {
1400 14380 : for (int j = 0; j < numberColorConverterRecords; ++j)
1401 : {
1402 14314 : if (j == CADRG_MAX_COLOR_ENTRY_COUNT)
1403 : {
1404 : // It is not specified what we should do about the transparent
1405 : // entry...
1406 58 : poBuffer->AppendUInt32(0);
1407 : }
1408 : else
1409 : {
1410 28512 : poBuffer->AppendUInt32((iCvt == 0)
1411 14256 : ? copyContext.anMapCT1ToCT2[j]
1412 7128 : : copyContext.anMapCT1ToCT3[j]);
1413 : }
1414 : }
1415 : }
1416 33 : }
1417 :
1418 : /************************************************************************/
1419 : /* Perform_CADRG_VQ_Compression() */
1420 : /************************************************************************/
1421 :
1422 33 : static bool Perform_CADRG_VQ_Compression(
1423 : GDALDataset *poSrcDS,
1424 : std::vector<BucketItem<ColorTableBased4x4Pixels>> &codebook,
1425 : std::vector<short> &VQImage, bool &bHasTransparentPixels)
1426 : {
1427 33 : const int nY = poSrcDS->GetRasterYSize();
1428 33 : const int nX = poSrcDS->GetRasterXSize();
1429 33 : CPLAssert((nY % SUBSAMPLING) == 0);
1430 33 : CPLAssert((nX % SUBSAMPLING) == 0);
1431 33 : CPLAssert(nX < INT_MAX / nY);
1432 :
1433 33 : auto poBand = poSrcDS->GetRasterBand(1);
1434 :
1435 66 : std::vector<GByte> pixels;
1436 33 : if (poBand->ReadRaster(pixels) != CE_None)
1437 0 : return false;
1438 :
1439 33 : const auto poCT = poBand->GetColorTable();
1440 33 : CPLAssert(poCT);
1441 66 : std::vector<GByte> vR, vG, vB;
1442 : const int nColorCount =
1443 33 : std::min(CADRG_MAX_COLOR_ENTRY_COUNT, poCT->GetColorEntryCount());
1444 : const bool bHasTransparentEntry =
1445 33 : poCT->GetColorEntryCount() >= CADRG_MAX_COLOR_ENTRY_COUNT + 1 &&
1446 32 : poCT->GetColorEntry(TRANSPARENT_COLOR_TABLE_ENTRY)->c1 == 0 &&
1447 32 : poCT->GetColorEntry(TRANSPARENT_COLOR_TABLE_ENTRY)->c2 == 0 &&
1448 97 : poCT->GetColorEntry(TRANSPARENT_COLOR_TABLE_ENTRY)->c3 == 0 &&
1449 32 : poCT->GetColorEntry(TRANSPARENT_COLOR_TABLE_ENTRY)->c4 == 0;
1450 7161 : for (int i = 0; i < nColorCount; ++i)
1451 : {
1452 7128 : const auto entry = poCT->GetColorEntry(i);
1453 7128 : vR.push_back(static_cast<GByte>(entry->c1));
1454 7128 : vG.push_back(static_cast<GByte>(entry->c2));
1455 7128 : vB.push_back(static_cast<GByte>(entry->c3));
1456 : }
1457 33 : ColorTableBased4x4Pixels ctxt(vR, vG, vB);
1458 :
1459 : struct Occurrences
1460 : {
1461 : // number of 4x4 pixel blocks using those 4x4 pixel values
1462 : int nCount = 0;
1463 : // Point to indices in the output image that use that 4x4 pixel blocks
1464 : std::vector<int> anIndicesToOutputImage{};
1465 : };
1466 :
1467 33 : const int nVQImgHeight = nY / SUBSAMPLING;
1468 33 : const int nVQImgWidth = nX / SUBSAMPLING;
1469 33 : CPLAssert(nVQImgWidth > 0);
1470 33 : CPLAssert(nVQImgHeight > 0);
1471 33 : CPLAssert(nVQImgWidth < INT_MAX / nVQImgHeight);
1472 33 : VQImage.resize(nVQImgHeight * nVQImgWidth);
1473 :
1474 : // Collect all the occurrences of 4x4 pixel values into a map indexed by them
1475 66 : std::map<Vector<ColorTableBased4x4Pixels>, Occurrences> vectorMap;
1476 33 : bHasTransparentPixels = false;
1477 : std::array<GByte, SUBSAMPLING * SUBSAMPLING> vals;
1478 33 : std::fill(vals.begin(), vals.end(), static_cast<GByte>(0));
1479 33 : int nTransparentPixels = 0;
1480 12705 : for (int j = 0, nOutputIdx = 0; j < nVQImgHeight; ++j)
1481 : {
1482 4878720 : for (int i = 0; i < nVQImgWidth; ++i, ++nOutputIdx)
1483 : {
1484 24330200 : for (int y = 0; y < SUBSAMPLING; ++y)
1485 : {
1486 97321000 : for (int x = 0; x < SUBSAMPLING; ++x)
1487 : {
1488 77856800 : const GByte val = pixels[(j * SUBSAMPLING + y) * nX +
1489 77856800 : (i * SUBSAMPLING + x)];
1490 77856800 : if (bHasTransparentEntry &&
1491 : val == TRANSPARENT_COLOR_TABLE_ENTRY)
1492 : {
1493 : // As soon as one of the pixels in the 4x4 block is
1494 : // transparent, the whole block is flagged as transparent
1495 56641900 : bHasTransparentPixels = true;
1496 56641900 : VQImage[nOutputIdx] = TRANSPARENT_CODEBOOK_CODE;
1497 : }
1498 : else
1499 : {
1500 21214900 : if (val >= nColorCount)
1501 : {
1502 0 : CPLError(CE_Failure, CPLE_AppDefined,
1503 : "Out of range pixel value found: %d", val);
1504 0 : return false;
1505 : }
1506 21214900 : vals[SUBSAMPLING * y + x] = val;
1507 : }
1508 : }
1509 : }
1510 4866050 : if (VQImage[nOutputIdx] == TRANSPARENT_CODEBOOK_CODE)
1511 : {
1512 3544380 : nTransparentPixels += SUBSAMPLING * SUBSAMPLING;
1513 : }
1514 : else
1515 : {
1516 1321660 : auto &elt = vectorMap[Vector<ColorTableBased4x4Pixels>(vals)];
1517 1321660 : ++elt.nCount;
1518 1321660 : elt.anIndicesToOutputImage.push_back(nOutputIdx);
1519 : }
1520 : }
1521 : }
1522 :
1523 : // Convert that map into a std::vector
1524 66 : std::vector<BucketItem<ColorTableBased4x4Pixels>> vectors;
1525 33 : vectors.reserve(vectorMap.size());
1526 154153 : for (auto &[key, value] : vectorMap)
1527 : {
1528 154120 : vectors.emplace_back(key, value.nCount,
1529 154120 : std::move(value.anIndicesToOutputImage));
1530 : }
1531 33 : vectorMap.clear();
1532 :
1533 : // Create the KD-Tree
1534 66 : PNNKDTree<ColorTableBased4x4Pixels> kdtree;
1535 :
1536 : // Insert the initial items
1537 33 : const bool bEmptyImage = vectors.empty();
1538 33 : int nCodeCount = kdtree.insert(std::move(vectors), ctxt);
1539 33 : if (!bEmptyImage && nCodeCount == 0)
1540 0 : return false;
1541 :
1542 : // Reduce to the maximum target
1543 33 : const int nMaxCodes =
1544 33 : bHasTransparentPixels ? CODEBOOK_MAX_SIZE - 1 : CODEBOOK_MAX_SIZE;
1545 33 : if (nCodeCount > nMaxCodes)
1546 : {
1547 7 : const int nNewCodeCount = kdtree.cluster(nCodeCount, nMaxCodes, ctxt);
1548 7 : if (nNewCodeCount == 0)
1549 0 : return false;
1550 7 : CPLDebug("NITF", "VQ compression: reducing from %d codes to %d",
1551 : nCodeCount, nNewCodeCount);
1552 : }
1553 : else
1554 : {
1555 26 : CPLDebug("NITF",
1556 : "Already less than %d codes. VQ compression is lossless",
1557 : nMaxCodes);
1558 : }
1559 :
1560 : // Create the code book and the target VQ-compressed image.
1561 33 : codebook.reserve(CODEBOOK_MAX_SIZE);
1562 33 : kdtree.iterateOverLeaves(
1563 1426900 : [&codebook, &VQImage](PNNKDTree<ColorTableBased4x4Pixels> &node)
1564 : {
1565 57444 : for (auto &item : node.bucketItems())
1566 : {
1567 47791 : const int i = static_cast<int>(codebook.size());
1568 1369450 : for (const auto idx : item.m_origVectorIndices)
1569 : {
1570 1321660 : VQImage[idx] = static_cast<short>(i);
1571 : }
1572 47791 : codebook.push_back(std::move(item));
1573 : }
1574 9653 : });
1575 :
1576 : // Add dummy entries until CODEBOOK_MAX_SIZE is reached. In theory we
1577 : // could provide less code if we don't reach it, but for broader
1578 : // compatibility it seems best to go up to the typical max value.
1579 : // Furthermore when there is transparency, the CADRG spec mentions 4095
1580 : // to be reserved for a transparent 4x4 kernel pointing to color table entry
1581 : // 216.
1582 87410 : while (codebook.size() < CODEBOOK_MAX_SIZE)
1583 : {
1584 : codebook.emplace_back(
1585 87377 : Vector<ColorTableBased4x4Pixels>(
1586 : filled_array<GByte,
1587 0 : Vector<ColorTableBased4x4Pixels>::PIX_COUNT>(0)),
1588 174754 : 0, std::vector<int>());
1589 : }
1590 :
1591 33 : if (bHasTransparentPixels)
1592 : {
1593 29 : codebook[TRANSPARENT_CODEBOOK_CODE]
1594 29 : .m_vec = Vector<ColorTableBased4x4Pixels>(
1595 29 : filled_array<GByte, Vector<ColorTableBased4x4Pixels>::PIX_COUNT>(
1596 29 : static_cast<GByte>(TRANSPARENT_COLOR_TABLE_ENTRY)));
1597 29 : codebook[TRANSPARENT_CODEBOOK_CODE].m_count = nTransparentPixels;
1598 29 : codebook[TRANSPARENT_CODEBOOK_CODE].m_origVectorIndices.clear();
1599 : }
1600 :
1601 33 : return true;
1602 : }
1603 :
1604 : /************************************************************************/
1605 : /* RPFFrameCreateCADRG_TREs() */
1606 : /************************************************************************/
1607 :
1608 : std::unique_ptr<CADRGInformation>
1609 33 : RPFFrameCreateCADRG_TREs(GDALOffsetPatcher::OffsetPatcher *offsetPatcher,
1610 : const std::string &osFilename, GDALDataset *poSrcDS,
1611 : CPLStringList &aosOptions,
1612 : const CADRGCreateCopyContext ©Context)
1613 : {
1614 66 : auto priv = std::make_unique<CADRGInformation::Private>();
1615 33 : if (!Perform_CADRG_VQ_Compression(poSrcDS, priv->codebook, priv->VQImage,
1616 33 : priv->bHasTransparentPixels))
1617 : {
1618 0 : return nullptr;
1619 : }
1620 :
1621 33 : Create_CADRG_RPFHDR(offsetPatcher, osFilename, aosOptions);
1622 :
1623 : // Create buffers that will be written into file by RPFFrameWriteCADRG_RPFIMG()s
1624 33 : Create_CADRG_LocationComponent(offsetPatcher);
1625 33 : if (!Create_CADRG_CoverageSection(offsetPatcher, poSrcDS, osFilename,
1626 33 : aosOptions, copyContext.nReciprocalScale))
1627 0 : return nullptr;
1628 33 : Create_CADRG_ColorGrayscaleSection(offsetPatcher);
1629 33 : Create_CADRG_ColormapSection(offsetPatcher, poSrcDS,
1630 33 : priv->bHasTransparentPixels, priv->codebook,
1631 : copyContext);
1632 33 : Create_CADRG_ColorConverterSubsection(
1633 33 : offsetPatcher, priv->bHasTransparentPixels, copyContext);
1634 33 : Create_CADRG_ImageDescriptionSection(offsetPatcher, poSrcDS,
1635 33 : priv->bHasTransparentPixels);
1636 33 : return std::make_unique<CADRGInformation>(std::move(priv));
1637 : }
1638 :
1639 : /************************************************************************/
1640 : /* RPFFrameWriteCADRG_RPFIMG() */
1641 : /************************************************************************/
1642 :
1643 32 : bool RPFFrameWriteCADRG_RPFIMG(GDALOffsetPatcher::OffsetPatcher *offsetPatcher,
1644 : VSILFILE *fp, int &nUDIDL)
1645 : {
1646 32 : std::vector<GDALOffsetPatcher::OffsetPatcherBuffer *> apoBuffers;
1647 32 : int nContentLength = 0;
1648 384 : for (const char *pszName :
1649 : {"LocationComponent", "CoverageSectionSubheader",
1650 : "ColorGrayscaleSectionSubheader", "ColormapSubsection",
1651 224 : "ColorConverterSubsection", "ImageDescriptionSubheader"})
1652 : {
1653 192 : const auto poBuffer = offsetPatcher->GetBufferFromName(pszName);
1654 192 : CPLAssert(poBuffer);
1655 192 : apoBuffers.push_back(poBuffer);
1656 192 : nContentLength += static_cast<int>(poBuffer->GetBuffer().size());
1657 : }
1658 :
1659 32 : CPLAssert(nContentLength <= 99999);
1660 32 : constexpr const char *pszUDOFL = "000";
1661 32 : const char *pszTREPrefix = CPLSPrintf("RPFIMG%05d", nContentLength);
1662 32 : nUDIDL = static_cast<int>(strlen(pszUDOFL) + strlen(pszTREPrefix) +
1663 : nContentLength);
1664 :
1665 : // UDIDL
1666 32 : bool bOK = fp->Write(CPLSPrintf("%05d", nUDIDL), 1, 5) == 5;
1667 :
1668 : // UDOFL
1669 32 : bOK = bOK && fp->Write(pszUDOFL, 1, strlen(pszUDOFL)) == strlen(pszUDOFL);
1670 :
1671 : // UDID
1672 64 : bOK = bOK && fp->Write(pszTREPrefix, 1, strlen(pszTREPrefix)) ==
1673 32 : strlen(pszTREPrefix);
1674 :
1675 224 : for (auto *poBuffer : apoBuffers)
1676 : {
1677 192 : poBuffer->DeclareBufferWrittenAtPosition(VSIFTellL(fp));
1678 576 : bOK = bOK && VSIFWriteL(poBuffer->GetBuffer().data(), 1,
1679 192 : poBuffer->GetBuffer().size(),
1680 192 : fp) == poBuffer->GetBuffer().size();
1681 : }
1682 :
1683 64 : return bOK;
1684 : }
1685 :
1686 : /************************************************************************/
1687 : /* Write_CADRG_MaskSubsection() */
1688 : /************************************************************************/
1689 :
1690 : static bool
1691 32 : Write_CADRG_MaskSubsection(GDALOffsetPatcher::OffsetPatcher *offsetPatcher,
1692 : VSILFILE *fp, GDALDataset *poSrcDS,
1693 : bool bHasTransparentPixels,
1694 : const std::vector<short> &VQImage)
1695 : {
1696 32 : auto poBuffer = offsetPatcher->CreateBuffer(
1697 : "MaskSubsection", /* bEndiannessIsLittle = */ false);
1698 32 : CPLAssert(poBuffer);
1699 32 : poBuffer->DeclareBufferWrittenAtPosition(fp->Tell());
1700 32 : poBuffer->DeclareOffsetAtCurrentPosition("MASK_SUBSECTION_LOCATION");
1701 32 : poBuffer->AppendUInt16(0); // SUBFRAME_SEQUENCE_RECORD_LENGTH
1702 32 : poBuffer->AppendUInt16(0); // TRANSPARENCY_SEQUENCE_RECORD_LENGTH
1703 32 : if (bHasTransparentPixels)
1704 : {
1705 28 : poBuffer->AppendUInt16(8); // TRANSPARENT_OUTPUT_PIXEL_CODE_LENGTH
1706 : // TRANSPARENT_OUTPUT_PIXEL_CODE
1707 28 : poBuffer->AppendByte(static_cast<GByte>(TRANSPARENT_COLOR_TABLE_ENTRY));
1708 :
1709 28 : const int nWidth = poSrcDS->GetRasterXSize();
1710 28 : const int nHeight = poSrcDS->GetRasterYSize();
1711 28 : CPLAssert((nWidth % BLOCK_SIZE) == 0);
1712 28 : CPLAssert((nHeight % BLOCK_SIZE) == 0);
1713 28 : const int nSubFramesPerRow = nWidth / BLOCK_SIZE;
1714 28 : const int nSubFramesPerCol = nHeight / BLOCK_SIZE;
1715 : static_assert((BLOCK_SIZE % SUBSAMPLING) == 0);
1716 28 : constexpr int SUBFRAME_YSIZE = BLOCK_SIZE / SUBSAMPLING;
1717 28 : constexpr int SUBFRAME_XSIZE = BLOCK_SIZE / SUBSAMPLING;
1718 28 : const int nPixelsPerRow = nSubFramesPerRow * SUBFRAME_XSIZE;
1719 28 : constexpr GByte CODE_WORD_BIT_LENGTH = 12;
1720 : static_assert((1 << CODE_WORD_BIT_LENGTH) == CODEBOOK_MAX_SIZE);
1721 : static_assert(
1722 : ((SUBFRAME_XSIZE * SUBFRAME_YSIZE * CODE_WORD_BIT_LENGTH) % 8) ==
1723 : 0);
1724 28 : constexpr int SIZEOF_SUBFRAME_IN_BYTES =
1725 : (SUBFRAME_XSIZE * SUBFRAME_YSIZE * CODE_WORD_BIT_LENGTH) / 8;
1726 :
1727 28 : CPLAssert(VQImage.size() == static_cast<size_t>(nSubFramesPerRow) *
1728 : nSubFramesPerCol * SUBFRAME_YSIZE *
1729 : SUBFRAME_XSIZE);
1730 :
1731 196 : for (int yBlock = 0, nIdxBlock = 0; yBlock < nSubFramesPerCol; ++yBlock)
1732 : {
1733 168 : int nOffsetBlock = yBlock * SUBFRAME_YSIZE * nPixelsPerRow;
1734 1176 : for (int xBlock = 0; xBlock < nSubFramesPerRow;
1735 1008 : ++xBlock, nOffsetBlock += SUBFRAME_XSIZE, ++nIdxBlock)
1736 : {
1737 1008 : bool bBlockHasTransparentPixels = false;
1738 65520 : for (int ySubBlock = 0; ySubBlock < SUBFRAME_YSIZE; ySubBlock++)
1739 : {
1740 64512 : int nOffset = nOffsetBlock + ySubBlock * nPixelsPerRow;
1741 4193280 : for (int xSubBlock = 0; xSubBlock < SUBFRAME_XSIZE;
1742 : ++xSubBlock, ++nOffset)
1743 : {
1744 4128770 : if (VQImage[nOffset] == TRANSPARENT_CODEBOOK_CODE)
1745 3398150 : bBlockHasTransparentPixels = true;
1746 : }
1747 : }
1748 : // Cf MIL-STD-2411 page 23
1749 1008 : if (!bBlockHasTransparentPixels)
1750 138 : poBuffer->AppendUInt32(UINT32_MAX);
1751 : else
1752 870 : poBuffer->AppendUInt32(nIdxBlock *
1753 : SIZEOF_SUBFRAME_IN_BYTES);
1754 : }
1755 : }
1756 : }
1757 : else
1758 : {
1759 4 : poBuffer->AppendUInt16(0); // TRANSPARENT_OUTPUT_PIXEL_CODE_LENGTH
1760 : }
1761 :
1762 32 : return fp->Write(poBuffer->GetBuffer().data(), 1,
1763 32 : poBuffer->GetBuffer().size()) ==
1764 32 : poBuffer->GetBuffer().size();
1765 : }
1766 :
1767 : /************************************************************************/
1768 : /* Write_CADRG_CompressionSection() */
1769 : /************************************************************************/
1770 :
1771 : constexpr uint16_t NUMBER_OF_COMPRESSION_LOOKUP_OFFSET_RECORDS = 4;
1772 :
1773 : static bool
1774 32 : Write_CADRG_CompressionSection(GDALOffsetPatcher::OffsetPatcher *offsetPatcher,
1775 : VSILFILE *fp)
1776 : {
1777 32 : auto poBuffer = offsetPatcher->CreateBuffer(
1778 : "CompressionSection", /* bEndiannessIsLittle = */ false);
1779 32 : CPLAssert(poBuffer);
1780 32 : poBuffer->DeclareBufferWrittenAtPosition(fp->Tell());
1781 32 : poBuffer->DeclareOffsetAtCurrentPosition("COMPRESSION_SECTION_LOCATION");
1782 32 : poBuffer->AppendUInt16(1); // COMPRESSION_ALGORITHM_ID = VQ
1783 32 : poBuffer->AppendUInt16(NUMBER_OF_COMPRESSION_LOOKUP_OFFSET_RECORDS);
1784 : // NUMBER_OF_COMPRESSION_PARAMETER_OFFSET_RECORDS
1785 32 : poBuffer->AppendUInt16(0);
1786 :
1787 32 : return fp->Write(poBuffer->GetBuffer().data(), 1,
1788 32 : poBuffer->GetBuffer().size()) ==
1789 32 : poBuffer->GetBuffer().size();
1790 : }
1791 :
1792 : /************************************************************************/
1793 : /* Write_CADRG_ImageDisplayParametersSection() */
1794 : /************************************************************************/
1795 :
1796 32 : static bool Write_CADRG_ImageDisplayParametersSection(
1797 : GDALOffsetPatcher::OffsetPatcher *offsetPatcher, VSILFILE *fp)
1798 : {
1799 32 : auto poBuffer = offsetPatcher->CreateBuffer(
1800 : "ImageDisplayParametersSubheader", /* bEndiannessIsLittle = */ false);
1801 32 : CPLAssert(poBuffer);
1802 32 : poBuffer->DeclareBufferWrittenAtPosition(fp->Tell());
1803 32 : poBuffer->DeclareOffsetAtCurrentPosition(
1804 : "IMAGE_DISPLAY_PARAMETERS_SECTION_LOCATION");
1805 32 : poBuffer->AppendUInt32(BLOCK_SIZE / SUBSAMPLING); // NUMBER_OF_IMAGE_ROWS
1806 32 : poBuffer->AppendUInt32(BLOCK_SIZE /
1807 : SUBSAMPLING); // NUMBER_OF_CODES_PER_ROW
1808 32 : constexpr GByte CODE_WORD_BIT_LENGTH = 12;
1809 : static_assert((1 << CODE_WORD_BIT_LENGTH) == CODEBOOK_MAX_SIZE);
1810 32 : poBuffer->AppendByte(CODE_WORD_BIT_LENGTH); // IMAGE_CODE_BIT_LENGTH
1811 :
1812 32 : return fp->Write(poBuffer->GetBuffer().data(), 1,
1813 32 : poBuffer->GetBuffer().size()) ==
1814 32 : poBuffer->GetBuffer().size();
1815 : }
1816 :
1817 : /************************************************************************/
1818 : /* Write_CADRG_CompressionLookupSubSection() */
1819 : /************************************************************************/
1820 :
1821 32 : static bool Write_CADRG_CompressionLookupSubSection(
1822 : GDALOffsetPatcher::OffsetPatcher *offsetPatcher, VSILFILE *fp,
1823 : const std::vector<BucketItem<ColorTableBased4x4Pixels>> &codebook)
1824 : {
1825 32 : auto poBuffer = offsetPatcher->CreateBuffer(
1826 : "CompressionLookupSubsection", /* bEndiannessIsLittle = */ false);
1827 32 : CPLAssert(poBuffer);
1828 32 : poBuffer->DeclareBufferWrittenAtPosition(fp->Tell());
1829 32 : poBuffer->DeclareOffsetAtCurrentPosition("COMPRESSION_LOOKUP_LOCATION");
1830 32 : poBuffer->AppendUInt32(6); // COMPRESSION_LOOKUP_OFFSET_TABLE_OFFSET
1831 : // COMPRESSION_LOOKUP_TABLE_OFFSET_RECORD_LENGTH
1832 32 : poBuffer->AppendUInt16(14);
1833 :
1834 32 : constexpr int OFFSET_OF_FIRST_LOOKUP_TABLE = 62;
1835 32 : constexpr uint16_t NUMBER_OF_VALUES_PER_COMPRESSION_RECORDS =
1836 : static_cast<uint16_t>(SUBSAMPLING);
1837 160 : for (int i = 0; i < NUMBER_OF_COMPRESSION_LOOKUP_OFFSET_RECORDS; ++i)
1838 : {
1839 : // COMPRESSION_LOOKUP_TABLE_ID
1840 128 : poBuffer->AppendUInt16(static_cast<uint16_t>(i + 1));
1841 128 : poBuffer->AppendUInt32(
1842 : CODEBOOK_MAX_SIZE); // NUMBER_OF_COMPRESSION_LOOKUP_RECORDS
1843 128 : poBuffer->AppendUInt16(NUMBER_OF_VALUES_PER_COMPRESSION_RECORDS);
1844 128 : poBuffer->AppendUInt16(8); // COMPRESSION_RECORD_VALUE_BIT_LENGTH
1845 128 : poBuffer->AppendUInt32(OFFSET_OF_FIRST_LOOKUP_TABLE +
1846 : CODEBOOK_MAX_SIZE *
1847 128 : NUMBER_OF_VALUES_PER_COMPRESSION_RECORDS *
1848 : i); // COMPRESSION_LOOKUP_TABLE_OFFSET
1849 : }
1850 :
1851 160 : for (int row = 0; row < SUBSAMPLING; ++row)
1852 : {
1853 128 : int i = 0;
1854 524416 : for (; i < static_cast<int>(codebook.size()); ++i)
1855 : {
1856 2621440 : for (int j = 0; j < SUBSAMPLING; ++j)
1857 : {
1858 2097150 : poBuffer->AppendByte(
1859 2097150 : codebook[i].m_vec.val(row * SUBSAMPLING + j));
1860 : }
1861 : }
1862 128 : for (; i < CODEBOOK_MAX_SIZE; ++i)
1863 : {
1864 0 : poBuffer->AppendUInt32(0);
1865 : }
1866 : }
1867 :
1868 32 : return fp->Write(poBuffer->GetBuffer().data(), 1,
1869 32 : poBuffer->GetBuffer().size()) ==
1870 32 : poBuffer->GetBuffer().size();
1871 : }
1872 :
1873 : /************************************************************************/
1874 : /* Write_CADRG_SpatialDataSubsection() */
1875 : /************************************************************************/
1876 :
1877 32 : static bool Write_CADRG_SpatialDataSubsection(
1878 : GDALOffsetPatcher::OffsetPatcher *offsetPatcher, VSILFILE *fp,
1879 : GDALDataset *poSrcDS, const std::vector<short> &VQImage)
1880 : {
1881 32 : auto poBuffer = offsetPatcher->CreateBuffer(
1882 : "SpatialDataSubsection", /* bEndiannessIsLittle = */ false);
1883 32 : CPLAssert(poBuffer);
1884 32 : poBuffer->DeclareBufferWrittenAtPosition(fp->Tell());
1885 32 : poBuffer->DeclareOffsetAtCurrentPosition(
1886 : "SPATIAL_DATA_SUBSECTION_LOCATION");
1887 :
1888 32 : const int nWidth = poSrcDS->GetRasterXSize();
1889 32 : const int nHeight = poSrcDS->GetRasterYSize();
1890 32 : CPLAssert((nWidth % BLOCK_SIZE) == 0);
1891 32 : CPLAssert((nHeight % BLOCK_SIZE) == 0);
1892 32 : const int nSubFramesPerRow = nWidth / BLOCK_SIZE;
1893 32 : const int nSubFramesPerCol = nHeight / BLOCK_SIZE;
1894 : static_assert((BLOCK_SIZE % SUBSAMPLING) == 0);
1895 32 : constexpr int SUBFRAME_YSIZE = BLOCK_SIZE / SUBSAMPLING;
1896 32 : constexpr int SUBFRAME_XSIZE = BLOCK_SIZE / SUBSAMPLING;
1897 32 : const int nPixelsPerRow = nSubFramesPerRow * SUBFRAME_XSIZE;
1898 224 : for (int yBlock = 0; yBlock < nSubFramesPerCol; ++yBlock)
1899 : {
1900 192 : int nOffsetBlock = yBlock * SUBFRAME_YSIZE * nPixelsPerRow;
1901 1344 : for (int xBlock = 0; xBlock < nSubFramesPerRow;
1902 1152 : ++xBlock, nOffsetBlock += SUBFRAME_XSIZE)
1903 : {
1904 74880 : for (int ySubBlock = 0; ySubBlock < SUBFRAME_YSIZE; ySubBlock++)
1905 : {
1906 73728 : int nOffset = nOffsetBlock + ySubBlock * nPixelsPerRow;
1907 : // Combine 2 codes of 12 bits each into 3 bytes
1908 : // This is the reverse of function NITFUncompressVQTile()
1909 2433020 : for (int xSubBlock = 0; xSubBlock < SUBFRAME_XSIZE;
1910 2359300 : xSubBlock += 2, nOffset += 2)
1911 : {
1912 2359300 : const int v1 = VQImage[nOffset + 0];
1913 2359300 : const int v2 = VQImage[nOffset + 1];
1914 2359300 : poBuffer->AppendByte(static_cast<GByte>(v1 >> 4));
1915 2359300 : poBuffer->AppendByte(
1916 2359300 : static_cast<GByte>(((v1 & 0xF) << 4) | (v2 >> 8)));
1917 2359300 : poBuffer->AppendByte(static_cast<GByte>(v2 & 0xFF));
1918 : }
1919 : }
1920 : }
1921 : }
1922 :
1923 32 : return fp->Write(poBuffer->GetBuffer().data(), 1,
1924 32 : poBuffer->GetBuffer().size()) ==
1925 32 : poBuffer->GetBuffer().size();
1926 : }
1927 :
1928 : /************************************************************************/
1929 : /* RPFFrameWriteCADRG_ImageContent() */
1930 : /************************************************************************/
1931 :
1932 32 : bool RPFFrameWriteCADRG_ImageContent(
1933 : GDALOffsetPatcher::OffsetPatcher *offsetPatcher, VSILFILE *fp,
1934 : GDALDataset *poSrcDS, CADRGInformation *info)
1935 : {
1936 32 : return fp->Seek(0, SEEK_END) == 0 &&
1937 32 : Write_CADRG_MaskSubsection(offsetPatcher, fp, poSrcDS,
1938 32 : info->m_private->bHasTransparentPixels,
1939 32 : info->m_private->VQImage) &&
1940 32 : Write_CADRG_CompressionSection(offsetPatcher, fp) &&
1941 32 : Write_CADRG_ImageDisplayParametersSection(offsetPatcher, fp) &&
1942 32 : Write_CADRG_CompressionLookupSubSection(offsetPatcher, fp,
1943 96 : info->m_private->codebook) &&
1944 32 : Write_CADRG_SpatialDataSubsection(offsetPatcher, fp, poSrcDS,
1945 64 : info->m_private->VQImage);
1946 : }
1947 :
1948 : /************************************************************************/
1949 : /* RPFAttribute */
1950 : /************************************************************************/
1951 :
1952 : namespace
1953 : {
1954 : struct RPFAttribute
1955 : {
1956 : uint16_t nAttrId = 0;
1957 : uint8_t nParamId = 0;
1958 : std::string osValue{};
1959 :
1960 130 : RPFAttribute(int nAttrIdIn, int nParamIdIn, const std::string &osValueIn)
1961 130 : : nAttrId(static_cast<uint16_t>(nAttrIdIn)),
1962 130 : nParamId(static_cast<uint8_t>(nParamIdIn)), osValue(osValueIn)
1963 : {
1964 130 : }
1965 : };
1966 : } // namespace
1967 :
1968 : /************************************************************************/
1969 : /* RPFFrameWriteGetDESHeader() */
1970 : /************************************************************************/
1971 :
1972 51 : const char *RPFFrameWriteGetDESHeader()
1973 : {
1974 :
1975 51 : constexpr const char *pszDESHeader =
1976 : "DE" // Segment type
1977 : "Registered Extensions " // DESID
1978 : "01" // DESVER
1979 : "U" // DECLAS
1980 : " " // DESCLSY
1981 : " " // DESCODE
1982 : " " // DESCTLH
1983 : " " // DESREL
1984 : " " // DESDCDT
1985 : " " // DESDCDT
1986 : " " // DESDCXM
1987 : " " // DESDG
1988 : " " // DESDGDT
1989 : " " // DESCLTX
1990 : " " // DESCATP
1991 : " " // DESCAUT
1992 : " " // DESCRSN
1993 : " " // DESSRDT
1994 : " " // DESCTLN
1995 : "UDID " // DESOVFL
1996 : "001" // DESITEM
1997 : "0000" // DESSHL
1998 : ;
1999 :
2000 51 : return pszDESHeader;
2001 : }
2002 :
2003 : /************************************************************************/
2004 : /* RPFFrameWriteCADRG_RPFDES() */
2005 : /************************************************************************/
2006 :
2007 32 : bool RPFFrameWriteCADRG_RPFDES(GDALOffsetPatcher::OffsetPatcher *offsetPatcher,
2008 : VSILFILE *fp, vsi_l_offset nOffsetLDSH,
2009 : const CPLStringList &aosOptions,
2010 : int nReciprocalScale)
2011 : {
2012 32 : bool bOK = fp->Seek(0, SEEK_END) == 0;
2013 :
2014 32 : const char *pszDESHeader = RPFFrameWriteGetDESHeader();
2015 32 : bOK &= fp->Write(pszDESHeader, 1, strlen(pszDESHeader)) ==
2016 32 : strlen(pszDESHeader);
2017 :
2018 64 : std::string osDESData("RPFDES");
2019 64 : std::string osDESDataPayload;
2020 : const auto nPosAttributeSectionSubheader =
2021 32 : fp->Tell() + osDESData.size() + strlen("XXXXX");
2022 :
2023 32 : auto poBufferASSH = offsetPatcher->CreateBuffer(
2024 : "AttributeSectionSubheader", /* bEndiannessIsLittle = */ false);
2025 32 : CPLAssert(poBufferASSH);
2026 32 : poBufferASSH->DeclareBufferWrittenAtPosition(nPosAttributeSectionSubheader);
2027 32 : poBufferASSH->DeclareOffsetAtCurrentPosition(
2028 : "ATTRIBUTE_SECTION_SUBHEADER_LOCATION");
2029 :
2030 32 : std::vector<RPFAttribute> asAttributes;
2031 :
2032 96 : const auto GetYYYMMDDDate = [](const std::string &osValue) -> std::string
2033 : {
2034 96 : if (EQUAL(osValue.c_str(), "NOW"))
2035 : {
2036 : time_t unixTime;
2037 0 : time(&unixTime);
2038 : struct tm brokenDownTime;
2039 0 : CPLUnixTimeToYMDHMS(unixTime, &brokenDownTime);
2040 0 : return CPLString().Printf(
2041 0 : "%04d%02d%02d", brokenDownTime.tm_year + 1900,
2042 0 : brokenDownTime.tm_mon + 1, brokenDownTime.tm_mday);
2043 : }
2044 : else
2045 : {
2046 96 : return osValue;
2047 : }
2048 : };
2049 :
2050 : {
2051 : const char *pszV =
2052 32 : aosOptions.FetchNameValueDef("CURRENCY_DATE", "20260101");
2053 32 : if (pszV && pszV[0])
2054 : {
2055 0 : asAttributes.emplace_back(1, 1,
2056 32 : StrPadTruncate(GetYYYMMDDDate(pszV), 8));
2057 : }
2058 : }
2059 :
2060 : {
2061 : const char *pszV =
2062 32 : aosOptions.FetchNameValueDef("PRODUCTION_DATE", "20260101");
2063 32 : if (pszV && pszV[0])
2064 : {
2065 0 : asAttributes.emplace_back(2, 1,
2066 32 : StrPadTruncate(GetYYYMMDDDate(pszV), 8));
2067 : }
2068 : }
2069 :
2070 : {
2071 : const char *pszV =
2072 32 : aosOptions.FetchNameValueDef("SIGNIFICANT_DATE", "20260101");
2073 32 : if (pszV && pszV[0])
2074 : {
2075 0 : asAttributes.emplace_back(3, 1,
2076 32 : StrPadTruncate(GetYYYMMDDDate(pszV), 8));
2077 : }
2078 : }
2079 :
2080 32 : if (const char *pszV = aosOptions.FetchNameValue("DATA_SERIES_DESIGNATION"))
2081 : {
2082 0 : asAttributes.emplace_back(4, 1, StrPadTruncate(pszV, 10));
2083 : }
2084 32 : else if (const char *pszSeriesCode =
2085 32 : aosOptions.FetchNameValue("SERIES_CODE"))
2086 : {
2087 19 : const auto *psEntry = NITFGetRPFSeriesInfoFromCode(pszSeriesCode);
2088 19 : if (psEntry && psEntry->abbreviation[0])
2089 : {
2090 : // If the data series abbreviation doesn't contain a scale indication,
2091 : // add it.
2092 2 : std::string osVal(psEntry->abbreviation);
2093 2 : if (osVal.find('0') != std::string::npos)
2094 : {
2095 1 : if (nReciprocalScale >= Million)
2096 0 : osVal += CPLSPrintf(" %dM", nReciprocalScale / Million);
2097 1 : else if (nReciprocalScale >= Kilo)
2098 1 : osVal += CPLSPrintf(" %dK", nReciprocalScale / Kilo);
2099 : }
2100 :
2101 2 : asAttributes.emplace_back(4, 1, StrPadTruncate(osVal, 10));
2102 : }
2103 : }
2104 :
2105 32 : if (const char *pszV = aosOptions.FetchNameValue("MAP_DESIGNATION"))
2106 : {
2107 0 : asAttributes.emplace_back(4, 2, StrPadTruncate(pszV, 8));
2108 : }
2109 :
2110 : // Horizontal datum code
2111 32 : asAttributes.emplace_back(7, 1, StrPadTruncate("WGE", 4));
2112 :
2113 32 : const auto nAttrCount = static_cast<uint16_t>(asAttributes.size());
2114 32 : poBufferASSH->AppendUInt16(nAttrCount);
2115 32 : poBufferASSH->AppendUInt16(0); // NUMBER_OF_EXPLICIT_AREAL_COVERAGE_RECORDS
2116 32 : poBufferASSH->AppendUInt32(0); // ATTRIBUTE_OFFSET_TABLE_OFFSET
2117 32 : constexpr uint16_t ATTRIBUTE_OFFSET_RECORD_LENGTH =
2118 : static_cast<uint16_t>(8);
2119 32 : poBufferASSH->AppendUInt16(ATTRIBUTE_OFFSET_RECORD_LENGTH);
2120 :
2121 : osDESDataPayload.insert(
2122 32 : osDESDataPayload.end(),
2123 32 : reinterpret_cast<const char *>(poBufferASSH->GetBuffer().data()),
2124 64 : reinterpret_cast<const char *>(poBufferASSH->GetBuffer().data() +
2125 64 : poBufferASSH->GetBuffer().size()));
2126 :
2127 32 : auto poBufferAS = offsetPatcher->CreateBuffer(
2128 : "AttributeSubsection", /* bEndiannessIsLittle = */ false);
2129 32 : CPLAssert(poBufferAS);
2130 32 : poBufferAS->DeclareBufferWrittenAtPosition(
2131 32 : nPosAttributeSectionSubheader + poBufferASSH->GetBuffer().size());
2132 32 : poBufferAS->DeclareOffsetAtCurrentPosition("ATTRIBUTE_SUBSECTION_LOCATION");
2133 :
2134 : size_t nAttrValueOffset =
2135 32 : ATTRIBUTE_OFFSET_RECORD_LENGTH * asAttributes.size();
2136 :
2137 : // Attribute definitions
2138 162 : for (const auto &sAttr : asAttributes)
2139 : {
2140 130 : poBufferAS->AppendUInt16(sAttr.nAttrId);
2141 130 : poBufferAS->AppendByte(sAttr.nParamId);
2142 130 : poBufferAS->AppendByte(0); // Areal coverage sequence number
2143 130 : poBufferAS->AppendUInt32(static_cast<uint32_t>(
2144 : nAttrValueOffset)); // Attribute record offset
2145 130 : nAttrValueOffset += sAttr.osValue.size();
2146 : }
2147 :
2148 : // Attribute values
2149 162 : for (const auto &sAttr : asAttributes)
2150 : {
2151 130 : poBufferAS->AppendString(sAttr.osValue);
2152 : }
2153 :
2154 : osDESDataPayload.insert(
2155 32 : osDESDataPayload.end(),
2156 32 : reinterpret_cast<const char *>(poBufferAS->GetBuffer().data()),
2157 64 : reinterpret_cast<const char *>(poBufferAS->GetBuffer().data() +
2158 64 : poBufferAS->GetBuffer().size()));
2159 :
2160 32 : CPLAssert(osDESDataPayload.size() <= 99999U);
2161 32 : osDESData += CPLSPrintf("%05d", static_cast<int>(osDESDataPayload.size()));
2162 32 : osDESData += osDESDataPayload;
2163 32 : bOK &=
2164 32 : fp->Write(osDESData.c_str(), 1, osDESData.size()) == osDESData.size();
2165 :
2166 : // Update LDSH and LD in the NITF Header
2167 32 : const int iDES = 0;
2168 32 : bOK &= fp->Seek(nOffsetLDSH + iDES * 13, SEEK_SET) == 0;
2169 32 : bOK &= fp->Write(CPLSPrintf("%04d", static_cast<int>(strlen(pszDESHeader))),
2170 32 : 1, 4) == 4;
2171 32 : bOK &= fp->Write(CPLSPrintf("%09d", static_cast<int>(osDESData.size())), 1,
2172 32 : 9) == 9;
2173 :
2174 64 : return bOK;
2175 : }
2176 :
2177 : /************************************************************************/
2178 : /* CADRGGetWarpedVRTDataset() */
2179 : /************************************************************************/
2180 :
2181 : static std::unique_ptr<GDALDataset>
2182 23 : CADRGGetWarpedVRTDataset(GDALDataset *poSrcDS, int nZone, double dfXMin,
2183 : double dfYMin, double dfXMax, double dfYMax,
2184 : double dfResX, double dfResY,
2185 : const char *pszResampling)
2186 : {
2187 46 : CPLStringList aosWarpArgs;
2188 23 : aosWarpArgs.push_back("-of");
2189 23 : aosWarpArgs.push_back("VRT");
2190 23 : aosWarpArgs.push_back("-t_srs");
2191 23 : if (nZone == MAX_ZONE_NORTHERN_HEMISPHERE)
2192 1 : aosWarpArgs.push_back(pszNorthPolarProjection);
2193 22 : else if (nZone == MAX_ZONE)
2194 1 : aosWarpArgs.push_back(pszSouthPolarProjection);
2195 : else
2196 21 : aosWarpArgs.push_back("EPSG:4326");
2197 23 : aosWarpArgs.push_back("-te");
2198 23 : aosWarpArgs.push_back(CPLSPrintf("%.17g", dfXMin));
2199 23 : aosWarpArgs.push_back(CPLSPrintf("%.17g", dfYMin));
2200 23 : aosWarpArgs.push_back(CPLSPrintf("%.17g", dfXMax));
2201 23 : aosWarpArgs.push_back(CPLSPrintf("%.17g", dfYMax));
2202 23 : aosWarpArgs.push_back("-tr");
2203 23 : aosWarpArgs.push_back(CPLSPrintf("%.17g", dfResX));
2204 23 : aosWarpArgs.push_back(CPLSPrintf("%.17g", dfResY));
2205 23 : if (poSrcDS->GetRasterBand(1)->GetColorTable())
2206 : {
2207 4 : aosWarpArgs.push_back("-dstnodata");
2208 4 : aosWarpArgs.push_back(CPLSPrintf("%d", TRANSPARENT_COLOR_TABLE_ENTRY));
2209 : }
2210 : else
2211 : {
2212 19 : aosWarpArgs.push_back("-r");
2213 19 : aosWarpArgs.push_back(CPLString(pszResampling).tolower());
2214 19 : if (poSrcDS->GetRasterCount() == 3)
2215 18 : aosWarpArgs.push_back("-dstalpha");
2216 : }
2217 : std::unique_ptr<GDALWarpAppOptions, decltype(&GDALWarpAppOptionsFree)>
2218 : psWarpOptions(GDALWarpAppOptionsNew(aosWarpArgs.List(), nullptr),
2219 46 : GDALWarpAppOptionsFree);
2220 23 : std::unique_ptr<GDALDataset> poWarpedDS;
2221 23 : if (psWarpOptions)
2222 : {
2223 23 : GDALDatasetH hSrcDS = GDALDataset::ToHandle(poSrcDS);
2224 23 : poWarpedDS.reset(GDALDataset::FromHandle(
2225 23 : GDALWarp("", nullptr, 1, &hSrcDS, psWarpOptions.get(), nullptr)));
2226 : }
2227 46 : return poWarpedDS;
2228 : }
2229 :
2230 : /************************************************************************/
2231 : /* CADRGGetClippedDataset() */
2232 : /************************************************************************/
2233 :
2234 : static std::unique_ptr<GDALDataset>
2235 31 : CADRGGetClippedDataset(GDALDataset *poSrcDS, double dfXMin, double dfYMin,
2236 : double dfXMax, double dfYMax)
2237 : {
2238 62 : CPLStringList aosTranslateArgs;
2239 31 : aosTranslateArgs.push_back("-of");
2240 31 : aosTranslateArgs.push_back("MEM");
2241 31 : aosTranslateArgs.push_back("-projwin");
2242 31 : aosTranslateArgs.push_back(CPLSPrintf("%.17g", dfXMin));
2243 31 : aosTranslateArgs.push_back(CPLSPrintf("%.17g", dfYMax));
2244 31 : aosTranslateArgs.push_back(CPLSPrintf("%.17g", dfXMax));
2245 31 : aosTranslateArgs.push_back(CPLSPrintf("%.17g", dfYMin));
2246 : std::unique_ptr<GDALTranslateOptions, decltype(&GDALTranslateOptionsFree)>
2247 : psTranslateOptions(
2248 : GDALTranslateOptionsNew(aosTranslateArgs.List(), nullptr),
2249 62 : GDALTranslateOptionsFree);
2250 31 : std::unique_ptr<GDALDataset> poClippedDS;
2251 31 : if (psTranslateOptions)
2252 : {
2253 31 : poClippedDS.reset(GDALDataset::FromHandle(
2254 : GDALTranslate("", GDALDataset::ToHandle(poSrcDS),
2255 31 : psTranslateOptions.get(), nullptr)));
2256 : }
2257 62 : return poClippedDS;
2258 : }
2259 :
2260 : /************************************************************************/
2261 : /* CADRGGetPalettedDataset() */
2262 : /************************************************************************/
2263 :
2264 : static std::unique_ptr<GDALDataset>
2265 21 : CADRGGetPalettedDataset(GDALDataset *poSrcDS, GDALColorTable *poCT,
2266 : int nColorQuantizationBits)
2267 : {
2268 21 : CPLAssert(poSrcDS->GetRasterCount() == 3 || poSrcDS->GetRasterCount() == 4);
2269 21 : auto poMemDrv = GetGDALDriverManager()->GetDriverByName("MEM");
2270 : std::unique_ptr<GDALDataset> poPalettedDS(
2271 : poMemDrv->Create("", poSrcDS->GetRasterXSize(),
2272 21 : poSrcDS->GetRasterYSize(), 1, GDT_UInt8, nullptr));
2273 21 : if (poPalettedDS)
2274 : {
2275 21 : poPalettedDS->SetSpatialRef(poSrcDS->GetSpatialRef());
2276 21 : GDALGeoTransform gt;
2277 21 : if (poSrcDS->GetGeoTransform(gt) == CE_None)
2278 21 : poPalettedDS->SetGeoTransform(gt);
2279 21 : poPalettedDS->GetRasterBand(1)->SetColorTable(poCT);
2280 21 : poPalettedDS->GetRasterBand(1)->SetNoDataValue(
2281 21 : TRANSPARENT_COLOR_TABLE_ENTRY);
2282 21 : GDALDitherRGB2PCTInternal(
2283 : GDALRasterBand::ToHandle(poSrcDS->GetRasterBand(1)),
2284 : GDALRasterBand::ToHandle(poSrcDS->GetRasterBand(2)),
2285 : GDALRasterBand::ToHandle(poSrcDS->GetRasterBand(3)),
2286 : GDALRasterBand::ToHandle(poPalettedDS->GetRasterBand(1)),
2287 : GDALColorTable::ToHandle(poCT), nColorQuantizationBits, nullptr,
2288 : /* dither = */ false, nullptr, nullptr);
2289 : }
2290 21 : return poPalettedDS;
2291 : }
2292 :
2293 : /************************************************************************/
2294 : /* CADRG_RGB_Type */
2295 : /************************************************************************/
2296 :
2297 : struct CADRG_RGB_Type
2298 : {
2299 : };
2300 :
2301 : /************************************************************************/
2302 : /* Vector<CADRG_RGB_Type> */
2303 : /************************************************************************/
2304 :
2305 : template <> class Vector<CADRG_RGB_Type>
2306 : {
2307 : private:
2308 : GByte m_R = 0;
2309 : GByte m_G = 0;
2310 : GByte m_B = 0;
2311 :
2312 : Vector() = default;
2313 :
2314 : public:
2315 1601 : explicit Vector(GByte R, GByte G, GByte B) : m_R(R), m_G(G), m_B(B)
2316 : {
2317 1601 : }
2318 :
2319 : static constexpr int DIM_COUNT /* specialize */ = 3;
2320 :
2321 : static constexpr bool getReturnUInt8 /* specialize */ = true;
2322 :
2323 : static constexpr bool hasComputeFourSquaredDistances = false;
2324 :
2325 117510 : inline int get(int i, const CADRG_RGB_Type &) const /* specialize */
2326 : {
2327 117510 : return i == 0 ? m_R : i == 1 ? m_G : m_B;
2328 : }
2329 :
2330 671 : inline GByte r() const
2331 : {
2332 671 : return m_R;
2333 : }
2334 :
2335 671 : inline GByte g() const
2336 : {
2337 671 : return m_G;
2338 : }
2339 :
2340 671 : inline GByte b() const
2341 : {
2342 671 : return m_B;
2343 : }
2344 :
2345 : /************************************************************************/
2346 : /* squared_distance() */
2347 : /************************************************************************/
2348 :
2349 11437 : int squared_distance(const Vector &other,
2350 : const CADRG_RGB_Type &) const /* specialize */
2351 : {
2352 11437 : const int nSqDist1 = square(m_R - other.m_R);
2353 11437 : const int nSqDist2 = square(m_G - other.m_G);
2354 11437 : const int nSqDist3 = square(m_B - other.m_B);
2355 11437 : return nSqDist1 + nSqDist2 + nSqDist3;
2356 : }
2357 :
2358 : /************************************************************************/
2359 : /* centroid() */
2360 : /************************************************************************/
2361 :
2362 1305 : static Vector centroid(const Vector &a, int nA, const Vector &b, int nB,
2363 : const CADRG_RGB_Type &) /* specialize */
2364 : {
2365 1305 : Vector res;
2366 1305 : res.m_R = static_cast<GByte>((static_cast<uint64_t>(a.m_R) * nA +
2367 1305 : static_cast<uint64_t>(b.m_R) * nB +
2368 1305 : (nA + nB) / 2) /
2369 1305 : (nA + nB));
2370 1305 : res.m_G = static_cast<GByte>((static_cast<uint64_t>(a.m_G) * nA +
2371 1305 : static_cast<uint64_t>(b.m_G) * nB +
2372 1305 : (nA + nB) / 2) /
2373 1305 : (nA + nB));
2374 1305 : res.m_B = static_cast<GByte>((static_cast<uint64_t>(a.m_B) * nA +
2375 1305 : static_cast<uint64_t>(b.m_B) * nB +
2376 1305 : (nA + nB) / 2) /
2377 1305 : (nA + nB));
2378 1305 : return res;
2379 : }
2380 :
2381 : /************************************************************************/
2382 : /* operator == () */
2383 : /************************************************************************/
2384 :
2385 6988 : inline bool operator==(const Vector &other) const
2386 : {
2387 6988 : return m_R == other.m_R && m_G == other.m_G && m_B == other.m_B;
2388 : }
2389 :
2390 : /************************************************************************/
2391 : /* operator < () */
2392 : /************************************************************************/
2393 :
2394 : // Purely arbitrary for the purpose of distinguishing a vector from
2395 : // another one
2396 74378 : inline bool operator<(const Vector &other) const
2397 : {
2398 74378 : const int nA = (m_R) | (m_G << 8) | (m_B << 16);
2399 74378 : const int nB = (m_R) | (other.m_G << 8) | (other.m_B << 16);
2400 74378 : return nA < nB;
2401 : }
2402 : };
2403 :
2404 : /************************************************************************/
2405 : /* ComputeColorTables() */
2406 : /************************************************************************/
2407 :
2408 26 : static bool ComputeColorTables(GDALDataset *poSrcDS, GDALColorTable &oCT,
2409 : int nColorQuantizationBits,
2410 : GDALProgressFunc pfnProgress,
2411 : void *pProgressData,
2412 : CADRGCreateCopyContext ©Context)
2413 :
2414 : {
2415 52 : std::vector<GUIntBig> anPixelCountPerColorTableEntry;
2416 26 : if (poSrcDS->GetRasterCount() >= 3)
2417 : {
2418 21 : if (GDALComputeMedianCutPCT(
2419 : GDALRasterBand::ToHandle(poSrcDS->GetRasterBand(1)),
2420 : GDALRasterBand::ToHandle(poSrcDS->GetRasterBand(2)),
2421 : GDALRasterBand::ToHandle(poSrcDS->GetRasterBand(3)), nullptr,
2422 : nullptr, nullptr, nullptr, CADRG_MAX_COLOR_ENTRY_COUNT,
2423 : nColorQuantizationBits, static_cast<GUIntBig *>(nullptr),
2424 : GDALColorTable::ToHandle(&oCT), pfnProgress, pProgressData,
2425 21 : &anPixelCountPerColorTableEntry) != CE_None)
2426 : {
2427 0 : return false;
2428 : }
2429 : }
2430 : else
2431 : {
2432 : // Compute histogram of the source dataset
2433 5 : const auto poCT = poSrcDS->GetRasterBand(1)->GetColorTable();
2434 5 : CPLAssert(poCT);
2435 : const int nColors =
2436 5 : std::min(poCT->GetColorEntryCount(), CADRG_MAX_COLOR_ENTRY_COUNT);
2437 870 : for (int i = 0; i < nColors; ++i)
2438 : {
2439 865 : oCT.SetColorEntry(i, poCT->GetColorEntry(i));
2440 : }
2441 5 : anPixelCountPerColorTableEntry.resize(nColors);
2442 5 : const int nXSize = poSrcDS->GetRasterXSize();
2443 5 : const int nYSize = poSrcDS->GetRasterYSize();
2444 5 : std::vector<GByte> abyLine(nXSize);
2445 6164 : for (int iY = 0; iY < nYSize; ++iY)
2446 : {
2447 12318 : if (poSrcDS->GetRasterBand(1)->RasterIO(
2448 6159 : GF_Read, 0, iY, nXSize, 1, abyLine.data(), nXSize, 1,
2449 6159 : GDT_UInt8, 0, 0, nullptr) != CE_None)
2450 : {
2451 0 : return false;
2452 : }
2453 9443370 : for (int iX = 0; iX < nXSize; ++iX)
2454 : {
2455 9437210 : const int nVal = abyLine[iX];
2456 9437210 : if (nVal < nColors)
2457 7309470 : ++anPixelCountPerColorTableEntry[nVal];
2458 : }
2459 : }
2460 : }
2461 :
2462 52 : std::vector<BucketItem<CADRG_RGB_Type>> vectors;
2463 : const int nColors =
2464 26 : std::min(oCT.GetColorEntryCount(), CADRG_MAX_COLOR_ENTRY_COUNT);
2465 26 : CPLAssert(anPixelCountPerColorTableEntry.size() >=
2466 : static_cast<size_t>(nColors));
2467 :
2468 26 : GUIntBig nTotalCount = 0;
2469 1698 : for (int i = 0; i < nColors; ++i)
2470 : {
2471 1672 : nTotalCount += anPixelCountPerColorTableEntry[i];
2472 : }
2473 :
2474 : // 3 for R,G,B
2475 52 : std::map<std::array<GByte, 3>, std::vector<int>> oMapUniqueEntries;
2476 1698 : for (int i = 0; i < nColors; ++i)
2477 : {
2478 1672 : const auto entry = oCT.GetColorEntry(i);
2479 1672 : if (entry->c4 && anPixelCountPerColorTableEntry[i])
2480 : {
2481 1601 : const auto R = static_cast<GByte>(entry->c1);
2482 1601 : const auto G = static_cast<GByte>(entry->c2);
2483 1601 : const auto B = static_cast<GByte>(entry->c3);
2484 1601 : oMapUniqueEntries[std::array<GByte, 3>{R, G, B}].push_back(i);
2485 : }
2486 : }
2487 :
2488 26 : const int nUniqueColors = static_cast<int>(oMapUniqueEntries.size());
2489 1627 : for (auto &[RGB, indices] : oMapUniqueEntries)
2490 : {
2491 1601 : uint64_t nThisEntryCount = 0;
2492 3202 : for (int idx : indices)
2493 1601 : nThisEntryCount += anPixelCountPerColorTableEntry[idx];
2494 :
2495 : // Rescale pixel counts for primary color table so that their sum
2496 : // does not exceed INT_MAX, as this is the type of m_count in the
2497 : // BucketItem class.
2498 1601 : const int nCountRescaled = std::max(
2499 3202 : 1, static_cast<int>(
2500 3202 : static_cast<double>(nThisEntryCount) /
2501 1601 : static_cast<double>(std::max<GUIntBig>(1, nTotalCount)) *
2502 1601 : (INT_MAX - nUniqueColors)));
2503 1601 : Vector<CADRG_RGB_Type> v(RGB[0], RGB[1], RGB[2]);
2504 1601 : vectors.emplace_back(v, nCountRescaled, std::move(indices));
2505 : }
2506 :
2507 : // Create the KD-Tree
2508 52 : PNNKDTree<CADRG_RGB_Type> kdtree;
2509 : CADRG_RGB_Type ctxt;
2510 :
2511 26 : int nCodeCount = kdtree.insert(std::move(vectors), ctxt);
2512 :
2513 : // Compute 32 entry color table
2514 26 : if (nCodeCount > CADRG_SECOND_CT_COUNT)
2515 : {
2516 7 : nCodeCount = kdtree.cluster(nCodeCount, CADRG_SECOND_CT_COUNT, ctxt);
2517 7 : if (nCodeCount == 0)
2518 0 : return false;
2519 : }
2520 :
2521 26 : copyContext.oCT2 = GDALColorTable();
2522 26 : copyContext.anMapCT1ToCT2.clear();
2523 26 : copyContext.anMapCT1ToCT2.resize(CADRG_MAX_COLOR_ENTRY_COUNT);
2524 26 : kdtree.iterateOverLeaves(
2525 2178 : [&oCT, ©Context](PNNKDTree<CADRG_RGB_Type> &node)
2526 : {
2527 87 : CPL_IGNORE_RET_VAL(oCT);
2528 87 : int i = copyContext.oCT2.GetColorEntryCount();
2529 490 : for (auto &item : node.bucketItems())
2530 : {
2531 2004 : for (const auto idx : item.m_origVectorIndices)
2532 : {
2533 : #ifdef DEBUG_VERBOSE
2534 : const auto psSrcEntry = oCT.GetColorEntry(idx);
2535 : CPLDebugOnly(
2536 : "CADRG", "Second CT: %d: %d %d %d <-- %d: %d %d %d", i,
2537 : item.m_vec.r(), item.m_vec.g(), item.m_vec.b(), idx,
2538 : psSrcEntry->c1, psSrcEntry->c2, psSrcEntry->c3);
2539 : #endif
2540 1601 : copyContext.anMapCT1ToCT2[idx] = i;
2541 : }
2542 : GDALColorEntry sEntry;
2543 403 : sEntry.c1 = item.m_vec.r();
2544 403 : sEntry.c2 = item.m_vec.g();
2545 403 : sEntry.c3 = item.m_vec.b();
2546 403 : sEntry.c4 = 255;
2547 403 : copyContext.oCT2.SetColorEntry(i, &sEntry);
2548 403 : ++i;
2549 : }
2550 87 : });
2551 :
2552 : // Compute 16 entry color table
2553 26 : if (nCodeCount > CADRG_THIRD_CT_COUNT)
2554 : {
2555 15 : nCodeCount = kdtree.cluster(nCodeCount, CADRG_THIRD_CT_COUNT, ctxt);
2556 15 : if (nCodeCount == 0)
2557 0 : return false;
2558 : }
2559 :
2560 26 : copyContext.oCT3 = GDALColorTable();
2561 26 : copyContext.anMapCT1ToCT3.clear();
2562 26 : copyContext.anMapCT1ToCT3.resize(CADRG_MAX_COLOR_ENTRY_COUNT);
2563 26 : kdtree.iterateOverLeaves(
2564 2007 : [&oCT, ©Context](PNNKDTree<CADRG_RGB_Type> &node)
2565 : {
2566 69 : CPL_IGNORE_RET_VAL(oCT);
2567 69 : int i = copyContext.oCT3.GetColorEntryCount();
2568 337 : for (auto &item : node.bucketItems())
2569 : {
2570 1869 : for (const auto idx : item.m_origVectorIndices)
2571 : {
2572 : #ifdef DEBUG_VERBOSE
2573 : const auto psSrcEntry = oCT.GetColorEntry(idx);
2574 : CPLDebugOnly(
2575 : "CADRG", "Third CT: %d: %d %d %d <-- %d: %d %d %d", i,
2576 : item.m_vec.r(), item.m_vec.g(), item.m_vec.b(), idx,
2577 : psSrcEntry->c1, psSrcEntry->c2, psSrcEntry->c3);
2578 : #endif
2579 1601 : copyContext.anMapCT1ToCT3[idx] = i;
2580 : }
2581 : GDALColorEntry sEntry;
2582 268 : sEntry.c1 = item.m_vec.r();
2583 268 : sEntry.c2 = item.m_vec.g();
2584 268 : sEntry.c3 = item.m_vec.b();
2585 268 : sEntry.c4 = 255;
2586 268 : copyContext.oCT3.SetColorEntry(i, &sEntry);
2587 268 : ++i;
2588 : }
2589 69 : });
2590 :
2591 : // Add transparency entry to primary color table
2592 26 : GDALColorEntry sEntry = {0, 0, 0, 0};
2593 26 : oCT.SetColorEntry(TRANSPARENT_COLOR_TABLE_ENTRY, &sEntry);
2594 :
2595 26 : return true;
2596 : }
2597 :
2598 : /************************************************************************/
2599 : /* CADRGCreateCopy() */
2600 : /************************************************************************/
2601 :
2602 : #ifndef NITFDUMP_BUILD
2603 : std::variant<bool, std::unique_ptr<GDALDataset>>
2604 71 : CADRGCreateCopy(const char *pszFilename, GDALDataset *poSrcDS, int bStrict,
2605 : CSLConstList papszOptions, GDALProgressFunc pfnProgress,
2606 : void *pProgressData, int nRecLevel,
2607 : CADRGCreateCopyContext *copyContext)
2608 : {
2609 71 : int &nReciprocalScale = copyContext->nReciprocalScale;
2610 :
2611 71 : if (poSrcDS->GetRasterBand(1)->GetRasterDataType() != GDT_UInt8)
2612 : {
2613 1 : CPLError(CE_Failure, CPLE_AppDefined,
2614 : "CADRG only supports datasets of UInt8 data type");
2615 1 : return false;
2616 : }
2617 :
2618 70 : const auto poSrcSRS = poSrcDS->GetSpatialRef();
2619 70 : if (!poSrcSRS)
2620 : {
2621 1 : CPLError(CE_Failure, CPLE_AppDefined,
2622 : "CADRG only supports datasets with a CRS");
2623 1 : return false;
2624 : }
2625 :
2626 69 : GDALGeoTransform srcGT;
2627 69 : if (poSrcDS->GetGeoTransform(srcGT) != CE_None)
2628 : {
2629 1 : CPLError(CE_Failure, CPLE_AppDefined,
2630 : "CADRG only supports datasets with a geotransform");
2631 1 : return false;
2632 : }
2633 :
2634 68 : double dfDPIOverride = 0;
2635 68 : const char *pszDPI = CSLFetchNameValue(papszOptions, "DPI");
2636 68 : if (pszDPI)
2637 : {
2638 1 : dfDPIOverride = CPLAtof(pszDPI);
2639 1 : if (!(dfDPIOverride >= 1 && dfDPIOverride <= 7200))
2640 : {
2641 1 : CPLError(CE_Failure, CPLE_AppDefined, "Invalid value for DPI: %s",
2642 : pszDPI);
2643 1 : return false;
2644 : }
2645 : }
2646 :
2647 : const char *pszSeriesCode =
2648 67 : CSLFetchNameValueDef(papszOptions, "SERIES_CODE", "MM");
2649 :
2650 : const std::string osResampling =
2651 134 : CSLFetchNameValueDef(papszOptions, "RESAMPLING", "CUBIC");
2652 :
2653 67 : const char *pszScale = CSLFetchNameValue(papszOptions, "SCALE");
2654 67 : if (pszScale)
2655 : {
2656 51 : if (EQUAL(pszScale, "GUESS"))
2657 : {
2658 0 : bool bGotDPI = false;
2659 0 : nReciprocalScale = RPFGetCADRGClosestReciprocalScale(
2660 : poSrcDS, dfDPIOverride, bGotDPI);
2661 0 : if (nReciprocalScale <= 0)
2662 : {
2663 : // Error message emitted by RPFGetCADRGClosestReciprocalScale()
2664 0 : return false;
2665 : }
2666 : else
2667 : {
2668 0 : CPLDebug("CADRG", "Guessed reciprocal scale: %d",
2669 : nReciprocalScale);
2670 : }
2671 : }
2672 : else
2673 : {
2674 51 : nReciprocalScale = atoi(pszScale);
2675 51 : if (!(nReciprocalScale >= Kilo && nReciprocalScale <= 20 * Million))
2676 : {
2677 1 : CPLError(CE_Failure, CPLE_AppDefined,
2678 : "Invalid value for SCALE: %s", pszScale);
2679 1 : return false;
2680 : }
2681 : }
2682 :
2683 : const int nTheoreticalScale =
2684 50 : RPFCADRGGetScaleFromDataSeriesCode(pszSeriesCode);
2685 50 : if (nTheoreticalScale != 0 && nTheoreticalScale != nReciprocalScale &&
2686 : nRecLevel == 0)
2687 : {
2688 0 : CPLError(CE_Warning, CPLE_AppDefined,
2689 : "Theoretical reciprocal scale from data series code %s is "
2690 : "%d, whereas SCALE has been specified to %d",
2691 : pszSeriesCode, nTheoreticalScale, nReciprocalScale);
2692 : }
2693 : }
2694 : else
2695 : {
2696 16 : nReciprocalScale = RPFCADRGGetScaleFromDataSeriesCode(pszSeriesCode);
2697 : }
2698 :
2699 66 : if (pszDPI && !(pszScale && EQUAL(pszScale, "GUESS")))
2700 : {
2701 0 : CPLError(CE_Warning, CPLE_AppDefined,
2702 : "DPI is ignored when SCALE is not set to GUESS");
2703 : }
2704 :
2705 : const CPLString osExtUpperCase(
2706 198 : CPLString(CPLGetExtensionSafe(pszFilename)).toupper());
2707 66 : bool bLooksLikeCADRGFilename = (osExtUpperCase == "NTF");
2708 99 : if (!bLooksLikeCADRGFilename && osExtUpperCase.size() == 3 &&
2709 33 : RPFCADRGZoneCharToNum(osExtUpperCase[2]) > 0)
2710 : {
2711 : bLooksLikeCADRGFilename =
2712 33 : RPFCADRGIsKnownDataSeriesCode(osExtUpperCase.substr(0, 2).c_str());
2713 : }
2714 :
2715 66 : const bool bColorTablePerFrame = CPLTestBool(
2716 : CSLFetchNameValueDef(papszOptions, "COLOR_TABLE_PER_FRAME", "NO"));
2717 :
2718 106 : if (!(poSrcDS->GetRasterXSize() == CADRG_FRAME_PIXEL_COUNT &&
2719 40 : poSrcDS->GetRasterYSize() == CADRG_FRAME_PIXEL_COUNT))
2720 : {
2721 : VSIStatBufL sStat;
2722 26 : if (VSIStatL(pszFilename, &sStat) == 0)
2723 : {
2724 19 : if (VSI_ISREG(sStat.st_mode))
2725 : {
2726 1 : CPLError(CE_Failure, CPLE_AppDefined,
2727 : "Given that source dataset dimension do not match "
2728 : "a %dx%d frame, several frames will be generated "
2729 : "and thus the output filename should be a "
2730 : "directory name",
2731 : CADRG_FRAME_PIXEL_COUNT, CADRG_FRAME_PIXEL_COUNT);
2732 1 : return false;
2733 : }
2734 : }
2735 : else
2736 : {
2737 7 : if (bLooksLikeCADRGFilename)
2738 : {
2739 1 : CPLError(CE_Failure, CPLE_AppDefined,
2740 : "Given that source dataset dimension do not match "
2741 : "a %dx%d frame, several frames will be "
2742 : "generated and thus the output filename "
2743 : "should be a directory name (without a NITF or "
2744 : "CADRG file extension)",
2745 : CADRG_FRAME_PIXEL_COUNT, CADRG_FRAME_PIXEL_COUNT);
2746 1 : return false;
2747 : }
2748 : }
2749 : }
2750 :
2751 64 : const auto poCT = poSrcDS->GetRasterBand(1)->GetColorTable();
2752 64 : if (poCT && poCT->GetColorEntryCount() > CADRG_MAX_COLOR_ENTRY_COUNT)
2753 : {
2754 149 : for (int i = CADRG_MAX_COLOR_ENTRY_COUNT;
2755 149 : i < poCT->GetColorEntryCount(); ++i)
2756 : {
2757 114 : const auto psEntry = poCT->GetColorEntry(i);
2758 114 : if (psEntry->c1 != 0 || psEntry->c2 != 0 || psEntry->c3 != 0)
2759 : {
2760 1 : CPLError(CE_Failure, CPLE_AppDefined,
2761 : "CADRG only supports up to %d entries in color table "
2762 : "(and an extra optional transparent one)",
2763 : CADRG_MAX_COLOR_ENTRY_COUNT);
2764 1 : return false;
2765 : }
2766 : }
2767 : }
2768 :
2769 126 : GDALColorTable oCT;
2770 63 : double dfLastPct = 0;
2771 :
2772 63 : constexpr int DEFAULT_QUANTIZATION_BITS = 5;
2773 : const char *pszBits =
2774 63 : CSLFetchNameValue(papszOptions, "COLOR_QUANTIZATION_BITS");
2775 63 : int nColorQuantizationBits = DEFAULT_QUANTIZATION_BITS;
2776 63 : if (pszBits)
2777 : {
2778 4 : nColorQuantizationBits = atoi(pszBits);
2779 4 : if (nColorQuantizationBits < 5 || nColorQuantizationBits > 8)
2780 : {
2781 0 : CPLError(CE_Failure, CPLE_AppDefined,
2782 : "COLOR_QUANTIZATION_BITS value must be between 5 and 8");
2783 0 : return false;
2784 : }
2785 : }
2786 :
2787 : class NITFDummyDataset final : public GDALDataset
2788 : {
2789 : public:
2790 22 : NITFDummyDataset() = default;
2791 : };
2792 :
2793 : const bool bInputIsShapedAndNameForCADRG =
2794 63 : poSrcDS->GetRasterXSize() == CADRG_FRAME_PIXEL_COUNT &&
2795 63 : poSrcDS->GetRasterYSize() == CADRG_FRAME_PIXEL_COUNT &&
2796 63 : bLooksLikeCADRGFilename;
2797 :
2798 : // Check if it is a whole blank frame, and do not generate it.
2799 63 : if (bInputIsShapedAndNameForCADRG)
2800 : {
2801 38 : bool bIsBlank = false;
2802 :
2803 38 : if (poSrcDS->GetRasterCount() == 4)
2804 : {
2805 1 : std::vector<GByte> abyData;
2806 1 : if (poSrcDS->GetRasterBand(4)->ReadRaster(abyData) != CE_None)
2807 0 : return false;
2808 :
2809 1 : bIsBlank = GDALBufferHasOnlyNoData(
2810 1 : abyData.data(), 0, poSrcDS->GetRasterXSize(),
2811 1 : poSrcDS->GetRasterYSize(),
2812 1 : /* stride = */ poSrcDS->GetRasterXSize(),
2813 : /* nComponents = */ 1,
2814 : /* nBitsPerSample = */ 8, GSF_UNSIGNED_INT);
2815 : }
2816 37 : else if (poSrcDS->GetRasterCount() == 3)
2817 : {
2818 2 : std::vector<GByte> abyData;
2819 : try
2820 : {
2821 2 : abyData.resize(poSrcDS->GetRasterCount() *
2822 2 : poSrcDS->GetRasterXSize() *
2823 2 : poSrcDS->GetRasterYSize());
2824 : }
2825 0 : catch (const std::exception &)
2826 : {
2827 0 : CPLError(CE_Failure, CPLE_OutOfMemory,
2828 : "Out of memory when allocating temporary buffer");
2829 0 : return false;
2830 : }
2831 4 : if (poSrcDS->RasterIO(
2832 : GF_Read, 0, 0, poSrcDS->GetRasterXSize(),
2833 2 : poSrcDS->GetRasterYSize(), abyData.data(),
2834 : poSrcDS->GetRasterXSize(), poSrcDS->GetRasterYSize(),
2835 : GDT_UInt8, poSrcDS->GetRasterCount(), nullptr,
2836 2 : poSrcDS->GetRasterCount(),
2837 2 : poSrcDS->GetRasterXSize() * poSrcDS->GetRasterCount(), 1,
2838 2 : nullptr) != CE_None)
2839 : {
2840 0 : return false;
2841 : }
2842 :
2843 : // Both (0,0,0) or (255,255,255) are considered blank
2844 4 : bIsBlank = GDALBufferHasOnlyNoData(
2845 2 : abyData.data(), 0, poSrcDS->GetRasterXSize(),
2846 2 : poSrcDS->GetRasterYSize(),
2847 2 : /* stride = */ poSrcDS->GetRasterXSize(),
2848 2 : /* nComponents = */ poSrcDS->GetRasterCount(),
2849 3 : /* nBitsPerSample = */ 8, GSF_UNSIGNED_INT) ||
2850 1 : GDALBufferHasOnlyNoData(
2851 1 : abyData.data(), 255, poSrcDS->GetRasterXSize(),
2852 1 : poSrcDS->GetRasterYSize(),
2853 1 : /* stride = */ poSrcDS->GetRasterXSize(),
2854 1 : /* nComponents = */ poSrcDS->GetRasterCount(),
2855 : /* nBitsPerSample = */ 8, GSF_UNSIGNED_INT);
2856 : }
2857 70 : else if (poCT &&
2858 35 : poCT->GetColorEntryCount() > CADRG_MAX_COLOR_ENTRY_COUNT)
2859 : {
2860 34 : std::vector<GByte> abyData;
2861 34 : if (poSrcDS->GetRasterBand(1)->ReadRaster(abyData) != CE_None)
2862 0 : return false;
2863 :
2864 34 : bIsBlank = GDALBufferHasOnlyNoData(
2865 34 : abyData.data(), TRANSPARENT_COLOR_TABLE_ENTRY,
2866 34 : poSrcDS->GetRasterXSize(), poSrcDS->GetRasterYSize(),
2867 34 : /* stride = */ poSrcDS->GetRasterXSize(),
2868 : /* nComponents = */ 1,
2869 : /* nBitsPerSample = */ 8, GSF_UNSIGNED_INT);
2870 : }
2871 :
2872 38 : if (bIsBlank)
2873 : {
2874 4 : CPLDebug("CADRG", "Skipping generation of %s as it is blank",
2875 : pszFilename);
2876 4 : return std::make_unique<NITFDummyDataset>();
2877 : }
2878 : }
2879 :
2880 59 : if (poSrcDS->GetRasterCount() == 1)
2881 : {
2882 36 : if (!poCT)
2883 : {
2884 1 : CPLError(CE_Failure, CPLE_AppDefined,
2885 : "CADRG only supports single band input datasets that "
2886 : "have an associated color table");
2887 1 : return false;
2888 : }
2889 35 : if (copyContext->oCT2.GetColorEntryCount() == 0)
2890 : {
2891 : // First time we go through this code path, we need to compute
2892 : // second and third color table from primary one.
2893 5 : dfLastPct = 0.1;
2894 : std::unique_ptr<void, decltype(&GDALDestroyScaledProgress)>
2895 : pScaledData(GDALCreateScaledProgress(0, dfLastPct, pfnProgress,
2896 : pProgressData),
2897 5 : GDALDestroyScaledProgress);
2898 5 : if (!ComputeColorTables(poSrcDS, oCT, nColorQuantizationBits,
2899 : GDALScaledProgress, pScaledData.get(),
2900 : *(copyContext)))
2901 : {
2902 0 : return false;
2903 : }
2904 : }
2905 : }
2906 23 : else if (poSrcDS->GetRasterCount() >= 3 &&
2907 23 : poSrcDS->GetRasterBand(1)->GetColorInterpretation() ==
2908 22 : GCI_RedBand &&
2909 22 : poSrcDS->GetRasterBand(2)->GetColorInterpretation() ==
2910 22 : GCI_GreenBand &&
2911 22 : poSrcDS->GetRasterBand(3)->GetColorInterpretation() ==
2912 46 : GCI_BlueBand &&
2913 22 : (poSrcDS->GetRasterCount() == 3 ||
2914 2 : (poSrcDS->GetRasterCount() == 4 &&
2915 2 : poSrcDS->GetRasterBand(4)->GetColorInterpretation() ==
2916 : GCI_AlphaBand)))
2917 : {
2918 22 : if (!bColorTablePerFrame || bInputIsShapedAndNameForCADRG)
2919 : {
2920 21 : dfLastPct = 0.1;
2921 : std::unique_ptr<void, decltype(&GDALDestroyScaledProgress)>
2922 : pScaledData(GDALCreateScaledProgress(0, dfLastPct, pfnProgress,
2923 : pProgressData),
2924 21 : GDALDestroyScaledProgress);
2925 21 : if (!ComputeColorTables(poSrcDS, oCT, nColorQuantizationBits,
2926 : GDALScaledProgress, pScaledData.get(),
2927 : *(copyContext)))
2928 : {
2929 0 : return false;
2930 : }
2931 : }
2932 : }
2933 : else
2934 : {
2935 1 : CPLError(CE_Failure, CPLE_AppDefined,
2936 : "CADRG only supports single-band paletted, RGB or RGBA "
2937 : "input datasets");
2938 1 : return false;
2939 : }
2940 :
2941 57 : if (bInputIsShapedAndNameForCADRG)
2942 : {
2943 34 : if (!poCT)
2944 : {
2945 : auto poPalettedDS =
2946 2 : CADRGGetPalettedDataset(poSrcDS, &oCT, nColorQuantizationBits);
2947 1 : if (!poPalettedDS)
2948 0 : return false;
2949 : std::unique_ptr<void, decltype(&GDALDestroyScaledProgress)>
2950 : pScaledData(GDALCreateScaledProgress(
2951 : dfLastPct, 1.0, pfnProgress, pProgressData),
2952 1 : GDALDestroyScaledProgress);
2953 2 : return NITFDataset::CreateCopy(
2954 : pszFilename, poPalettedDS.get(), bStrict, papszOptions,
2955 : GDALScaledProgress, pScaledData.get(), nRecLevel + 1,
2956 1 : copyContext);
2957 : }
2958 :
2959 33 : return true;
2960 : }
2961 : else
2962 : {
2963 23 : if (nReciprocalScale == 0)
2964 : {
2965 1 : CPLError(CE_Failure, CPLE_AppDefined, "SCALE must be defined");
2966 1 : return nullptr;
2967 : }
2968 :
2969 : const char *pszProducerCodeId =
2970 22 : CSLFetchNameValueDef(papszOptions, "PRODUCER_CODE_ID", "0");
2971 : const char *pszVersionNumber =
2972 22 : CSLFetchNameValueDef(papszOptions, "VERSION_NUMBER", "01");
2973 :
2974 22 : OGREnvelope sExtentNativeCRS;
2975 22 : OGREnvelope sExtentWGS84;
2976 44 : if (poSrcDS->GetExtent(&sExtentNativeCRS) != CE_None ||
2977 22 : poSrcDS->GetExtentWGS84LongLat(&sExtentWGS84) != CE_None)
2978 : {
2979 0 : CPLError(CE_Failure, CPLE_AppDefined, "Cannot get dataset extent");
2980 0 : return false;
2981 : }
2982 :
2983 22 : const char *pszZone = CSLFetchNameValue(papszOptions, "ZONE");
2984 : const int nZoneHintCandidate =
2985 43 : pszZone && strlen(pszZone) == 1 ? RPFCADRGZoneCharToNum(pszZone[0])
2986 21 : : pszZone && strlen(pszZone) == 2 ? atoi(pszZone)
2987 22 : : 0;
2988 : const int nZoneHint =
2989 22 : RPFCADRGIsValidZone(nZoneHintCandidate) ? nZoneHintCandidate : 0;
2990 : auto frameDefinitions =
2991 44 : RPFGetCADRGFramesForEnvelope(nZoneHint, nReciprocalScale, poSrcDS);
2992 :
2993 22 : if (nZoneHint)
2994 : {
2995 : // Only/mostly for debugging
2996 1 : const char *pszFrameX = CSLFetchNameValue(papszOptions, "FRAME_X");
2997 1 : const char *pszFrameY = CSLFetchNameValue(papszOptions, "FRAME_Y");
2998 1 : if (pszFrameX && pszFrameY)
2999 : {
3000 : RPFFrameDef frameDef;
3001 0 : frameDef.nReciprocalScale = nReciprocalScale;
3002 0 : frameDef.nZone = nZoneHint;
3003 0 : frameDef.nFrameMinX = atoi(pszFrameX);
3004 0 : frameDef.nFrameMinY = atoi(pszFrameY);
3005 0 : frameDef.nFrameMaxX = atoi(pszFrameX);
3006 0 : frameDef.nFrameMaxY = atoi(pszFrameY);
3007 0 : frameDef.dfResX = GetXPixelSize(nZoneHint, nReciprocalScale);
3008 0 : frameDef.dfResY = GetYPixelSize(nZoneHint, nReciprocalScale);
3009 :
3010 0 : frameDefinitions.clear();
3011 0 : frameDefinitions.push_back(std::move(frameDef));
3012 : }
3013 : }
3014 :
3015 22 : if (frameDefinitions.empty())
3016 : {
3017 2 : CPLError(CE_Failure, CPLE_AppDefined,
3018 : "Cannot establish CADRG frames intersecting dataset "
3019 : "extent");
3020 2 : return false;
3021 : }
3022 20 : else if (poSrcDS->GetRasterXSize() == CADRG_FRAME_PIXEL_COUNT &&
3023 4 : poSrcDS->GetRasterYSize() == CADRG_FRAME_PIXEL_COUNT &&
3024 2 : frameDefinitions.size() == 1 &&
3025 0 : frameDefinitions[0].nFrameMinX ==
3026 22 : frameDefinitions[0].nFrameMaxX &&
3027 0 : frameDefinitions[0].nFrameMinY ==
3028 0 : frameDefinitions[0].nFrameMaxY)
3029 : {
3030 : // poSrcDS already properly aligned
3031 0 : if (poCT && CPLGetExtensionSafe(pszFilename).size() == 3)
3032 0 : return true;
3033 :
3034 0 : std::string osFilename(pszFilename);
3035 0 : if (!bLooksLikeCADRGFilename)
3036 : {
3037 0 : const int nZone = frameDefinitions[0].nZone;
3038 0 : VSIMkdir(pszFilename, 0755);
3039 : const std::string osRPFDir =
3040 0 : CPLFormFilenameSafe(pszFilename, "RPF", nullptr);
3041 0 : VSIMkdir(osRPFDir.c_str(), 0755);
3042 : const std::string osZoneDir = CPLFormFilenameSafe(
3043 : osRPFDir.c_str(),
3044 0 : CPLSPrintf("ZONE%c", RPFCADRGZoneNumToChar(nZone)),
3045 0 : nullptr);
3046 0 : VSIMkdir(osZoneDir.c_str(), 0755);
3047 : VSIStatBufL sStat;
3048 0 : if (VSIStatL(osZoneDir.c_str(), &sStat) != 0 ||
3049 0 : !VSI_ISDIR(sStat.st_mode))
3050 : {
3051 0 : CPLError(CE_Failure, CPLE_AppDefined,
3052 : "Cannot create directory %s", osZoneDir.c_str());
3053 0 : return false;
3054 : }
3055 0 : osFilename = CPLFormFilenameSafe(
3056 : osZoneDir.c_str(),
3057 0 : RPFGetCADRGFrameNumberAsString(
3058 0 : nZone, nReciprocalScale, frameDefinitions[0].nFrameMinX,
3059 0 : frameDefinitions[0].nFrameMinY)
3060 : .c_str(),
3061 0 : nullptr);
3062 0 : osFilename += pszVersionNumber;
3063 0 : osFilename += pszProducerCodeId;
3064 0 : osFilename += '.';
3065 0 : osFilename += pszSeriesCode;
3066 0 : osFilename += RPFCADRGZoneNumToChar(nZone);
3067 0 : CPLDebug("CADRG", "Creating file %s", osFilename.c_str());
3068 : }
3069 :
3070 0 : if (bColorTablePerFrame)
3071 : {
3072 0 : return NITFDataset::CreateCopy(
3073 : osFilename.c_str(), poSrcDS, bStrict, papszOptions,
3074 0 : pfnProgress, pProgressData, nRecLevel + 1, copyContext);
3075 : }
3076 0 : else if (poCT)
3077 : {
3078 0 : return NITFDataset::CreateCopy(
3079 : osFilename.c_str(), poSrcDS, bStrict, papszOptions,
3080 0 : pfnProgress, pProgressData, nRecLevel + 1, copyContext);
3081 : }
3082 : else
3083 : {
3084 : auto poPalettedDS = CADRGGetPalettedDataset(
3085 0 : poSrcDS, &oCT, nColorQuantizationBits);
3086 0 : if (!poPalettedDS)
3087 0 : return false;
3088 0 : return NITFDataset::CreateCopy(
3089 : osFilename.c_str(), poPalettedDS.get(), bStrict,
3090 : papszOptions, pfnProgress, pProgressData, nRecLevel + 1,
3091 0 : copyContext);
3092 : }
3093 : }
3094 : else
3095 : {
3096 20 : int nTotalFrameCount = 0;
3097 43 : for (const auto &frameDef : frameDefinitions)
3098 : {
3099 23 : nTotalFrameCount +=
3100 23 : (frameDef.nFrameMaxX - frameDef.nFrameMinX + 1) *
3101 23 : (frameDef.nFrameMaxY - frameDef.nFrameMinY + 1);
3102 : }
3103 20 : CPLDebug("CADRG", "%d frame(s) to generate", nTotalFrameCount);
3104 :
3105 20 : int nCurFrameCounter = 0;
3106 :
3107 20 : VSIMkdir(pszFilename, 0755);
3108 : const std::string osRPFDir =
3109 40 : CPLFormFilenameSafe(pszFilename, "RPF", nullptr);
3110 20 : VSIMkdir(osRPFDir.c_str(), 0755);
3111 :
3112 20 : bool bMissingFramesFound = false;
3113 :
3114 41 : for (const auto &frameDef : frameDefinitions)
3115 : {
3116 23 : double dfXMin = 0, dfYMin = 0, dfXMax = 0, dfYMax = 0;
3117 : double dfTmpX, dfTmpY;
3118 23 : RPFGetCADRGFrameExtent(frameDef.nZone,
3119 23 : frameDef.nReciprocalScale,
3120 23 : frameDef.nFrameMinX, frameDef.nFrameMinY,
3121 : dfXMin, dfYMin, dfTmpX, dfTmpY);
3122 23 : RPFGetCADRGFrameExtent(frameDef.nZone,
3123 23 : frameDef.nReciprocalScale,
3124 23 : frameDef.nFrameMaxX, frameDef.nFrameMaxY,
3125 : dfTmpX, dfTmpY, dfXMax, dfYMax);
3126 :
3127 23 : CPLDebugOnly("CADRG", "Zone %d extent: %f,%f,%f,%f",
3128 : frameDef.nZone, dfXMin, dfYMin, dfXMax, dfYMax);
3129 :
3130 : auto poWarpedDS = CADRGGetWarpedVRTDataset(
3131 23 : poSrcDS, frameDef.nZone, dfXMin, dfYMin, dfXMax, dfYMax,
3132 23 : frameDef.dfResX, frameDef.dfResY, osResampling.c_str());
3133 23 : if (!poWarpedDS)
3134 0 : return false;
3135 :
3136 23 : if ((poWarpedDS->GetRasterXSize() % CADRG_FRAME_PIXEL_COUNT) !=
3137 46 : 0 ||
3138 23 : (poWarpedDS->GetRasterYSize() % CADRG_FRAME_PIXEL_COUNT) !=
3139 : 0)
3140 : {
3141 : // Should not happen unless there's a bug in our code
3142 0 : CPLError(CE_Failure, CPLE_AppDefined,
3143 : "Warped dataset dimension is %dx%d, which is "
3144 : "not a multiple of %dx%d",
3145 : poWarpedDS->GetRasterXSize(),
3146 : poWarpedDS->GetRasterYSize(),
3147 : CADRG_FRAME_PIXEL_COUNT, CADRG_FRAME_PIXEL_COUNT);
3148 0 : return false;
3149 : }
3150 :
3151 : const std::string osZoneDir = CPLFormFilenameSafe(
3152 : osRPFDir.c_str(),
3153 23 : CPLSPrintf("ZONE%c", RPFCADRGZoneNumToChar(frameDef.nZone)),
3154 23 : nullptr);
3155 : VSIStatBufL sStat;
3156 : const bool bDirectoryAlreadyExists =
3157 23 : (VSIStatL(osZoneDir.c_str(), &sStat) == 0);
3158 23 : if (!bDirectoryAlreadyExists)
3159 : {
3160 22 : VSIMkdir(osZoneDir.c_str(), 0755);
3161 43 : if (VSIStatL(osZoneDir.c_str(), &sStat) != 0 ||
3162 21 : !VSI_ISDIR(sStat.st_mode))
3163 : {
3164 1 : CPLError(CE_Failure, CPLE_AppDefined,
3165 : "Cannot create directory %s",
3166 : osZoneDir.c_str());
3167 1 : return false;
3168 : }
3169 : }
3170 :
3171 22 : CPLStringList aosOptions(papszOptions);
3172 : aosOptions.SetNameValue("ZONE",
3173 22 : CPLSPrintf("%d", frameDef.nZone));
3174 :
3175 22 : int nFrameCountThisZone =
3176 22 : (frameDef.nFrameMaxY - frameDef.nFrameMinY + 1) *
3177 22 : (frameDef.nFrameMaxX - frameDef.nFrameMinX + 1);
3178 :
3179 : #ifdef __COVERITY__
3180 : // Coverity has false positives about lambda captures by references
3181 : // being done without the lock being taken
3182 : CPL_IGNORE_RET_VAL(nFrameCountThisZone);
3183 : #else
3184 : // Limit to 4 as it doesn't scale beyond
3185 : const int nNumThreads =
3186 22 : std::min(GDALGetNumThreads(/* nMaxThreads = */ 4,
3187 : /* bDefaultToAllCPUs = */ true),
3188 22 : nFrameCountThisZone);
3189 :
3190 0 : CPLJobQueuePtr jobQueue;
3191 0 : std::unique_ptr<CPLWorkerThreadPool> poTP;
3192 22 : if (nNumThreads > 1)
3193 : {
3194 6 : poTP = std::make_unique<CPLWorkerThreadPool>();
3195 6 : if (poTP->Setup(nNumThreads, nullptr, nullptr))
3196 : {
3197 6 : jobQueue = poTP->CreateJobQueue();
3198 : }
3199 : }
3200 : #endif
3201 22 : std::mutex oMutex;
3202 22 : bool bError = false;
3203 22 : bool bErrorLocal = false;
3204 22 : int nCurFrameThisZone = 0;
3205 22 : std::string osErrorMsg;
3206 22 : nFrameCountThisZone = 0;
3207 22 : int nNonEmptyFrameCountThisZone = 0;
3208 :
3209 0 : std::unique_ptr<OGRCoordinateTransformation> poReprojCTToSrcSRS;
3210 22 : OGRSpatialReference oSRS;
3211 22 : if (frameDef.nZone == MAX_ZONE_NORTHERN_HEMISPHERE)
3212 1 : oSRS.importFromWkt(pszNorthPolarProjection);
3213 21 : else if (frameDef.nZone == MAX_ZONE)
3214 1 : oSRS.importFromWkt(pszSouthPolarProjection);
3215 : else
3216 20 : oSRS.importFromEPSG(4326);
3217 22 : oSRS.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
3218 22 : if (!poSrcSRS->IsSame(&oSRS) ||
3219 34 : frameDef.nZone == MAX_ZONE_NORTHERN_HEMISPHERE ||
3220 12 : frameDef.nZone == MAX_ZONE)
3221 : {
3222 10 : poReprojCTToSrcSRS.reset(
3223 : OGRCreateCoordinateTransformation(&oSRS, poSrcSRS));
3224 : }
3225 :
3226 48 : for (int iY = frameDef.nFrameMinY;
3227 48 : !bErrorLocal && iY <= frameDef.nFrameMaxY; ++iY)
3228 : {
3229 89 : for (int iX = frameDef.nFrameMinX;
3230 89 : !bErrorLocal && iX <= frameDef.nFrameMaxX; ++iX)
3231 : {
3232 : std::string osFilename = CPLFormFilenameSafe(
3233 : osZoneDir.c_str(),
3234 0 : RPFGetCADRGFrameNumberAsString(
3235 63 : frameDef.nZone, nReciprocalScale, iX, iY)
3236 : .c_str(),
3237 63 : nullptr);
3238 63 : osFilename += pszVersionNumber;
3239 63 : osFilename += pszProducerCodeId;
3240 63 : osFilename += '.';
3241 63 : osFilename += pszSeriesCode;
3242 63 : osFilename += RPFCADRGZoneNumToChar(frameDef.nZone);
3243 :
3244 63 : double dfFrameMinX = 0;
3245 63 : double dfFrameMinY = 0;
3246 63 : double dfFrameMaxX = 0;
3247 63 : double dfFrameMaxY = 0;
3248 63 : RPFGetCADRGFrameExtent(
3249 63 : frameDef.nZone, frameDef.nReciprocalScale, iX, iY,
3250 : dfFrameMinX, dfFrameMinY, dfFrameMaxX, dfFrameMaxY);
3251 :
3252 : bool bFrameFullyInSrcDS;
3253 63 : if (poReprojCTToSrcSRS)
3254 : {
3255 47 : OGREnvelope sFrameExtentInSrcSRS;
3256 47 : poReprojCTToSrcSRS->TransformBounds(
3257 : dfFrameMinX, dfFrameMinY, dfFrameMaxX,
3258 : dfFrameMaxY, &sFrameExtentInSrcSRS.MinX,
3259 : &sFrameExtentInSrcSRS.MinY,
3260 : &sFrameExtentInSrcSRS.MaxX,
3261 : &sFrameExtentInSrcSRS.MaxY,
3262 47 : DEFAULT_DENSIFY_PTS);
3263 :
3264 47 : if (!sFrameExtentInSrcSRS.Intersects(
3265 : sExtentNativeCRS))
3266 : {
3267 32 : CPLDebug(
3268 : "CADRG",
3269 : "Skipping creation of file %s as it does "
3270 : "not intersect source raster extent",
3271 : osFilename.c_str());
3272 32 : continue;
3273 : }
3274 15 : bFrameFullyInSrcDS =
3275 15 : sExtentNativeCRS.Contains(sFrameExtentInSrcSRS);
3276 : }
3277 : else
3278 : {
3279 16 : bFrameFullyInSrcDS =
3280 20 : dfFrameMinX >= sExtentWGS84.MinX &&
3281 4 : dfFrameMinY >= sExtentWGS84.MinY &&
3282 22 : dfFrameMaxX <= sExtentWGS84.MaxX &&
3283 2 : dfFrameMaxY <= sExtentWGS84.MaxY;
3284 : }
3285 :
3286 31 : ++nFrameCountThisZone;
3287 :
3288 : const auto task =
3289 31 : [osFilename = std::move(osFilename), copyContext,
3290 : dfFrameMinX, dfFrameMinY, dfFrameMaxX, dfFrameMaxY,
3291 : bFrameFullyInSrcDS, poCT, nRecLevel,
3292 : nColorQuantizationBits, bStrict,
3293 : bColorTablePerFrame, &oMutex, &poWarpedDS,
3294 : &nCurFrameThisZone, &nCurFrameCounter, &bError,
3295 : &oCT, &aosOptions, &bMissingFramesFound,
3296 424 : &osErrorMsg, &nNonEmptyFrameCountThisZone]()
3297 : {
3298 : #ifdef __COVERITY__
3299 : #define LOCK() \
3300 : do \
3301 : { \
3302 : (void)oMutex; \
3303 : } while (0)
3304 : #else
3305 : #define LOCK() std::lock_guard oLock(oMutex)
3306 : #endif
3307 : {
3308 31 : LOCK();
3309 31 : if (bError)
3310 0 : return;
3311 : }
3312 31 : CPLDebug("CADRG", "Creating file %s",
3313 : osFilename.c_str());
3314 :
3315 31 : CPLDebugOnly("CADRG", "Frame extent: %f %f %f %f",
3316 : dfFrameMinX, dfFrameMinY, dfFrameMaxX,
3317 : dfFrameMaxY);
3318 :
3319 0 : std::unique_ptr<GDALDataset> poClippedDS;
3320 : {
3321 : // Lock because poWarpedDS is not thread-safe
3322 31 : LOCK();
3323 62 : poClippedDS = CADRGGetClippedDataset(
3324 : poWarpedDS.get(), dfFrameMinX, dfFrameMinY,
3325 31 : dfFrameMaxX, dfFrameMaxY);
3326 : }
3327 31 : if (poClippedDS)
3328 : {
3329 31 : if (poCT)
3330 : {
3331 10 : if (!bFrameFullyInSrcDS)
3332 : {
3333 : // If the CADRG frame is not strictly inside
3334 : // the zone extent, add a transparent color
3335 9 : poClippedDS->GetRasterBand(1)
3336 9 : ->SetNoDataValue(
3337 9 : TRANSPARENT_COLOR_TABLE_ENTRY);
3338 9 : GDALColorEntry sEntry = {0, 0, 0, 0};
3339 9 : poClippedDS->GetRasterBand(1)
3340 9 : ->GetColorTable()
3341 9 : ->SetColorEntry(
3342 : TRANSPARENT_COLOR_TABLE_ENTRY,
3343 : &sEntry);
3344 : }
3345 : }
3346 21 : else if (!bColorTablePerFrame)
3347 : {
3348 40 : poClippedDS = CADRGGetPalettedDataset(
3349 : poClippedDS.get(), &oCT,
3350 20 : nColorQuantizationBits);
3351 : }
3352 : }
3353 :
3354 0 : std::unique_ptr<GDALDataset> poDS_CADRG;
3355 31 : if (poClippedDS)
3356 : {
3357 : CADRGCreateCopyContext copyContextCopy(
3358 31 : *copyContext);
3359 62 : poDS_CADRG = NITFDataset::CreateCopy(
3360 : osFilename.c_str(), poClippedDS.get(),
3361 31 : bStrict, aosOptions.List(), nullptr,
3362 31 : nullptr, nRecLevel + 1, copyContext);
3363 : }
3364 31 : if (!poDS_CADRG)
3365 : {
3366 1 : LOCK();
3367 1 : if (osErrorMsg.empty())
3368 1 : osErrorMsg = CPLGetLastErrorMsg();
3369 1 : bError = true;
3370 1 : return;
3371 : }
3372 : const bool bIsEmpty =
3373 30 : dynamic_cast<NITFDummyDataset *>(
3374 60 : poDS_CADRG.get()) != nullptr;
3375 30 : poDS_CADRG.reset();
3376 30 : VSIUnlink((osFilename + ".aux.xml").c_str());
3377 :
3378 60 : LOCK();
3379 30 : ++nCurFrameThisZone;
3380 30 : ++nCurFrameCounter;
3381 30 : if (bIsEmpty)
3382 : {
3383 1 : bMissingFramesFound = true;
3384 : }
3385 : else
3386 : {
3387 29 : ++nNonEmptyFrameCountThisZone;
3388 : }
3389 31 : };
3390 :
3391 : #ifndef __COVERITY__
3392 31 : if (jobQueue)
3393 : {
3394 15 : if (!jobQueue->SubmitJob(task))
3395 : {
3396 : {
3397 0 : std::lock_guard oLock(oMutex);
3398 0 : bError = true;
3399 : }
3400 0 : bErrorLocal = true;
3401 : }
3402 : }
3403 : else
3404 : #endif
3405 : {
3406 16 : task();
3407 16 : if (bError)
3408 1 : return false;
3409 30 : if (pfnProgress &&
3410 15 : !pfnProgress(
3411 : dfLastPct +
3412 30 : (1.0 - dfLastPct) * nCurFrameCounter /
3413 30 : std::max(1, nTotalFrameCount),
3414 : "", pProgressData))
3415 : {
3416 0 : return false;
3417 : }
3418 : }
3419 : }
3420 : }
3421 :
3422 : #ifndef __COVERITY__
3423 21 : if (jobQueue)
3424 : {
3425 : while (true)
3426 : {
3427 : {
3428 21 : std::lock_guard oLock(oMutex);
3429 42 : if (pfnProgress &&
3430 21 : !pfnProgress(
3431 : dfLastPct +
3432 42 : (1.0 - dfLastPct) * nCurFrameCounter /
3433 42 : std::max(1, nTotalFrameCount),
3434 : "", pProgressData))
3435 : {
3436 0 : bError = true;
3437 : }
3438 21 : if (bError ||
3439 21 : nCurFrameThisZone == nFrameCountThisZone)
3440 : break;
3441 : }
3442 15 : jobQueue->WaitEvent();
3443 15 : }
3444 :
3445 6 : jobQueue->WaitCompletion();
3446 6 : std::lock_guard oLock(oMutex);
3447 6 : if (bError)
3448 : {
3449 0 : if (!osErrorMsg.empty())
3450 : {
3451 0 : CPLError(CE_Failure, CPLE_AppDefined, "%s",
3452 : osErrorMsg.c_str());
3453 : }
3454 0 : return false;
3455 : }
3456 : }
3457 : #endif
3458 :
3459 21 : if (!bDirectoryAlreadyExists &&
3460 21 : nNonEmptyFrameCountThisZone == 0)
3461 : {
3462 1 : VSIRmdir(osZoneDir.c_str());
3463 : }
3464 : }
3465 :
3466 : const char *pszClassification =
3467 18 : CSLFetchNameValueDef(papszOptions, "FSCLAS", "U");
3468 18 : const char chIndexClassification =
3469 18 : pszClassification && pszClassification[0] ? pszClassification[0]
3470 : : 'U';
3471 36 : if (nCurFrameCounter > 0 &&
3472 54 : !RPFTOCCreate(
3473 : osRPFDir,
3474 36 : CPLFormFilenameSafe(osRPFDir.c_str(), "A.TOC", nullptr),
3475 : chIndexClassification, nReciprocalScale,
3476 : CSLFetchNameValueDef(papszOptions, "OSTAID",
3477 : ""), // producer id
3478 : CSLFetchNameValueDef(papszOptions, "ONAME",
3479 : ""), // producer name
3480 : CSLFetchNameValueDef(papszOptions, "SECURITY_COUNTRY_CODE",
3481 : ""),
3482 : /* bDoNotCreateIfNoFrame =*/bMissingFramesFound))
3483 : {
3484 0 : return nullptr;
3485 : }
3486 :
3487 18 : return std::make_unique<NITFDummyDataset>();
3488 : }
3489 : }
3490 : }
3491 : #endif
|
