1 /******************************************************************************
2  *
3  * Project:  JPEG-2000
4  * Purpose:  Implementation of the ISO/IEC 15444-1 standard based on Kakadu.
5  * Author:   Frank Warmerdam, warmerdam@pobox.com
6  *
7  ******************************************************************************
8  * Copyright (c) 2002, Frank Warmerdam <warmerdam@pobox.com>
9  * Copyright (c) 2007-2013, Even Rouault <even dot rouault at spatialys.com>
10  *
11  * Permission is hereby granted, free of charge, to any person obtaining a
12  * copy of this software and associated documentation files (the "Software"),
13  * to deal in the Software without restriction, including without limitation
14  * the rights to use, copy, modify, merge, publish, distribute, sublicense,
15  * and/or sell copies of the Software, and to permit persons to whom the
16  * Software is furnished to do so, subject to the following conditions:
17  *
18  * The above copyright notice and this permission notice shall be included
19  * in all copies or substantial portions of the Software.
20  *
21  * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
22  * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
23  * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
24  * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
25  * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
26  * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
27  * DEALINGS IN THE SOFTWARE.
28  ****************************************************************************/
29 
30 #ifdef DEBUG_BOOL
31 #define DO_NOT_USE_DEBUG_BOOL
32 #endif
33 
34 #include "cpl_port.h"
35 #include "jp2kakdataset.h"
36 
37 #include "cpl_multiproc.h"
38 #include "cpl_string.h"
39 #include "gdal_frmts.h"
40 #include "gdaljp2abstractdataset.h"
41 #include "gdaljp2metadata.h"
42 
43 #include "../mem/memdataset.h"
44 
45 #include "jp2kak_headers.h"
46 
47 #include "subfile_source.h"
48 #include "vsil_target.h"
49 
50 #include <cstdlib>
51 #include <cstring>
52 #include <algorithm>
53 #include <cmath>
54 #include <vector>
55 
56 CPL_CVSID("$Id: jp2kakdataset.cpp 126b0897e64c233ed06ca072549e110bb6b28ced 2021-04-20 16:42:23 +0200 Even Rouault $")
57 
58 // Before v7.5 Kakadu does not advertise its version well
59 // After v7.5 Kakadu has KDU_{MAJOR,MINOR,PATCH}_VERSION defines so it is easier
60 // For older releases compile with them manually specified
61 #ifndef KDU_MAJOR_VERSION
62 #  error Compile with eg. -DKDU_MAJOR_VERSION=7 -DKDU_MINOR_VERSION=3 -DKDU_PATCH_VERSION=2 to specify Kakadu library version
63 #endif
64 
65 #if KDU_MAJOR_VERSION > 7 || (KDU_MAJOR_VERSION == 7 && KDU_MINOR_VERSION >= 5)
66     using namespace kdu_core;
67     using namespace kdu_supp;
68 #endif
69 
70 // #define KAKADU_JPX 1
71 
72 static bool kakadu_initialized = false;
73 
74 constexpr unsigned char jp2_header[] =
75     {0x00, 0x00, 0x00, 0x0c, 0x6a, 0x50, 0x20, 0x20, 0x0d, 0x0a, 0x87, 0x0a};
76 
77 constexpr unsigned char jpc_header[] = {0xff, 0x4f};
78 
79 /* -------------------------------------------------------------------- */
80 /*      The number of tiles at a time we will push through the          */
81 /*      encoder per flush when writing jpeg2000 streams.                */
82 /* -------------------------------------------------------------------- */
83 constexpr int TILE_CHUNK_SIZE = 1024;
84 
85 /************************************************************************/
86 /* ==================================================================== */
87 /*                            JP2KAKRasterBand                          */
88 /* ==================================================================== */
89 /************************************************************************/
90 
91 /************************************************************************/
92 /*                           JP2KAKRasterBand()                         */
93 /************************************************************************/
94 
95 JP2KAKRasterBand::JP2KAKRasterBand( int nBandIn, int nDiscardLevelsIn,
96                                     kdu_codestream oCodeStreamIn,
97                                     int nResCount, kdu_client *jpip_clientIn,
98                                     jp2_channels oJP2Channels,
99                                     JP2KAKDataset *poBaseDSIn ) :
100     poBaseDS(poBaseDSIn),
101     nDiscardLevels(nDiscardLevelsIn),
102     nOverviewCount(0),
103     papoOverviewBand(nullptr),
104     jpip_client(jpip_clientIn),
105     oCodeStream(oCodeStreamIn),
106     eInterp(GCI_Undefined)
107 {
108     nBand = nBandIn;  // From GDALRasterBand.
109 
110     if( oCodeStream.get_bit_depth(nBand - 1) > 8 &&
111         oCodeStream.get_bit_depth(nBand - 1) <= 16 &&
112         oCodeStream.get_signed(nBand - 1) )
113         eDataType = GDT_Int16;
114     else if( oCodeStream.get_bit_depth(nBand - 1) > 8 &&
115              oCodeStream.get_bit_depth(nBand - 1) <= 16 &&
116              !oCodeStream.get_signed(nBand - 1) )
117         eDataType = GDT_UInt16;
118     else
119         eDataType = GDT_Byte;
120 
121     oCodeStream.apply_input_restrictions(0, 0, nDiscardLevels, 0, nullptr);
122     oCodeStream.get_dims(0, band_dims);
123 
124     nRasterXSize = band_dims.size.x;
125     nRasterYSize = band_dims.size.y;
126 
127     // Capture some useful metadata.
128     if( oCodeStream.get_bit_depth(nBand - 1) % 8 != 0 &&
129         !poBaseDSIn->bPromoteTo8Bit )
130     {
131         SetMetadataItem(
132             "NBITS",
133             CPLString().Printf("%d", oCodeStream.get_bit_depth(nBand - 1)),
134             "IMAGE_STRUCTURE");
135     }
136     SetMetadataItem("COMPRESSION", "JP2000", "IMAGE_STRUCTURE");
137 
138     // Use tile dimension as block size, unless it is too big
139     kdu_dims valid_tiles;
140     kdu_dims tile_dims;
141     oCodeStream.get_valid_tiles(valid_tiles);
142     oCodeStream.get_tile_dims(valid_tiles.pos, -1, tile_dims);
143     // Configuration option only for testing purposes
144     if( CPLTestBool(CPLGetConfigOption("USE_TILE_AS_BLOCK", "NO")) )
145     {
146         nBlockXSize = std::min(tile_dims.size.x, nRasterXSize);
147         nBlockYSize = std::min(tile_dims.size.y, nRasterYSize);
148     }
149     else
150     {
151         nBlockXSize = std::min(std::min(tile_dims.size.x, 2048), nRasterXSize);
152         nBlockYSize = std::min(std::min(tile_dims.size.y, 2048), nRasterYSize);
153     }
154     CPLDebug( "JP2KAK", "JP2KAKRasterBand::JP2KAKRasterBand() : "
155             "Tile dimension : %d X %d\n",
156             nBlockXSize, nBlockYSize);
157 
158     // Figure out the color interpretation for this band.
159     eInterp = GCI_Undefined;
160 
161     if( oJP2Channels.exists() )
162     {
163         int nRedIndex = -1;
164         int nGreenIndex = -1;
165         int nBlueIndex = -1;
166         int nLutIndex = 0;
167         int nCSI = 0;
168 
169 #if KDU_MAJOR_VERSION > 7 || (KDU_MAJOR_VERSION == 7 && KDU_MINOR_VERSION >= 8)
170         int nFMT = 0;
171         if( oJP2Channels.get_num_colours() == 3 )
172         {
173             oJP2Channels.get_colour_mapping(0, nRedIndex, nLutIndex, nCSI,
174                                             nFMT);
175             oJP2Channels.get_colour_mapping(1, nGreenIndex, nLutIndex, nCSI,
176                                             nFMT);
177             oJP2Channels.get_colour_mapping(2, nBlueIndex, nLutIndex, nCSI,
178                                             nFMT);
179         }
180         else
181         {
182             oJP2Channels.get_colour_mapping(0, nRedIndex, nLutIndex, nCSI,
183                                             nFMT);
184             if( nBand == 1 )
185                 eInterp = GCI_GrayIndex;
186         }
187 #else
188         if( oJP2Channels.get_num_colours() == 3 )
189         {
190             oJP2Channels.get_colour_mapping(0, nRedIndex, nLutIndex, nCSI);
191             oJP2Channels.get_colour_mapping(1, nGreenIndex, nLutIndex, nCSI);
192             oJP2Channels.get_colour_mapping(2, nBlueIndex, nLutIndex, nCSI);
193         }
194         else
195         {
196             oJP2Channels.get_colour_mapping(0, nRedIndex, nLutIndex, nCSI);
197             if( nBand == 1 )
198                 eInterp = GCI_GrayIndex;
199         }
200 #endif
201 
202         if( eInterp != GCI_Undefined )
203             /* nothing to do */;
204         // If we have LUT info, it is a palette image.
205         else if( nLutIndex != -1 )
206             eInterp = GCI_PaletteIndex;
207         // Establish color band this is.
208         else if( nRedIndex == nBand - 1 )
209             eInterp = GCI_RedBand;
210         else if( nGreenIndex == nBand - 1 )
211             eInterp = GCI_GreenBand;
212         else if( nBlueIndex == nBand - 1 )
213             eInterp = GCI_BlueBand;
214         else
215             eInterp = GCI_Undefined;
216 
217         // Could this band be an alpha band?
218         if( eInterp == GCI_Undefined )
219         {
220             for( int color_idx = 0;
221                  color_idx < oJP2Channels.get_num_colours(); color_idx++ )
222             {
223                 int opacity_idx = 0;
224                 int lut_idx = 0;
225 
226                 // get_opacity_mapping sets that last 3 args by non-const refs.
227 #if KDU_MAJOR_VERSION > 7 || (KDU_MAJOR_VERSION == 7 && KDU_MINOR_VERSION >= 8)
228                 if(oJP2Channels.get_opacity_mapping(color_idx, opacity_idx,
229                                                     lut_idx, nCSI, nFMT))
230 #else
231                 if(oJP2Channels.get_opacity_mapping(color_idx, opacity_idx,
232                                                     lut_idx, nCSI))
233 #endif
234                 {
235                     if( opacity_idx == nBand - 1 )
236                         eInterp = GCI_AlphaBand;
237                 }
238 #if KDU_MAJOR_VERSION > 7 || (KDU_MAJOR_VERSION == 7 && KDU_MINOR_VERSION >= 8)
239                 if( oJP2Channels.get_premult_mapping(color_idx, opacity_idx,
240                                                      lut_idx, nCSI, nFMT) )
241 #else
242                 if( oJP2Channels.get_premult_mapping(color_idx, opacity_idx,
243                                                      lut_idx, nCSI) )
244 #endif
245                 {
246                     if( opacity_idx == nBand - 1 )
247                         eInterp = GCI_AlphaBand;
248                 }
249             }
250         }
251     }
252     else if( nBand == 1 )
253     {
254         eInterp = GCI_RedBand;
255     }
256     else if( nBand == 2 )
257     {
258         eInterp = GCI_GreenBand;
259     }
260     else if( nBand == 3 )
261     {
262         eInterp = GCI_BlueBand;
263     }
264     else
265     {
266         eInterp = GCI_GrayIndex;
267     }
268 
269     // Do we have any overviews?  Only check if we are the full res image.
270     if( nDiscardLevels == 0 && GDALPamRasterBand::GetOverviewCount() == 0 )
271     {
272         int nXSize = nRasterXSize;
273         int nYSize = nRasterYSize;
274 
275         for( int nDiscard = 1; nDiscard < nResCount; nDiscard++ )
276         {
277             nXSize = (nXSize + 1) / 2;
278             nYSize = (nYSize + 1) / 2;
279 
280             if( (nXSize + nYSize) < 128 || nXSize < 4 || nYSize < 4 )
281                 continue;  // Skip super reduced resolution layers.
282 
283             oCodeStream.apply_input_restrictions(0, 0, nDiscard, 0, nullptr);
284             kdu_dims dims;  // Struct with default constructor.
285             oCodeStream.get_dims(0, dims);
286 
287             if( (dims.size.x == nXSize || dims.size.x == nXSize - 1) &&
288                 (dims.size.y == nYSize || dims.size.y == nYSize - 1) )
289             {
290                 nOverviewCount++;
291                 papoOverviewBand = static_cast<JP2KAKRasterBand **>(CPLRealloc(
292                     papoOverviewBand, sizeof(void *) * nOverviewCount));
293                 papoOverviewBand[nOverviewCount - 1] =
294                     new JP2KAKRasterBand(nBand, nDiscard, oCodeStream, 0,
295                                          jpip_client, oJP2Channels, poBaseDS);
296             }
297             else
298             {
299                 CPLDebug("GDAL", "Discard %dx%d JPEG2000 overview layer,\n"
300                          "expected %dx%d.",
301                          dims.size.x, dims.size.y, nXSize, nYSize);
302             }
303         }
304     }
305 }
306 
307 /************************************************************************/
308 /*                         ~JP2KAKRasterBand()                          */
309 /************************************************************************/
310 
311 JP2KAKRasterBand::~JP2KAKRasterBand()
312 
313 {
314     for( int i = 0; i < nOverviewCount; i++ )
315         delete papoOverviewBand[i];
316 
317     CPLFree(papoOverviewBand);
318 }
319 
320 /************************************************************************/
321 /*                          GetOverviewCount()                          */
322 /************************************************************************/
323 
324 int JP2KAKRasterBand::GetOverviewCount()
325 
326 {
327     if( !poBaseDS->AreOverviewsEnabled() )
328         return 0;
329 
330     if( GDALPamRasterBand::GetOverviewCount() > 0 )
331         return GDALPamRasterBand::GetOverviewCount();
332 
333     return nOverviewCount;
334 }
335 
336 /************************************************************************/
337 /*                            GetOverview()                             */
338 /************************************************************************/
339 
340 GDALRasterBand *JP2KAKRasterBand::GetOverview( int iOverviewIndex )
341 
342 {
343     if( GDALPamRasterBand::GetOverviewCount() > 0 )
344         return GDALPamRasterBand::GetOverview(iOverviewIndex);
345 
346     if( iOverviewIndex < 0 || iOverviewIndex >= nOverviewCount )
347         return nullptr;
348 
349     return papoOverviewBand[iOverviewIndex];
350 }
351 
352 /************************************************************************/
353 /*                             IReadBlock()                             */
354 /************************************************************************/
355 
356 CPLErr JP2KAKRasterBand::IReadBlock( int nBlockXOff, int nBlockYOff,
357                                      void * pImage )
358 {
359     const int nWordSize = GDALGetDataTypeSizeBytes(eDataType);
360     int nOvMult = 1;
361     int nLevelsLeft = nDiscardLevels;
362     while( nLevelsLeft-- > 0 )
363         nOvMult *= 2;
364 
365     CPLDebug("JP2KAK", "IReadBlock(%d,%d) on band %d.",
366              nBlockXOff, nBlockYOff, nBand);
367 
368     // Compute the normal window, and buffer size.
369     const int nWXOff = nBlockXOff * nBlockXSize * nOvMult;
370     const int nWYOff = nBlockYOff * nBlockYSize * nOvMult;
371     int nWXSize = nBlockXSize * nOvMult;
372     int nWYSize = nBlockYSize * nOvMult;
373 
374     int nXSize = nBlockXSize;
375     int nYSize = nBlockYSize;
376 
377     // Adjust if we have a partial block on the right or bottom of
378     // the image.  Unfortunately despite some care I can't seem to
379     // always get partial tiles to come from the desired overview
380     // level depending on how various things round - hopefully not
381     // a big deal.
382     if( nWXOff + nWXSize > poBaseDS->GetRasterXSize() )
383     {
384         nWXSize = poBaseDS->GetRasterXSize() - nWXOff;
385         nXSize = nRasterXSize - nBlockXSize * nBlockXOff;
386     }
387 
388     if( nWYOff + nWYSize > poBaseDS->GetRasterYSize() )
389     {
390         nWYSize = poBaseDS->GetRasterYSize() - nWYOff;
391         nYSize = nRasterYSize - nBlockYSize * nBlockYOff;
392     }
393 
394     if( nXSize != nBlockXSize || nYSize != nBlockYSize )
395         memset(pImage, 0, static_cast<size_t>(nBlockXSize) * nBlockYSize * nWordSize);
396 
397     // By default we invoke just for the requested band, directly
398     // into the target buffer.
399     GDALRasterIOExtraArg sExtraArg;
400     INIT_RASTERIO_EXTRA_ARG(sExtraArg);
401 
402     if( !poBaseDS->bUseYCC )
403     {
404         return poBaseDS->DirectRasterIO(GF_Read,
405                                         nWXOff, nWYOff, nWXSize, nWYSize,
406                                         pImage, nXSize, nYSize,
407                                         eDataType, 1, &nBand,
408                                         nWordSize, nWordSize*nBlockXSize,
409                                         0, &sExtraArg);
410     }
411 
412     // But for YCC or possible other effectively pixel interleaved
413     // products, we read all bands into a single buffer, fetch out
414     // what we want, and push the rest into the block cache.
415     std::vector<int> anBands;
416 
417     for( int iBand = 0; iBand < poBaseDS->GetRasterCount(); iBand++ )
418     {
419         GDALRasterBand *poBand = poBaseDS->GetRasterBand(iBand + 1);
420         if ( poBand->GetRasterDataType() != eDataType )
421           continue;
422         anBands.push_back(iBand + 1);
423     }
424 
425     GByte *pabyWrkBuffer = static_cast<GByte *>(
426         VSIMalloc3(nWordSize * anBands.size(), nBlockXSize, nBlockYSize));
427     if( pabyWrkBuffer == nullptr )
428         return CE_Failure;
429 
430     const CPLErr eErr = poBaseDS->DirectRasterIO(
431         GF_Read, nWXOff, nWYOff, nWXSize, nWYSize, pabyWrkBuffer, nXSize,
432         nYSize, eDataType, static_cast<int>(anBands.size()), &anBands[0],
433         nWordSize, nWordSize * nBlockXSize,
434         static_cast<GSpacing>(nWordSize) * nBlockXSize * nBlockYSize, &sExtraArg);
435 
436     if( eErr == CE_None )
437     {
438         int nBandStart = 0;
439         const int nTotalBands = static_cast<int>(anBands.size());
440         for( int iBand = 0; iBand < nTotalBands; iBand++ )
441         {
442             if( anBands[iBand] == nBand )
443             {
444                 // Application requested band.
445                 memcpy(pImage, pabyWrkBuffer + nBandStart,
446                        static_cast<size_t>(nWordSize) * nBlockXSize * nBlockYSize);
447             }
448             else
449             {
450                 // All others are pushed into cache.
451                 GDALRasterBand *poBaseBand =
452                     poBaseDS->GetRasterBand(anBands[iBand]);
453                 JP2KAKRasterBand *poBand = nullptr;
454 
455                 if( nDiscardLevels == 0 )
456                 {
457                     poBand = cpl::down_cast<JP2KAKRasterBand *>(poBaseBand);
458                 }
459                 else
460                 {
461                     int iOver = 0;  // Used after for.
462 
463                     for( ; iOver < poBaseBand->GetOverviewCount(); iOver++ )
464                     {
465                         poBand = cpl::down_cast<JP2KAKRasterBand *>(
466                             poBaseBand->GetOverview( iOver ) );
467                         if( poBand->nDiscardLevels == nDiscardLevels )
468                             break;
469                     }
470                     if( iOver == poBaseBand->GetOverviewCount() )
471                     {
472                         CPLAssert(false);
473                     }
474                 }
475 
476                 GDALRasterBlock *poBlock = nullptr;
477 
478                 if( poBand != nullptr )
479                     poBlock =
480                         poBand->GetLockedBlockRef(nBlockXOff, nBlockYOff, TRUE);
481 
482                 if( poBlock )
483                 {
484                     memcpy(poBlock->GetDataRef(), pabyWrkBuffer + nBandStart,
485                            static_cast<size_t>(nWordSize) * nBlockXSize * nBlockYSize);
486                     poBlock->DropLock();
487                 }
488             }
489 
490             nBandStart += static_cast<size_t>(nWordSize) * nBlockXSize * nBlockYSize;
491         }
492     }
493 
494     VSIFree(pabyWrkBuffer);
495     return eErr;
496 }
497 
498 /************************************************************************/
499 /*                             IRasterIO()                              */
500 /************************************************************************/
501 
502 CPLErr
503 JP2KAKRasterBand::IRasterIO( GDALRWFlag eRWFlag,
504                              int nXOff, int nYOff, int nXSize, int nYSize,
505                              void * pData, int nBufXSize, int nBufYSize,
506                              GDALDataType eBufType,
507                              GSpacing nPixelSpace, GSpacing nLineSpace,
508                              GDALRasterIOExtraArg* psExtraArg)
509 
510 {
511     // We need various criteria to skip out to block based methods.
512     if( poBaseDS->TestUseBlockIO( nXOff, nYOff, nXSize, nYSize,
513                                   nBufXSize, nBufYSize,
514                                   eBufType, 1, &nBand ) )
515         return GDALPamRasterBand::IRasterIO(
516             eRWFlag, nXOff, nYOff, nXSize, nYSize,
517             pData, nBufXSize, nBufYSize, eBufType,
518             nPixelSpace, nLineSpace, psExtraArg );
519 
520     int nOverviewDiscard = nDiscardLevels;
521 
522     // Adjust request for overview level.
523     while( nOverviewDiscard > 0 )
524     {
525         nXOff *= 2;
526         nYOff *= 2;
527         nXSize *= 2;
528         nYSize *= 2;
529         nOverviewDiscard--;
530     }
531 
532     return poBaseDS->DirectRasterIO(
533         eRWFlag, nXOff, nYOff, nXSize, nYSize,
534         pData, nBufXSize, nBufYSize, eBufType,
535         1, &nBand, nPixelSpace, nLineSpace, 0, psExtraArg );
536 }
537 
538 /************************************************************************/
539 /*                            ApplyPalette()                            */
540 /************************************************************************/
541 
542 namespace {
543 
544 inline short GetColorValue(const float *pafLUT, int nPos)
545 {
546     const short nVal = static_cast<short>(pafLUT[nPos] * 256.0f + 128.0f);
547     const short nMin = 0;
548     const short nMax = 255;
549     return std::max(nMin, std::min(nMax, nVal));
550 }
551 
552 }  // namespace
553 
554 void JP2KAKRasterBand::ApplyPalette( jp2_palette oJP2Palette )
555 
556 {
557     // Do we have a reasonable LUT configuration?  RGB or RGBA?
558     if( !oJP2Palette.exists() )
559         return;
560 
561     if( oJP2Palette.get_num_luts() == 0 || oJP2Palette.get_num_entries() == 0 )
562         return;
563 
564     if( oJP2Palette.get_num_luts() < 3 )
565     {
566         CPLDebug("JP2KAK", "JP2KAKRasterBand::ApplyPalette()\n"
567                  "Odd get_num_luts() value (%d)",
568                  oJP2Palette.get_num_luts());
569         return;
570     }
571 
572     // Fetch lut entries.  They are normalized in the -0.5 to 0.5 range.                                              */
573     const int nCount = oJP2Palette.get_num_entries();
574 
575     float *const pafLUT =
576         static_cast<float *>(CPLCalloc(sizeof(float) * 4, nCount));
577 
578     const int nRed = 0;
579     const int nGreen = 1;
580     const int nBlue = 2;
581     const int nAlpha = 3;
582     oJP2Palette.get_lut(nRed, pafLUT + 0);
583     oJP2Palette.get_lut(nGreen, pafLUT + nCount);
584     oJP2Palette.get_lut(nBlue, pafLUT + nCount * 2);
585 
586     if( oJP2Palette.get_num_luts() == 4 )
587     {
588         oJP2Palette.get_lut(nAlpha, pafLUT + nCount * 3);
589     }
590     else
591     {
592         for( int iColor = 0; iColor < nCount; iColor++ )
593         {
594             pafLUT[nCount * 3 + iColor] = 0.5;
595         }
596     }
597 
598     // Apply to GDAL colortable.
599     const int nRedOffset = nCount * nRed;
600     const int nGreenOffset = nCount * nGreen;
601     const int nBlueOffset = nCount * nBlue;
602     const int nAlphaOffset = nCount * nAlpha;
603 
604     for( int iColor = 0; iColor < nCount; iColor++ )
605     {
606         const GDALColorEntry sEntry = {
607             GetColorValue(pafLUT, iColor + nRedOffset),
608             GetColorValue(pafLUT, iColor + nGreenOffset),
609             GetColorValue(pafLUT, iColor + nBlueOffset),
610             GetColorValue(pafLUT, iColor + nAlphaOffset)
611         };
612 
613         oCT.SetColorEntry(iColor, &sEntry);
614     }
615 
616     CPLFree(pafLUT);
617 
618     eInterp = GCI_PaletteIndex;
619 }
620 
621 /************************************************************************/
622 /*                       GetColorInterpretation()                       */
623 /************************************************************************/
624 
625 GDALColorInterp JP2KAKRasterBand::GetColorInterpretation() { return eInterp; }
626 
627 /************************************************************************/
628 /*                           GetColorTable()                            */
629 /************************************************************************/
630 
631 GDALColorTable *JP2KAKRasterBand::GetColorTable()
632 
633 {
634     if( oCT.GetColorEntryCount() > 0 )
635         return &oCT;
636 
637     return nullptr;
638 }
639 
640 /************************************************************************/
641 /* ==================================================================== */
642 /*                           JP2KAKDataset                              */
643 /* ==================================================================== */
644 /************************************************************************/
645 
646 /************************************************************************/
647 /*                           JP2KAKDataset()                           */
648 /************************************************************************/
649 
650 JP2KAKDataset::JP2KAKDataset()
651 
652 {
653     poDriver = static_cast<GDALDriver *>(GDALGetDriverByName("JP2KAK"));
654 }
655 
656 /************************************************************************/
657 /*                            ~JP2KAKDataset()                         */
658 /************************************************************************/
659 
660 JP2KAKDataset::~JP2KAKDataset()
661 
662 {
663     FlushCache();
664 
665     if( poInput != nullptr )
666     {
667         oCodeStream.destroy();
668         poInput->close();
669         delete poInput;
670         if( family )
671         {
672             family->close();
673             delete family;
674         }
675         if( poRawInput != nullptr )
676             delete poRawInput;
677 #ifdef USE_JPIP
678         if( jpip_client != NULL )
679         {
680             jpip_client->close();
681             delete jpip_client;
682         }
683 #endif
684     }
685 
686     if( poThreadEnv != nullptr )
687     {
688         poThreadEnv->terminate(nullptr, true);
689         poThreadEnv->destroy();
690         delete poThreadEnv;
691     }
692 }
693 
694 /************************************************************************/
695 /*                          IBuildOverviews()                           */
696 /************************************************************************/
697 
698 CPLErr JP2KAKDataset::IBuildOverviews( const char *pszResampling,
699                                        int nOverviews, int *panOverviewList,
700                                        int nListBands, int *panBandList,
701                                        GDALProgressFunc pfnProgress,
702                                        void *pProgressData )
703 
704 {
705     // In order for building external overviews to work properly, we
706     // discard any concept of internal overviews when the user
707     // first requests to build external overviews.
708     for( int iBand = 0; iBand < GetRasterCount(); iBand++ )
709     {
710         JP2KAKRasterBand *poBand =
711             cpl::down_cast<JP2KAKRasterBand *>(GetRasterBand(iBand + 1));
712         for( int i = 0; i < poBand->nOverviewCount; i++ )
713             delete poBand->papoOverviewBand[i];
714 
715         CPLFree(poBand->papoOverviewBand);
716         poBand->papoOverviewBand = nullptr;
717         poBand->nOverviewCount = 0;
718     }
719 
720     return GDALPamDataset::IBuildOverviews(pszResampling,
721                                            nOverviews, panOverviewList,
722                                            nListBands, panBandList,
723                                            pfnProgress, pProgressData);
724 }
725 
726 /************************************************************************/
727 /*                          KakaduInitialize()                          */
728 /************************************************************************/
729 
730 void JP2KAKDataset::KakaduInitialize()
731 
732 {
733     // Initialize Kakadu warning/error reporting subsystem.
734     if( kakadu_initialized )
735         return;
736 
737     kakadu_initialized = true;
738 
739     kdu_cpl_error_message oErrHandler(CE_Failure);
740     kdu_cpl_error_message oWarningHandler(CE_Warning);
741     CPL_IGNORE_RET_VAL(oErrHandler);
742     CPL_IGNORE_RET_VAL(oWarningHandler);
743 
744     kdu_customize_warnings(new kdu_cpl_error_message(CE_Warning));
745     kdu_customize_errors(new kdu_cpl_error_message(CE_Failure));
746 }
747 
748 /************************************************************************/
749 /*                              Identify()                              */
750 /************************************************************************/
751 
752 int JP2KAKDataset::Identify( GDALOpenInfo * poOpenInfo )
753 
754 {
755     // Check header.
756     if( poOpenInfo->nHeaderBytes < static_cast<int>(sizeof(jp2_header)) )
757     {
758         if( (STARTS_WITH_CI(poOpenInfo->pszFilename, "http://")
759              || STARTS_WITH_CI(poOpenInfo->pszFilename, "https://")
760              || STARTS_WITH_CI(poOpenInfo->pszFilename, "jpip://"))
761             && EQUAL(CPLGetExtension( poOpenInfo->pszFilename ), "jp2") )
762         {
763 #ifdef USE_JPIP
764             return TRUE;
765 #else
766             return FALSE;
767 #endif
768         }
769         else if( STARTS_WITH_CI(poOpenInfo->pszFilename, "J2K_SUBFILE:") )
770             return TRUE;
771         else
772             return FALSE;
773     }
774 
775 /* -------------------------------------------------------------------- */
776 /*      Any extension is supported for JP2 files.  Only selected        */
777 /*      extensions are supported for JPC files since the standard       */
778 /*      prefix is so short (two bytes).                                 */
779 /* -------------------------------------------------------------------- */
780     if( memcmp(poOpenInfo->pabyHeader, jp2_header, sizeof(jp2_header)) == 0 )
781         return TRUE;
782     else if( memcmp(poOpenInfo->pabyHeader, jpc_header,
783                     sizeof(jpc_header)) == 0 )
784     {
785         const char *const pszExtension =
786             CPLGetExtension(poOpenInfo->pszFilename);
787         if( EQUAL(pszExtension,"jpc")
788             || EQUAL(pszExtension, "j2k")
789             || EQUAL(pszExtension, "jp2")
790             || EQUAL(pszExtension, "jpx")
791             || EQUAL(pszExtension, "j2c") )
792            return TRUE;
793 
794         // We also want to handle jpc datastreams vis /vsisubfile.
795         if( strstr(poOpenInfo->pszFilename, "vsisubfile") != nullptr )
796             return TRUE;
797     }
798 
799     return FALSE;
800 }
801 
802 /************************************************************************/
803 /*                                Open()                                */
804 /************************************************************************/
805 
806 GDALDataset *JP2KAKDataset::Open( GDALOpenInfo * poOpenInfo )
807 
808 {
809 #ifndef FUZZING_BUILD_MODE_UNSAFE_FOR_PRODUCTION
810     // During fuzzing, do not use Identify to reject crazy content.
811     if( !Identify(poOpenInfo) )
812         return nullptr;
813 #endif
814 
815     subfile_source *poRawInput = nullptr;
816     bool bIsJPIP = false;
817     bool bIsSubfile = false;
818     const GByte *pabyHeader = nullptr;
819 
820     const bool bResilient =
821         CPLTestBool(CPLGetConfigOption("JP2KAK_RESILIENT", "NO"));
822 
823     // Doesn't seem to bring any real performance gain on Linux.
824     const bool bBuffered = CPLTestBool(CPLGetConfigOption("JP2KAK_BUFFERED",
825 #ifdef WIN32
826                                                           "YES"
827 #else
828                                                           "NO"
829 #endif
830                                                           ));
831 
832     KakaduInitialize();
833 
834     // Handle setting up datasource for JPIP.
835     const char *pszExtension = CPLGetExtension(poOpenInfo->pszFilename);
836     std::vector<GByte> abySubfileHeader(16); // leave in this scope
837     if( poOpenInfo->nHeaderBytes < 16 )
838     {
839         if( (STARTS_WITH_CI(poOpenInfo->pszFilename, "http://")
840              || STARTS_WITH_CI(poOpenInfo->pszFilename, "https://")
841              || STARTS_WITH_CI(poOpenInfo->pszFilename, "jpip://"))
842             && EQUAL(pszExtension, "jp2") )
843         {
844             bIsJPIP = true;
845         }
846         else if( STARTS_WITH_CI(poOpenInfo->pszFilename, "J2K_SUBFILE:") )
847         {
848             try
849             {
850                 poRawInput = new subfile_source;
851                 poRawInput->open(poOpenInfo->pszFilename, bResilient,
852                                  bBuffered);
853                 poRawInput->seek(0);
854 
855                 poRawInput->read(&abySubfileHeader[0], 16);
856                 poRawInput->seek(0);
857             }
858             catch( ... )
859             {
860                 return nullptr;
861             }
862 
863             pabyHeader = abySubfileHeader.data();
864 
865             bIsSubfile = true;
866         }
867         else
868         {
869             return nullptr;
870         }
871     }
872     else
873     {
874         pabyHeader = poOpenInfo->pabyHeader;
875     }
876 
877     // If we think this should be access via vsil, then open it using
878     // subfile_source.  We do this if it does not seem to open normally
879     // or if we want to operate in resilient (sequential) mode.
880     VSIStatBuf sStat;
881     if( poRawInput == nullptr && !bIsJPIP &&
882         (bBuffered || bResilient ||
883          VSIStat(poOpenInfo->pszFilename, &sStat) != 0) )
884     {
885         try
886         {
887             poRawInput = new subfile_source;
888             poRawInput->open(poOpenInfo->pszFilename, bResilient, bBuffered);
889             poRawInput->seek(0);
890         }
891         catch( ... )
892         {
893             delete poRawInput;
894             return nullptr;
895         }
896     }
897 
898     // If the header is a JP2 header, mark this as a JP2 dataset.
899     if( pabyHeader && memcmp(pabyHeader, jp2_header, sizeof(jp2_header)) == 0 )
900         pszExtension = "jp2";
901 
902     // Try to open the file in a manner depending on the extension.
903     kdu_compressed_source *poInput = nullptr;
904     kdu_client *jpip_client = nullptr;
905     jp2_palette oJP2Palette;
906     jp2_channels oJP2Channels;
907 
908     jp2_family_src *family = nullptr;
909 
910     try
911     {
912         if( bIsJPIP )
913         {
914 #ifdef USE_JPIP
915             char *pszWrk =
916                 CPLStrdup(strstr(poOpenInfo->pszFilename, "://") + 3);
917             char *pszRequest = strstr(pszWrk, "/");
918 
919             if( pszRequest == NULL )
920             {
921                 CPLDebug("JP2KAK", "Failed to parse JPIP server and request.");
922                 CPLFree(pszWrk);
923                 return NULL;
924             }
925 
926             *(pszRequest++) = '\0';
927 
928             CPLDebug("JP2KAK", "server=%s, request=%s", pszWrk, pszRequest);
929 
930             CPLSleep(15.0);
931             jpip_client = new kdu_client;
932             jpip_client->connect(pszWrk, NULL, pszRequest, "http-tcp", "");
933 
934             CPLDebug("JP2KAK", "After connect()");
935 
936             bool bin0_complete = false;
937 
938             while( jpip_client->get_databin_length(KDU_META_DATABIN, 0, 0,
939                                                    &bin0_complete) <= 0
940                    || !bin0_complete )
941                 CPLSleep( 0.25 );
942 
943             family = new jp2_family_src;
944             family->open(jpip_client);
945 
946             // TODO(schwehr): Check for memory leaks.
947             jp2_source *jp2_src = new jp2_source;
948             jp2_src->open(family);
949             jp2_src->read_header();
950 
951             while( !jpip_client->is_idle() )
952                 CPLSleep(0.25);
953 
954             if( jpip_client->is_alive() )
955             {
956                 CPLDebug("JP2KAK", "connect() seems to be complete.");
957             }
958             else
959             {
960                 CPLDebug("JP2KAK", "connect() seems to have failed.");
961                 return NULL;
962             }
963 
964             oJP2Channels = jp2_src->access_channels();
965 
966             poInput = jp2_src;
967 #else
968             CPLError(CE_Failure, CPLE_OpenFailed,
969                      "JPIP Protocol not supported by GDAL with "
970                      "Kakadu 3.4 or on Unix.");
971             return nullptr;
972 #endif
973         }
974         else if( pszExtension != nullptr &&
975                  (EQUAL(pszExtension, "jp2") || EQUAL(pszExtension, "jpx")) )
976         {
977             family = new jp2_family_src;
978             if( poRawInput != nullptr )
979                 family->open(poRawInput);
980             else
981                 family->open(poOpenInfo->pszFilename, true);
982             jp2_source *jp2_src = new jp2_source;
983             poInput = jp2_src;
984             if( !jp2_src->open(family) || !jp2_src->read_header() )
985             {
986                 CPLDebug("JP2KAK", "Cannot read JP2 boxes");
987                 delete jp2_src;
988                 delete family;
989                 delete poRawInput;
990                 return nullptr;
991             }
992 
993             oJP2Palette = jp2_src->access_palette();
994             oJP2Channels = jp2_src->access_channels();
995 
996             jp2_colour oColors = jp2_src->access_colour();
997             if( oColors.get_space() != JP2_sRGB_SPACE &&
998                 oColors.get_space() != JP2_sLUM_SPACE )
999             {
1000                 CPLDebug("JP2KAK",
1001                          "Unusual ColorSpace=%d, not further interpreted.",
1002                          static_cast<int>(oColors.get_space()));
1003             }
1004         }
1005         else if( poRawInput == nullptr )
1006         {
1007             poInput = new kdu_simple_file_source(poOpenInfo->pszFilename);
1008         }
1009         else
1010         {
1011             poInput = poRawInput;
1012             poRawInput = nullptr;
1013         }
1014     }
1015     catch( ... )
1016     {
1017         CPLDebug("JP2KAK", "Trapped Kakadu exception.");
1018         delete family;
1019         delete poRawInput;
1020         delete poInput;
1021         return nullptr;
1022     }
1023 
1024     // Create a corresponding GDALDataset.
1025     JP2KAKDataset *poDS = nullptr;
1026 
1027     try
1028     {
1029         poDS = new JP2KAKDataset();
1030 
1031         poDS->poInput = poInput;
1032         poDS->poRawInput = poRawInput;
1033         poDS->family = family;
1034         poDS->oCodeStream.create(poInput);
1035         poDS->oCodeStream.set_persistent();
1036 
1037         poDS->bCached = bBuffered;
1038         poDS->bResilient = bResilient;
1039         poDS->bFussy = CPLTestBool(CPLGetConfigOption("JP2KAK_FUSSY", "NO"));
1040 
1041         if( poDS->bFussy )
1042             poDS->oCodeStream.set_fussy();
1043         if( poDS->bResilient )
1044             poDS->oCodeStream.set_resilient();
1045 
1046         poDS->jpip_client = jpip_client;
1047 
1048         // Get overall image size.
1049         poDS->oCodeStream.get_dims(0, poDS->dims);
1050 
1051         poDS->nRasterXSize = poDS->dims.size.x;
1052         poDS->nRasterYSize = poDS->dims.size.y;
1053 
1054         // Ensure that all the components have the same dimensions.  If
1055         // not, just process the first dimension.
1056         poDS->nBands = poDS->oCodeStream.get_num_components();
1057 
1058         if (poDS->nBands > 1 )
1059         {
1060             for( int iDim = 1; iDim < poDS->nBands; iDim++ )
1061             {
1062                 kdu_dims dim_this_comp;
1063 
1064                 poDS->oCodeStream.get_dims(iDim, dim_this_comp);
1065 
1066                 if( dim_this_comp != poDS->dims )
1067                 {
1068                     CPLError(CE_Warning, CPLE_AppDefined,
1069                              "Some components have mismatched dimensions, "
1070                              "ignoring all but first.");
1071                     poDS->nBands = 1;
1072                     break;
1073                 }
1074             }
1075         }
1076 
1077         // Setup the thread environment.
1078 
1079         int nNumThreads = atoi(CPLGetConfigOption("JP2KAK_THREADS", "-1"));
1080         if( nNumThreads == -1 )
1081             nNumThreads = kdu_get_num_processors() - 1;
1082         if( nNumThreads > 1024 )
1083             nNumThreads = 1024;
1084 
1085         if( nNumThreads > 0 )
1086         {
1087             poDS->poThreadEnv = new kdu_thread_env;
1088             poDS->poThreadEnv->create();
1089 
1090             for( int iThread=0; iThread < nNumThreads; iThread++ )
1091             {
1092                 if( !poDS->poThreadEnv->add_thread() )
1093                 {
1094                     CPLError(CE_Warning, CPLE_AppDefined,
1095                              "JP2KAK_THREADS: Unable to create thread.");
1096                     break;
1097                 }
1098             }
1099             CPLDebug("JP2KAK", "Using %d threads.", nNumThreads);
1100         }
1101         else
1102         {
1103             CPLDebug("JP2KAK", "Operating in singlethreaded mode.");
1104         }
1105 
1106         // Is this a file with poor internal navigation that will end
1107         // up using a great deal of memory if we use keep persistent
1108         // parsed information around?  (#3295)
1109         siz_params *siz = poDS->oCodeStream.access_siz();
1110         kdu_params *cod = siz->access_cluster(COD_params);
1111         bool use_precincts = false;
1112 
1113         cod->get(Cuse_precincts, 0, 0, use_precincts);
1114 
1115         const char *pszPersist = CPLGetConfigOption("JP2KAK_PERSIST", "AUTO");
1116         if( EQUAL(pszPersist, "AUTO") )
1117         {
1118             if( !use_precincts && !bIsJPIP &&
1119                 (poDS->nRasterXSize * static_cast<double>(poDS->nRasterYSize)) >
1120                     100000000.0 )
1121                 poDS->bPreferNPReads = true;
1122         }
1123         else
1124         {
1125             poDS->bPreferNPReads = !CPLTestBool(pszPersist);
1126         }
1127 
1128         CPLDebug("JP2KAK", "Cuse_precincts=%d, PreferNonPersistentReads=%d",
1129                  use_precincts ? 1 : 0, poDS->bPreferNPReads ? 1 : 0);
1130 
1131         //  Deduce some other info about the dataset.
1132         int order = 0;
1133 
1134         cod->get(Corder, 0, 0, order);
1135 
1136         // TODO(schwehr): Why not a switch statement?
1137         if( order == Corder_LRCP )
1138         {
1139             CPLDebug("JP2KAK", "order=LRCP");
1140             poDS->SetMetadataItem("Corder", "LRCP");
1141         }
1142         else if( order == Corder_RLCP )
1143         {
1144             CPLDebug("JP2KAK", "order=RLCP");
1145             poDS->SetMetadataItem("Corder", "RLCP");
1146         }
1147         else if( order == Corder_RPCL )
1148         {
1149             CPLDebug("JP2KAK", "order=RPCL");
1150             poDS->SetMetadataItem("Corder", "RPCL");
1151         }
1152         else if( order == Corder_PCRL )
1153         {
1154             CPLDebug("JP2KAK", "order=PCRL");
1155             poDS->SetMetadataItem("Corder", "PCRL");
1156         }
1157         else if( order == Corder_CPRL )
1158         {
1159             CPLDebug("JP2KAK", "order=CPRL");
1160             poDS->SetMetadataItem("Corder", "CPRL");
1161         }
1162         else
1163         {
1164             CPLDebug("JP2KAK", "order=%d, not recognized.", order);
1165         }
1166 
1167         poDS->bUseYCC = false;
1168         cod->get(Cycc, 0, 0, poDS->bUseYCC);
1169         if( poDS->bUseYCC )
1170             CPLDebug("JP2KAK", "ycc=true");
1171 
1172         // Find out how many resolutions levels are available.
1173         kdu_dims tile_indices;
1174         poDS->oCodeStream.get_valid_tiles(tile_indices);
1175 
1176         kdu_tile tile = poDS->oCodeStream.open_tile(tile_indices.pos);
1177         poDS->nResCount = tile.access_component(0).get_num_resolutions();
1178         tile.close();
1179 
1180         CPLDebug("JP2KAK", "nResCount=%d", poDS->nResCount);
1181 
1182         // Should we promote alpha channel to 8 bits?
1183         poDS->bPromoteTo8Bit =
1184             poDS->nBands == 4 &&
1185             poDS->oCodeStream.get_bit_depth(0) == 8 &&
1186             poDS->oCodeStream.get_bit_depth(1) == 8 &&
1187             poDS->oCodeStream.get_bit_depth(2) == 8 &&
1188             poDS->oCodeStream.get_bit_depth(3) == 1 &&
1189             CPLFetchBool(poOpenInfo->papszOpenOptions,
1190                          "1BIT_ALPHA_PROMOTION", true);
1191         if( poDS->bPromoteTo8Bit )
1192             CPLDebug("JP2KAK",
1193                      "Fourth (alpha) band is promoted from 1 bit to 8 bit");
1194 
1195         // Create band information objects.
1196         for( int iBand = 1; iBand <= poDS->nBands; iBand++ )
1197         {
1198             JP2KAKRasterBand *poBand = new JP2KAKRasterBand(
1199                 iBand, 0, poDS->oCodeStream, poDS->nResCount, jpip_client,
1200                 oJP2Channels, poDS);
1201 
1202             if( iBand == 1 && oJP2Palette.exists() )
1203                 poBand->ApplyPalette( oJP2Palette );
1204 
1205             poDS->SetBand( iBand, poBand );
1206         }
1207 
1208         // Look for supporting coordinate system information.
1209         if( poOpenInfo->nHeaderBytes != 0 )
1210         {
1211             poDS->LoadJP2Metadata(poOpenInfo);
1212         }
1213 
1214         // Establish our corresponding physical file.
1215         CPLString osPhysicalFilename = poOpenInfo->pszFilename;
1216 
1217         if( bIsSubfile ||
1218             STARTS_WITH_CI(poOpenInfo->pszFilename, "/vsisubfile/") )
1219         {
1220             const char *comma = strstr(poOpenInfo->pszFilename, ",");
1221             if( comma != nullptr )
1222                 osPhysicalFilename = comma + 1;
1223         }
1224 
1225 /* -------------------------------------------------------------------- */
1226 /*      Initialize any PAM information.                                 */
1227 /* -------------------------------------------------------------------- */
1228         poDS->SetDescription(poOpenInfo->pszFilename);
1229         if( !bIsSubfile )
1230             poDS->TryLoadXML();
1231         else
1232             poDS->nPamFlags |= GPF_NOSAVE;
1233 
1234         // Check for external overviews.
1235         poDS->oOvManager.Initialize(poDS, osPhysicalFilename);
1236 
1237         // Confirm the requested access is supported.
1238         if( poOpenInfo->eAccess == GA_Update )
1239         {
1240             delete poDS;
1241             CPLError(CE_Failure, CPLE_NotSupported,
1242                      "The JP2KAK driver does not support "
1243                      "update access to existing datasets.");
1244             return nullptr;
1245         }
1246 
1247         // Vector layers.
1248         if( poOpenInfo->nOpenFlags & GDAL_OF_VECTOR )
1249         {
1250             poDS->LoadVectorLayers(CPLFetchBool(poOpenInfo->papszOpenOptions,
1251                                                 "OPEN_REMOTE_GML", false));
1252 
1253             // If file opened in vector-only mode and there's no vector,
1254             // return.
1255             if( (poOpenInfo->nOpenFlags & GDAL_OF_RASTER) == 0 &&
1256                 poDS->GetLayerCount() == 0 )
1257             {
1258                 delete poDS;
1259                 return nullptr;
1260             }
1261         }
1262 
1263         return poDS;
1264     }
1265     catch( ... )
1266     {
1267         CPLDebug("JP2KAK", "JP2KAKDataset::Open() - caught exception.");
1268         if( poDS != nullptr )
1269             delete poDS;
1270 
1271         return nullptr;
1272     }
1273 }
1274 
1275 /************************************************************************/
1276 /*                           DirectRasterIO()                           */
1277 /************************************************************************/
1278 
1279 CPLErr
1280 JP2KAKDataset::DirectRasterIO( GDALRWFlag /* eRWFlag */,
1281                                int nXOff, int nYOff, int nXSize, int nYSize,
1282                                void * pData, int nBufXSize, int nBufYSize,
1283                                GDALDataType eBufType,
1284                                int nBandCount, int *panBandMap,
1285                                GSpacing nPixelSpace, GSpacing nLineSpace,
1286                                GSpacing nBandSpace,
1287                                GDALRasterIOExtraArg* psExtraArg)
1288 
1289 {
1290     CPLAssert(eBufType == GDT_Byte || eBufType == GDT_Int16 ||
1291               eBufType == GDT_UInt16);
1292 
1293     kdu_codestream *poCodeStream = &oCodeStream;
1294     const char *pszPersistency = "";
1295 
1296     // Do we want to do this non-persistently?  If so, we need to
1297     // open the file, and establish a local codestream.
1298     subfile_source subfile_src;
1299     jp2_source wrk_jp2_src;
1300     jp2_family_src wrk_family;
1301     kdu_codestream oWCodeStream;
1302 
1303     if( bPreferNPReads )
1304     {
1305         subfile_src.open(GetDescription(), bResilient, bCached);
1306 
1307         if( family != nullptr )
1308         {
1309             wrk_family.open(&subfile_src);
1310             wrk_jp2_src.open(&wrk_family);
1311             wrk_jp2_src.read_header();
1312 
1313             oWCodeStream.create(&wrk_jp2_src, poThreadEnv);
1314         }
1315         else
1316         {
1317             oWCodeStream.create(&subfile_src, poThreadEnv);
1318         }
1319 
1320         if( bFussy )
1321             oWCodeStream.set_fussy();
1322         if( bResilient )
1323             oWCodeStream.set_resilient();
1324 
1325         poCodeStream = &oWCodeStream;
1326 
1327         pszPersistency = "(non-persistent)";
1328     }
1329 
1330     // Select optimal resolution level.
1331     int nDiscardLevels = 0;
1332     int nResMult = 1;
1333 
1334     if( AreOverviewsEnabled() )
1335     {
1336         while( nDiscardLevels < nResCount - 1 &&
1337                nBufXSize * nResMult * 2 < nXSize * 1.01 &&
1338                nBufYSize * nResMult * 2 < nYSize * 1.01 )
1339         {
1340             nDiscardLevels++;
1341             nResMult *= 2;
1342         }
1343     }
1344 
1345     // Prepare component indices list.
1346     CPLErr eErr = CE_None;
1347 
1348     std::vector<int> component_indices(nBandCount);
1349     std::vector<int> stripe_heights(nBandCount);
1350     std::vector<int> sample_offsets(nBandCount);
1351     std::vector<int> sample_gaps(nBandCount);
1352     std::vector<int> row_gaps(nBandCount);
1353     std::vector<int> precisions(nBandCount);
1354 
1355     // std::vector<bool> is Unfortunately a specialized implementation such as &v[0] doesn't work
1356     // https://codereview.stackexchange.com/questions/241629/stdvectorbool-workaround-in-c
1357     class vector_safe_bool {
1358         bool value;
1359     public:
1360         vector_safe_bool() = default;
1361         // cppcheck-suppress noExplicitConstructor
1362         vector_safe_bool(bool b) : value{b} {}
1363 
1364         bool *operator&() noexcept { return &value; }
1365         const bool *operator&() const noexcept { return &value; }
1366 
1367         operator const bool &() const noexcept { return value; }
1368         operator bool &() noexcept { return value; }
1369     };
1370     std::vector<vector_safe_bool> is_signed(nBandCount);
1371 
1372     for( int i = 0; i < nBandCount; i++ )
1373         component_indices[i] = panBandMap[i] - 1;
1374 
1375     // Setup a ROI matching the block requested, and select desired
1376     // bands (components).
1377     try
1378     {
1379         poCodeStream->apply_input_restrictions(0, 0, nDiscardLevels, 0, nullptr);
1380         kdu_dims l_dims;
1381         poCodeStream->get_dims(0, l_dims);
1382         const int nOvrCanvasXSize = l_dims.pos.x + l_dims.size.x;
1383         const int nOvrCanvasYSize = l_dims.pos.y + l_dims.size.y;
1384 
1385         l_dims.pos.x = l_dims.pos.x + nXOff / nResMult;
1386         l_dims.pos.y = l_dims.pos.y + nYOff / nResMult;
1387         l_dims.size.x = nXSize / nResMult;
1388         l_dims.size.y = nYSize / nResMult;
1389 
1390         // Check if rounding helps detecting when data is being requested
1391         // exactly at the current resolution.
1392         if( nBufXSize != l_dims.size.x &&
1393             static_cast<int>(0.5 + static_cast<double>(nXSize) / nResMult)
1394             == nBufXSize )
1395         {
1396             l_dims.size.x = nBufXSize;
1397         }
1398         if( nBufYSize != l_dims.size.y &&
1399             static_cast<int>(0.5 + static_cast<double>(nYSize) / nResMult)
1400             == nBufYSize )
1401         {
1402             l_dims.size.y = nBufYSize;
1403         }
1404         if( l_dims.pos.x + l_dims.size.x > nOvrCanvasXSize )
1405             l_dims.size.x = nOvrCanvasXSize - l_dims.pos.x;
1406         if( l_dims.pos.y + l_dims.size.y > nOvrCanvasYSize )
1407             l_dims.size.y = nOvrCanvasYSize - l_dims.pos.y;
1408 
1409         kdu_dims l_dims_roi;
1410 
1411         poCodeStream->map_region(0, l_dims, l_dims_roi);
1412         poCodeStream->apply_input_restrictions(nBandCount, component_indices.data(),
1413                                                nDiscardLevels, 0, &l_dims_roi,
1414                                                KDU_WANT_OUTPUT_COMPONENTS);
1415 
1416         // Special case where the data is being requested exactly at
1417         // this resolution.  Avoid any extra sampling pass.
1418         if( nBufXSize == l_dims.size.x && nBufYSize == l_dims.size.y &&
1419             (nBandCount - 1) * nBandSpace / GDALGetDataTypeSizeBytes(eBufType) < INT_MAX )
1420         {
1421             kdu_stripe_decompressor decompressor;
1422             decompressor.start(*poCodeStream, false, false, poThreadEnv);
1423 
1424             CPLDebug("JP2KAK", "DirectRasterIO() for %d,%d,%d,%d -> %dx%d "
1425                      "(no intermediate) %s",
1426                      nXOff, nYOff, nXSize, nYSize, nBufXSize, nBufYSize,
1427                      pszPersistency);
1428 
1429             for( int i = 0; i < nBandCount; i++ )
1430             {
1431                 stripe_heights[i] = l_dims.size.y;
1432                 precisions[i] = poCodeStream->get_bit_depth(i);
1433                 if( eBufType == GDT_Byte )
1434                 {
1435                     is_signed[i] = false;
1436                     sample_offsets[i] = static_cast<int>(i * nBandSpace);
1437                     sample_gaps[i] = static_cast<int>(nPixelSpace);
1438                     row_gaps[i] = static_cast<int>(nLineSpace);
1439                 }
1440                 else if( eBufType == GDT_Int16 )
1441                 {
1442                     is_signed[i] = true;
1443                     sample_offsets[i] = static_cast<int>(i * nBandSpace / 2);
1444                     sample_gaps[i] = static_cast<int>(nPixelSpace / 2);
1445                     row_gaps[i] = static_cast<int>(nLineSpace / 2);
1446                 }
1447                 else if( eBufType == GDT_UInt16 )
1448                 {
1449                     is_signed[i] = false;
1450                     sample_offsets[i] = static_cast<int>(i * nBandSpace / 2);
1451                     sample_gaps[i] = static_cast<int>(nPixelSpace / 2);
1452                     row_gaps[i] = static_cast<int>(nLineSpace / 2);
1453                     /* Introduced in r25136 with an unrelated commit message.
1454                     Reverted per ticket #5328
1455                     if( precisions[i] == 12 )
1456                     {
1457                       CPLDebug( "JP2KAK", "16bit extend 12 bit data." );
1458                       precisions[i] = 16;
1459                     }*/
1460                 }
1461             }
1462 
1463             if( eBufType == GDT_Byte )
1464                 decompressor.pull_stripe(
1465                     static_cast<kdu_byte *>(pData), &stripe_heights[0],
1466                     &sample_offsets[0], &sample_gaps[0], &row_gaps[0],
1467                     &precisions[0]);
1468             else
1469                 decompressor.pull_stripe(
1470                     static_cast<kdu_int16 *>(pData), &stripe_heights[0],
1471                     &sample_offsets[0], &sample_gaps[0], &row_gaps[0],
1472                     &precisions[0], &is_signed[0]);
1473             decompressor.finish();
1474         }
1475         else
1476         {
1477             // More general case - first pull into working buffer.
1478 
1479             const int nDataTypeSize = GDALGetDataTypeSizeBytes(eBufType);
1480             GByte *pabyIntermediate = static_cast<GByte *>(VSI_MALLOC3_VERBOSE(
1481                 l_dims.size.x, l_dims.size.y, nDataTypeSize * nBandCount));
1482             if( pabyIntermediate == nullptr )
1483             {
1484                 return CE_Failure;
1485             }
1486 
1487             CPLDebug("JP2KAK",
1488                      "DirectRasterIO() for %d,%d,%d,%d -> %dx%d -> %dx%d %s",
1489                      nXOff, nYOff, nXSize, nYSize, l_dims.size.x, l_dims.size.y,
1490                      nBufXSize, nBufYSize, pszPersistency);
1491 
1492             kdu_stripe_decompressor decompressor;
1493             decompressor.start(*poCodeStream, false, false, poThreadEnv);
1494 
1495             for( int i = 0; i < nBandCount; i++ )
1496             {
1497                 stripe_heights[i] = l_dims.size.y;
1498                 precisions[i] = poCodeStream->get_bit_depth(i);
1499 
1500                 if( eBufType == GDT_Int16 || eBufType == GDT_UInt16 )
1501                 {
1502                     is_signed[i] = eBufType == GDT_Int16;
1503                 }
1504             }
1505 
1506             if( eBufType == GDT_Byte )
1507                 decompressor.pull_stripe(
1508                     reinterpret_cast<kdu_byte *>(pabyIntermediate),
1509                     &stripe_heights[0], nullptr, nullptr, nullptr, &precisions[0]);
1510             else
1511                 decompressor.pull_stripe(
1512                     reinterpret_cast<kdu_int16 *>(pabyIntermediate),
1513                     &stripe_heights[0], nullptr, nullptr, nullptr, &precisions[0],
1514                     &is_signed[0]);
1515 
1516             decompressor.finish();
1517 
1518             if( psExtraArg->eResampleAlg == GRIORA_NearestNeighbour )
1519             {
1520                 // Then resample (normally downsample) from the intermediate
1521                 // buffer into the final buffer in the desired output layout.
1522                 const double dfYRatio =
1523                     l_dims.size.y / static_cast<double>(nBufYSize);
1524                 const double dfXRatio =
1525                     l_dims.size.x / static_cast<double>(nBufXSize);
1526 
1527                 for( int iY = 0; iY < nBufYSize; iY++ )
1528                 {
1529                     const int iSrcY = std::min(
1530                         static_cast<int>(floor((iY + 0.5) * dfYRatio)),
1531                         l_dims.size.y - 1);
1532 
1533                     for( int iX = 0; iX < nBufXSize; iX++ )
1534                     {
1535                         const int iSrcX = std::min(
1536                             static_cast<int>(floor((iX + 0.5) * dfXRatio)),
1537                             l_dims.size.x - 1);
1538 
1539                         for( int i = 0; i < nBandCount; i++ )
1540                         {
1541                             // TODO(schwehr): Cleanup this block.
1542                             if( eBufType == GDT_Byte )
1543                                 ((GByte *) pData)[iX*nPixelSpace
1544                                                 + iY*nLineSpace
1545                                                 + i*nBandSpace] =
1546                                     pabyIntermediate[iSrcX*nBandCount
1547                                                     + static_cast<GPtrDiff_t>(iSrcY)*l_dims.size.x*nBandCount
1548                                                     + i];
1549                             else if( eBufType == GDT_Int16
1550                                     || eBufType == GDT_UInt16 )
1551                                 ((GUInt16 *) pData)[iX*nPixelSpace/2
1552                                                 + iY*nLineSpace/2
1553                                                 + i*nBandSpace/2] =
1554                                     ((GUInt16 *)pabyIntermediate)[
1555                                         iSrcX*nBandCount
1556                                         + static_cast<GPtrDiff_t>(iSrcY)*l_dims.size.x*nBandCount
1557                                         + i];
1558                         }
1559                     }
1560                 }
1561             }
1562             else
1563             {
1564                 // Create a MEM dataset that wraps the input buffer.
1565                 GDALDataset* poMEMDS = MEMDataset::Create( "", l_dims.size.x,
1566                                                            l_dims.size.y, 0,
1567                                                            eBufType, nullptr);
1568                 char szBuffer[64] = { '\0' };
1569                 int nRet = 0;
1570 
1571                 for( int i = 0; i < nBandCount; i++ )
1572                 {
1573 
1574                     nRet = CPLPrintPointer(
1575                         szBuffer, pabyIntermediate + i * nDataTypeSize,
1576                         sizeof(szBuffer) );
1577                     szBuffer[nRet] = 0;
1578                     char **papszOptions =
1579                         CSLSetNameValue(nullptr, "DATAPOINTER", szBuffer);
1580 
1581                     papszOptions = CSLSetNameValue(papszOptions, "PIXELOFFSET",
1582                         CPLSPrintf(
1583                             CPL_FRMT_GIB,
1584                             static_cast<GIntBig>(nDataTypeSize) * nBandCount));
1585 
1586                     papszOptions = CSLSetNameValue(papszOptions, "LINEOFFSET",
1587                         CPLSPrintf(
1588                             CPL_FRMT_GIB,
1589                             static_cast<GIntBig>(nDataTypeSize)
1590                             * nBandCount * l_dims.size.x) );
1591 
1592                     poMEMDS->AddBand(eBufType, papszOptions);
1593                     CSLDestroy(papszOptions);
1594 
1595                     const char *pszNBITS =
1596                         GetRasterBand(i + 1)
1597                             ->GetMetadataItem("NBITS", "IMAGE_STRUCTURE");
1598                     if( pszNBITS )
1599                         poMEMDS->GetRasterBand(i + 1)->SetMetadataItem(
1600                             "NBITS", pszNBITS, "IMAGE_STRUCTURE");
1601                 }
1602 
1603                 GDALRasterIOExtraArg sExtraArgTmp;
1604                 INIT_RASTERIO_EXTRA_ARG(sExtraArgTmp);
1605                 sExtraArgTmp.eResampleAlg = psExtraArg->eResampleAlg;
1606 
1607                 CPL_IGNORE_RET_VAL(
1608                     poMEMDS->RasterIO(
1609                         GF_Read, 0, 0, l_dims.size.x, l_dims.size.y,
1610                         pData, nBufXSize, nBufYSize,
1611                         eBufType,
1612                         nBandCount, nullptr,
1613                         nPixelSpace, nLineSpace, nBandSpace,
1614                         &sExtraArgTmp ) );
1615 
1616                 GDALClose(poMEMDS);
1617             }
1618 
1619             CPLFree(pabyIntermediate);
1620         }
1621     }
1622     catch( ... )
1623     {
1624         // Catch internal Kakadu errors.
1625         eErr = CE_Failure;
1626     }
1627 
1628     // 1-bit alpha promotion.
1629     if( nBandCount == 4 && bPromoteTo8Bit )
1630     {
1631         for( int j = 0; j < nBufYSize; j++ )
1632         {
1633             for( int i = 0; i < nBufXSize; i++ )
1634             {
1635                 static_cast<GByte *>(
1636                     pData)[j * nLineSpace + i * nPixelSpace + 3 * nBandSpace] *=
1637                     255;
1638             }
1639         }
1640     }
1641 /* -------------------------------------------------------------------- */
1642 /*      Cleanup                                                         */
1643 /* -------------------------------------------------------------------- */
1644     if( poCodeStream == &oWCodeStream )
1645     {
1646         oWCodeStream.destroy();
1647         wrk_jp2_src.close();
1648         wrk_family.close();
1649         subfile_src.close();
1650     }
1651 
1652     return eErr;
1653 }
1654 
1655 /************************************************************************/
1656 /*                           TestUseBlockIO()                           */
1657 /*                                                                      */
1658 /*      Check whether we should use blocked IO (true) or direct io      */
1659 /*      (FALSE) for a given request configuration and environment.      */
1660 /************************************************************************/
1661 
1662 bool
1663 JP2KAKDataset::TestUseBlockIO( int nXOff, int nYOff, int nXSize, int nYSize,
1664                                int nBufXSize, int nBufYSize,
1665                                GDALDataType eDataType,
1666                                int nBandCount, int *panBandList )
1667 
1668 {
1669     // Due to limitations in DirectRasterIO() we can only handle
1670     // 8bit and with no duplicates in the band list.
1671     if( eDataType != GetRasterBand(1)->GetRasterDataType()
1672         || ( eDataType != GDT_Byte
1673              && eDataType != GDT_Int16
1674              && eDataType != GDT_UInt16 ) )
1675         return true;
1676 
1677     for( int i = 0; i < nBandCount; i++ )
1678     {
1679         for( int j = i + 1; j < nBandCount; j++ )
1680             if( panBandList[j] == panBandList[i] )
1681                 return true;
1682     }
1683 
1684     // If we have external overviews built and they could be used to satisfy
1685     // this request, we will avoid DirectRasterIO() which would ignore them.
1686     if( GetRasterCount() == 0 )
1687         return true;
1688 
1689     JP2KAKRasterBand *poWrkBand =
1690         dynamic_cast<JP2KAKRasterBand *>(GetRasterBand(1));
1691     if( poWrkBand == nullptr )
1692     {
1693         CPLError( CE_Fatal, CPLE_AppDefined, "Dynamic cast failed" );
1694         return false;
1695     }
1696     if( poWrkBand->HasExternalOverviews() )
1697     {
1698         int nXOff2 = nXOff;
1699         int nYOff2 = nYOff;
1700         int nXSize2 = nXSize;
1701         int nYSize2 = nYSize;
1702 
1703         const int nOverview =
1704             GDALBandGetBestOverviewLevel2(poWrkBand, nXOff2, nYOff2, nXSize2,
1705                                           nYSize2, nBufXSize, nBufYSize, nullptr);
1706         if( nOverview >= 0 )
1707             return true;
1708     }
1709 
1710     // The rest of the rules are io strategy stuff and configuration checks.
1711     bool bUseBlockedIO = bForceCachedIO;
1712 
1713     if( nYSize == 1 || nXSize * static_cast<double>(nYSize) < 100.0 )
1714         bUseBlockedIO = true;
1715 
1716     if( nBufYSize == 1 || nBufXSize * static_cast<double>(nBufYSize) < 100.0 )
1717         bUseBlockedIO = true;
1718 
1719     if( strlen(CPLGetConfigOption("GDAL_ONE_BIG_READ", "")) > 0 )
1720         bUseBlockedIO =
1721             !CPLTestBool(CPLGetConfigOption("GDAL_ONE_BIG_READ", ""));
1722 
1723     return bUseBlockedIO;
1724 }
1725 
1726 /************************************************************************/
1727 /*                             IRasterIO()                              */
1728 /************************************************************************/
1729 
1730 CPLErr JP2KAKDataset::IRasterIO( GDALRWFlag eRWFlag,
1731                                  int nXOff, int nYOff, int nXSize, int nYSize,
1732                                  void * pData, int nBufXSize, int nBufYSize,
1733                                  GDALDataType eBufType,
1734                                  int nBandCount, int *panBandMap,
1735                                  GSpacing nPixelSpace, GSpacing nLineSpace,
1736                                  GSpacing nBandSpace,
1737                                  GDALRasterIOExtraArg* psExtraArg)
1738 
1739 {
1740     // We need various criteria to skip out to block based methods.
1741     if( TestUseBlockIO(nXOff, nYOff, nXSize, nYSize, nBufXSize, nBufYSize,
1742                        eBufType, nBandCount, panBandMap) )
1743         return GDALPamDataset::IRasterIO(
1744             eRWFlag, nXOff, nYOff, nXSize, nYSize,
1745             pData, nBufXSize, nBufYSize, eBufType,
1746             nBandCount, panBandMap, nPixelSpace, nLineSpace, nBandSpace,
1747             psExtraArg );
1748 
1749     return DirectRasterIO(
1750         eRWFlag, nXOff, nYOff, nXSize, nYSize, pData, nBufXSize, nBufYSize,
1751         eBufType, nBandCount, panBandMap, nPixelSpace, nLineSpace,
1752         nBandSpace, psExtraArg );
1753 }
1754 
1755 /************************************************************************/
1756 /*                           JP2KAKWriteBox()                           */
1757 /*                                                                      */
1758 /*      Write out the passed box and delete it.                         */
1759 /************************************************************************/
1760 
1761 static void JP2KAKWriteBox( jp2_family_tgt *jp2_family, GDALJP2Box *poBox )
1762 
1763 {
1764     if( poBox == nullptr )
1765         return;
1766 
1767     jp2_output_box jp2_out;
1768 
1769     GUInt32 nBoxType = 0;
1770     memcpy(&nBoxType, poBox->GetType(), sizeof(nBoxType));
1771     CPL_MSBPTR32(&nBoxType);
1772 
1773     int length = static_cast<int>(poBox->GetDataLength());
1774 
1775     // Write to a box on the JP2 file.
1776     jp2_out.open(jp2_family, nBoxType);
1777     jp2_out.set_target_size(length);
1778     jp2_out.write(const_cast<kdu_byte *>(poBox->GetWritableData()),
1779                    length);
1780     jp2_out.close();
1781 
1782     delete poBox;
1783 }
1784 
1785 /************************************************************************/
1786 /*                     JP2KAKCreateCopy_WriteTile()                     */
1787 /************************************************************************/
1788 
1789 static bool
1790 JP2KAKCreateCopy_WriteTile( GDALDataset *poSrcDS, kdu_tile &oTile,
1791                             kdu_roi_image *poROIImage,
1792                             int nXOff, int nYOff, int nXSize, int nYSize,
1793                             bool bReversible, int nBits, GDALDataType eType,
1794                             kdu_codestream &oCodeStream, int bFlushEnabled,
1795                             kdu_long *layer_bytes, int layer_count,
1796                             GDALProgressFunc pfnProgress, void * pProgressData,
1797                             bool bComseg )
1798 
1799 {
1800     // Create one big tile, a compressing engine, and a line buffer for each
1801     // component.
1802     const int num_components = oTile.get_num_components();
1803     kdu_push_ifc *engines = new kdu_push_ifc[num_components];
1804     kdu_line_buf *lines = new kdu_line_buf[num_components];
1805     kdu_sample_allocator allocator;
1806 
1807     // Ticket #4050 patch: Use a 32 bits kdu_line_buf for GDT_UInt16 reversible
1808     // compression.
1809     // TODO: Test for GDT_UInt16?
1810     bool bUseShorts = bReversible;
1811     if( eType == GDT_UInt16 && bReversible )
1812         bUseShorts = false;
1813 
1814     for( int c = 0; c < num_components; c++ )
1815     {
1816         kdu_resolution res = oTile.access_component(c).access_resolution();
1817         kdu_roi_node *roi_node = nullptr;
1818 
1819         if( poROIImage != nullptr )
1820         {
1821             kdu_dims dims;
1822 
1823             res.get_dims(dims);
1824             roi_node = poROIImage->acquire_node(c, dims);
1825         }
1826 #if KDU_MAJOR_VERSION >= 7
1827         lines[c].pre_create(&allocator, nXSize, bReversible, bUseShorts, 0, 0);
1828 #else
1829         lines[c].pre_create(&allocator, nXSize, bReversible, bUseShorts);
1830 #endif
1831         engines[c] = kdu_analysis(res, &allocator, bUseShorts, 1.0F, roi_node);
1832     }
1833 
1834     try
1835     {
1836 #if KDU_MAJOR_VERSION > 7 || (KDU_MAJOR_VERSION == 7 && (KDU_MINOR_VERSION > 3 || KDU_MINOR_VERSION == 3 && KDU_PATCH_VERSION >= 1))
1837         allocator.finalize(oCodeStream);
1838 #else
1839         allocator.finalize();
1840 #endif
1841 
1842         for( int c = 0; c < num_components; c++ )
1843             lines[c].create();
1844     }
1845     catch( ... )
1846     {
1847         // TODO(schwehr): Should this block do any cleanup?
1848         CPLError(CE_Failure, CPLE_AppDefined,
1849                  "allocate.finalize() failed, likely out of memory for "
1850                  "compression information.");
1851         return false;
1852     }
1853 
1854     // Write whole image.  Write 1024 lines of each component, then
1855     // go back to the first, and do again.  This gives the rate
1856     // computing machine all components to make good estimates.
1857     int iLinesWritten = 0;
1858 
1859     // TODO(schwehr): Making pabyBuffer void* should simplify the casting below.
1860     GByte *pabyBuffer = static_cast<GByte *>(
1861         CPLMalloc(nXSize * GDALGetDataTypeSizeBytes(eType)));
1862 
1863     CPLAssert(!oTile.get_ycc());
1864 
1865     bool bRet = true;
1866     for( int iLine = 0; iLine < nYSize && bRet; iLine += TILE_CHUNK_SIZE )
1867     {
1868         for( int c = 0; c < num_components && bRet; c++ )
1869         {
1870             GDALRasterBand *poBand = poSrcDS->GetRasterBand(c + 1);
1871 
1872             for( int iSubline = iLine;
1873                  iSubline < iLine + TILE_CHUNK_SIZE && iSubline < nYSize;
1874                  iSubline++ )
1875             {
1876                 if( poBand->RasterIO( GF_Read,
1877                                       nXOff, nYOff+iSubline, nXSize, 1,
1878                                       pabyBuffer, nXSize, 1, eType,
1879                                       0, 0, nullptr ) == CE_Failure )
1880                 {
1881                     bRet = false;
1882                     break;
1883                 }
1884 
1885                 if( bReversible && eType == GDT_Byte )
1886                 {
1887                     kdu_sample16 *dest = lines[c].get_buf16();
1888                     kdu_byte *sp = pabyBuffer;
1889                     const kdu_int16 nOffset = 1 << (nBits - 1);
1890 
1891                     for( int n = nXSize; n > 0; n--, dest++, sp++ )
1892                         dest->ival = ((kdu_int16)(*sp)) - nOffset;
1893                 }
1894                 else if( bReversible && eType == GDT_Int16 )
1895                 {
1896                     kdu_sample16 *dest = lines[c].get_buf16();
1897                     GInt16 *sp = reinterpret_cast<GInt16 *>(pabyBuffer);
1898 
1899                     for( int n = nXSize; n > 0; n--, dest++, sp++ )
1900                         dest->ival = *sp;
1901                 }
1902                 else if( bReversible && eType == GDT_UInt16 )
1903                 {
1904                     // Ticket #4050 patch : use a 32 bits kdu_line_buf for
1905                     // GDT_UInt16 reversible compression.
1906                     kdu_sample32 *dest = lines[c].get_buf32();
1907                     GUInt16 *sp = reinterpret_cast<GUInt16 *>(pabyBuffer);
1908                     const kdu_int32 nOffset = 1 << (nBits - 1);
1909 
1910                     for( int n=nXSize; n > 0; n--, dest++, sp++ )
1911                         dest->ival = static_cast<kdu_int32>(*sp) - nOffset;
1912                 }
1913                 else if( eType == GDT_Byte )
1914                 {
1915                     kdu_sample32 *dest = lines[c].get_buf32();
1916                     kdu_byte *sp = pabyBuffer;
1917                     const int nOffset = 1 << (nBits - 1);
1918                     const float fScale = 1.0f / (1 << nBits);
1919 
1920                     for (int n = nXSize; n > 0; n--, dest++, sp++)
1921                         dest->fval =
1922                             (static_cast<kdu_int16>(*sp) - nOffset) * fScale;
1923                 }
1924                 else if( eType == GDT_Int16 )
1925                 {
1926                     kdu_sample32 *dest = lines[c].get_buf32();
1927                     GInt16 *sp = reinterpret_cast<GInt16 *>(pabyBuffer);
1928                     const float fScale = 1.0f / (1 << nBits);
1929 
1930                     for( int n = nXSize; n > 0; n--, dest++, sp++ )
1931                         dest->fval = static_cast<kdu_int16>(*sp) * fScale;
1932                 }
1933                 else if( eType == GDT_UInt16 )
1934                 {
1935                     kdu_sample32 *dest = lines[c].get_buf32();
1936                     GUInt16 *sp = reinterpret_cast<GUInt16 *>(pabyBuffer);
1937                     const int nOffset = 1 << (nBits - 1);
1938                     const float fScale = 1.0f / (1 << nBits);
1939 
1940                     for( int n = nXSize; n > 0; n--, dest++, sp++ )
1941                         dest->fval = (static_cast<int>(*sp) - nOffset) * fScale;
1942                 }
1943                 else if( eType == GDT_Float32 )
1944                 {
1945                     kdu_sample32 *dest = lines[c].get_buf32();
1946                     float *sp = reinterpret_cast<float *>(pabyBuffer);
1947 
1948                     for( int n = nXSize; n > 0; n--, dest++, sp++ )
1949                         dest->fval = *sp;  // Scale it?
1950                 }
1951 
1952 #if KDU_MAJOR_VERSION >= 7
1953                 engines[c].push(lines[c]);
1954 #else
1955                 engines[c].push(lines[c], true);
1956 #endif
1957                 iLinesWritten++;
1958 
1959                 if( !pfnProgress(iLinesWritten / static_cast<double>(
1960                                                      num_components * nYSize),
1961                                  nullptr, pProgressData) )
1962                 {
1963                     bRet = false;
1964                     break;
1965                 }
1966             }
1967         }
1968         if( !bRet )
1969             break;
1970 
1971         if( oCodeStream.ready_for_flush() && bFlushEnabled )
1972         {
1973             CPLDebug("JP2KAK", "Calling oCodeStream.flush() at line %d",
1974                      std::min(nYSize, iLine + TILE_CHUNK_SIZE));
1975             try
1976             {
1977                 oCodeStream.flush(layer_bytes, layer_count, nullptr,
1978                                   true, bComseg);
1979             }
1980             catch(...)
1981             {
1982                 bRet = false;
1983             }
1984         }
1985         else if( bFlushEnabled )
1986         {
1987             CPLDebug("JP2KAK", "read_for_flush() is false at line %d.", iLine);
1988         }
1989     }
1990 
1991     for( int c = 0; c < num_components; c++ )
1992         engines[c].destroy();
1993 
1994     delete[] engines;
1995     delete[] lines;
1996 
1997     CPLFree(pabyBuffer);
1998 
1999     if( poROIImage != nullptr )
2000         delete poROIImage;
2001 
2002     return bRet;
2003     // For some reason cppcheck thinks that engines and lines are leaking.
2004     // cppcheck-suppress memleak
2005 }
2006 
2007 /************************************************************************/
2008 /*                             CreateCopy()                             */
2009 /************************************************************************/
2010 
2011 static GDALDataset *
2012 JP2KAKCreateCopy( const char * pszFilename, GDALDataset *poSrcDS,
2013                   int bStrict, char ** papszOptions,
2014                   GDALProgressFunc pfnProgress, void * pProgressData )
2015 
2016 {
2017     if( poSrcDS->GetRasterCount() == 0 )
2018     {
2019         CPLError(CE_Failure, CPLE_AppDefined,
2020                  "Creating zero band files not supported by JP2KAK driver.");
2021         return nullptr;
2022     }
2023 
2024     // Initialize Kakadu warning/error reporting subsystem.
2025     if( !kakadu_initialized )
2026     {
2027         kakadu_initialized = true;
2028 
2029         kdu_cpl_error_message oErrHandler(CE_Failure);
2030         kdu_cpl_error_message oWarningHandler(CE_Warning);
2031 
2032         kdu_customize_warnings(new kdu_cpl_error_message(CE_Warning));
2033         kdu_customize_errors(new kdu_cpl_error_message(CE_Failure));
2034     }
2035 
2036     // What data type should we use?  We assume all datatypes match
2037     // the first band.
2038     GDALRasterBand *poPrototypeBand = poSrcDS->GetRasterBand(1);
2039 
2040     GDALDataType eType = poPrototypeBand->GetRasterDataType();
2041     if( eType != GDT_Byte && eType != GDT_Int16 && eType != GDT_UInt16 &&
2042         eType != GDT_Float32 )
2043     {
2044         if( bStrict )
2045         {
2046             CPLError(CE_Failure, CPLE_AppDefined,
2047                      "JP2KAK (JPEG2000) driver does not support data type %s.",
2048                      GDALGetDataTypeName(eType));
2049             return nullptr;
2050         }
2051 
2052         CPLError(CE_Warning, CPLE_AppDefined,
2053                  "JP2KAK (JPEG2000) driver does not support data type %s, "
2054                  "forcing to Float32.",
2055                  GDALGetDataTypeName(eType));
2056 
2057         eType = GDT_Float32;
2058     }
2059 
2060     // Do we want to write a pseudo-colored image?
2061     const int bHaveCT =
2062         poPrototypeBand->GetColorTable() != nullptr &&
2063         poSrcDS->GetRasterCount() == 1;
2064 
2065     // How many layers?
2066     int layer_count = 12;
2067 
2068     if( CSLFetchNameValue(papszOptions,"LAYERS") != nullptr )
2069         layer_count = atoi(CSLFetchNameValue(papszOptions, "LAYERS"));
2070     else if( CSLFetchNameValue(papszOptions,"Clayers") != nullptr )
2071         layer_count = atoi(CSLFetchNameValue(papszOptions, "Clayers"));
2072 
2073     // Establish how many bytes of data we want for each layer.
2074     // We take the quality as a percentage, so if QUALITY of 50 is
2075     // selected, we will set the base layer to 50% the default size.
2076     // We let the other layers be computed internally.
2077 
2078     const double dfQuality =
2079         CSLFetchNameValue(papszOptions,"QUALITY") != nullptr
2080         ? CPLAtof(CSLFetchNameValue(papszOptions, "QUALITY"))
2081         : 20.0;
2082 
2083     if( dfQuality < 0.01 || dfQuality > 100.0 )
2084     {
2085         CPLError(CE_Failure, CPLE_IllegalArg,
2086                  "QUALITY=%s is not a legal value in the range 0.01-100.",
2087                  CSLFetchNameValue(papszOptions, "QUALITY"));
2088         return nullptr;
2089     }
2090 
2091     std::vector<kdu_long> layer_bytes;
2092     layer_bytes.resize(layer_count);
2093 
2094     const int nXSize = poSrcDS->GetRasterXSize();
2095     const int nYSize = poSrcDS->GetRasterYSize();
2096 
2097     bool bReversible = false;
2098 
2099     if( dfQuality < 99.5 )
2100     {
2101         double dfLayerBytes =
2102             (nXSize * static_cast<double>(nYSize) * dfQuality / 100.0)
2103             * GDALGetDataTypeSizeBytes(eType)
2104             * GDALGetRasterCount(poSrcDS);
2105 
2106         if( dfLayerBytes > 2000000000.0 && sizeof(kdu_long) == 4 )
2107         {
2108             CPLError(CE_Warning, CPLE_AppDefined,
2109                      "Trimmming maximum size of file 2GB from %.1fGB\n"
2110                      "to avoid overflow of kdu_long layer size.",
2111                      dfLayerBytes / 1000000000.0);
2112             dfLayerBytes = 2000000000.0;
2113         }
2114 
2115         layer_bytes[layer_count - 1] = static_cast<kdu_long>(dfLayerBytes);
2116 
2117         CPLDebug("JP2KAK", "layer_bytes[] = %g\n",
2118                  static_cast<double>(layer_bytes[layer_count - 1]));
2119     }
2120     else
2121     {
2122         bReversible = true;
2123     }
2124 
2125     // Do we want to use more than one tile?
2126     int nTileXSize = nXSize;
2127     int nTileYSize = nYSize;
2128 
2129     if( nTileXSize > 25000 )
2130     {
2131         // Don't generate tiles that are terrible wide by default, as
2132         // they consume a lot of memory for the compression engine.
2133         nTileXSize = 20000;
2134     }
2135 
2136     // TODO(schwehr): Why 253 and not 255?
2137     if( (nTileYSize / TILE_CHUNK_SIZE) > 253)
2138     {
2139         // We don't want to process a tile in more than 255 chunks as there
2140         // is a limit on the number of tile parts in a tile and we are likely
2141         // to flush out a tile part for each processing chunk.  If we might
2142         // go over try trimming our Y tile size such that we will get about
2143         // 200 tile parts.
2144         nTileYSize = 200 * TILE_CHUNK_SIZE;
2145     }
2146 
2147     if( CSLFetchNameValue( papszOptions, "BLOCKXSIZE" ) != nullptr )
2148         nTileXSize = atoi(CSLFetchNameValue( papszOptions, "BLOCKXSIZE"));
2149 
2150     if( CSLFetchNameValue( papszOptions, "BLOCKYSIZE" ) != nullptr )
2151         nTileYSize = atoi(CSLFetchNameValue( papszOptions, "BLOCKYSIZE"));
2152     if( nTileXSize <= 0 || nTileYSize <= 0 )
2153     {
2154         CPLError(CE_Failure, CPLE_AppDefined,
2155                  "Wrong value for BLOCKXSIZE/BLOCKYSIZE");
2156         return nullptr;
2157     }
2158 
2159     // Avoid splitting into too many tiles - apparently limiting to 64K tiles.
2160     // There is a hard limit on the number of tiles allowed in JPEG2000.
2161     const double dfXbyY = static_cast<double>(nXSize) * nYSize / (1024 * 64);
2162     while( dfXbyY >= static_cast<double>(nTileXSize) * nTileYSize )
2163     {
2164         nTileXSize *= 2;
2165         nTileYSize *= 2;
2166     }
2167 
2168     if( nTileXSize > nXSize ) nTileXSize = nXSize;
2169     if( nTileYSize > nYSize ) nTileYSize = nYSize;
2170 
2171     CPLDebug("JP2KAK", "Final JPEG2000 Tile Size is %dP x %dL.",
2172              nTileXSize, nTileYSize);
2173 
2174     // Do we want a comment segment emitted?
2175     const bool bComseg = CPLFetchBool(papszOptions, "COMSEG", true);
2176 
2177     // Work out the precision.
2178     int nBits = 0;
2179 
2180     if( CSLFetchNameValue( papszOptions, "NBITS" ) != nullptr )
2181         nBits = atoi(CSLFetchNameValue(papszOptions, "NBITS"));
2182     else if( poPrototypeBand->GetMetadataItem("NBITS", "IMAGE_STRUCTURE")
2183              != nullptr )
2184         nBits =
2185             atoi(poPrototypeBand->GetMetadataItem("NBITS", "IMAGE_STRUCTURE"));
2186     else
2187         nBits = GDALGetDataTypeSize(eType);
2188 
2189     // Establish the general image parameters.
2190     siz_params oSizeParams;
2191 
2192     oSizeParams.set(Scomponents, 0, 0, poSrcDS->GetRasterCount());
2193     oSizeParams.set(Sdims, 0, 0, nYSize);
2194     oSizeParams.set(Sdims, 0, 1, nXSize);
2195     oSizeParams.set(Sprecision, 0, 0, nBits);
2196     if( eType == GDT_UInt16 || eType == GDT_Byte )
2197         oSizeParams.set(Ssigned, 0, 0, false);
2198     else
2199         oSizeParams.set(Ssigned, 0, 0, true);
2200 
2201     if( nTileXSize != nXSize || nTileYSize != nYSize )
2202     {
2203         oSizeParams.set(Stiles, 0, 0, nTileYSize);
2204         oSizeParams.set(Stiles, 0, 1, nTileXSize);
2205 
2206         CPLDebug("JP2KAK", "Stiles=%d,%d", nTileYSize, nTileXSize);
2207     }
2208 
2209     kdu_params *poSizeRef = &oSizeParams;
2210     poSizeRef->finalize();
2211 
2212     // Open output file, and setup codestream.
2213     if( !pfnProgress( 0.0, nullptr, pProgressData ) )
2214         return nullptr;
2215 
2216     jp2_family_tgt family;
2217 #ifdef KAKADU_JPX
2218     jpx_family_tgt jpx_family;
2219     jpx_target jpx_out;
2220     const bool bIsJPX = !EQUAL(CPLGetExtension(pszFilename), "jpf") &&
2221                         !EQUAL(CPLGetExtension(pszFilename), "jpc") &&
2222                         !EQUAL(CPLGetExtension(pszFilename), "j2k") &&
2223                         !(pszCodec != NULL && EQUAL(pszCodec, "J2K"));
2224 #endif
2225 
2226     kdu_compressed_target *poOutputFile = nullptr;
2227     jp2_target jp2_out;
2228     const char* pszCodec = CSLFetchNameValueDef(papszOptions, "CODEC", nullptr);
2229     const bool bIsJP2 = (!EQUAL(CPLGetExtension(pszFilename), "jpc") &&
2230                          !EQUAL(CPLGetExtension(pszFilename), "j2k") &&
2231 #ifdef KAKADU_JPX
2232                          !bIsJPX &&
2233 #endif
2234                          !(pszCodec != nullptr && EQUAL(pszCodec, "J2K")) ) ||
2235                         (pszCodec != nullptr && EQUAL(pszCodec, "JP2"));
2236     kdu_codestream oCodeStream;
2237 
2238     vsil_target oVSILTarget;
2239 
2240     try
2241     {
2242         oVSILTarget.open(pszFilename, "w");
2243 
2244         if( bIsJP2 )
2245         {
2246             // family.open( pszFilename );
2247             family.open(&oVSILTarget);
2248 
2249             jp2_out.open(&family);
2250             poOutputFile = &jp2_out;
2251         }
2252 #ifdef KAKADU_JPX
2253         else if( bIsJPX )
2254         {
2255             jpx_family.open(pszFilename);
2256 
2257             jpx_out.open(&jpx_family);
2258             jpx_out.add_codestream();
2259         }
2260 #endif
2261         else
2262         {
2263             poOutputFile = &oVSILTarget;
2264         }
2265 
2266         oCodeStream.create(&oSizeParams, poOutputFile);
2267     }
2268     catch( ... )
2269     {
2270         return nullptr;
2271     }
2272 
2273     // Do we have a high res region of interest?
2274     kdu_roi_image *poROIImage = nullptr;
2275     const char* pszROI = CSLFetchNameValue(papszOptions, "ROI");
2276     if( pszROI )
2277     {
2278         CPLStringList aosTokens(
2279             CSLTokenizeStringComplex(pszROI, ",", FALSE, FALSE));
2280 
2281         if( aosTokens.size() != 4 )
2282         {
2283             CPLError(CE_Warning, CPLE_AppDefined,
2284                      "Skipping corrupt ROI def = \n%s", pszROI);
2285         }
2286         else
2287         {
2288             kdu_dims region;
2289             region.pos.x = atoi(aosTokens[0]);
2290             region.pos.y = atoi(aosTokens[1]);
2291             region.size.x = atoi(aosTokens[2]);
2292             region.size.y = atoi(aosTokens[3]);
2293 
2294             poROIImage = new kdu_roi_rect(oCodeStream, region);
2295         }
2296     }
2297 
2298     // Set some particular parameters.
2299     oCodeStream.access_siz()->parse_string(
2300         CPLString().Printf("Clayers=%d", layer_count).c_str());
2301     oCodeStream.access_siz()->parse_string("Cycc=no");
2302     if( eType == GDT_Int16 || eType == GDT_UInt16 )
2303         oCodeStream.access_siz()->parse_string("Qstep=0.0000152588");
2304 
2305     if( bReversible )
2306         oCodeStream.access_siz()->parse_string("Creversible=yes");
2307     else
2308         oCodeStream.access_siz()->parse_string("Creversible=no");
2309 
2310     // Set some user-overridable parameters.
2311     const char * const apszParams[] =
2312         { "Corder", "PCRL",
2313           "Cprecincts",
2314           "{512,512},{256,512},{128,512},{64,512},{32,512},{16,512},{8,512},{4,512},{2,512}",
2315           "ORGgen_plt", "yes",
2316           "ORGgen_tlm", nullptr,
2317           "ORGtparts", nullptr,
2318           "Qguard", nullptr,
2319           "Cmodes", nullptr,
2320           "Clevels", nullptr,
2321           "Cblk", nullptr,
2322           "Rshift", nullptr,
2323           "Rlevels", nullptr,
2324           "Rweight", nullptr,
2325           "Sprofile", nullptr,
2326           nullptr, nullptr };
2327 
2328     for( int iParam = 0; apszParams[iParam] != nullptr; iParam += 2 )
2329     {
2330         const char *pszValue =
2331             CSLFetchNameValue(papszOptions, apszParams[iParam]);
2332 
2333         if( pszValue == nullptr )
2334             pszValue = apszParams[iParam + 1];
2335 
2336         if( pszValue != nullptr )
2337         {
2338             CPLString osOpt;
2339 
2340             osOpt.Printf("%s=%s", apszParams[iParam], pszValue);
2341             try
2342             {
2343                 oCodeStream.access_siz()->parse_string(osOpt);
2344             }
2345             catch( ... )
2346             {
2347                 if( bIsJP2 )
2348                 {
2349                     jp2_out.close();
2350                     family.close();
2351                 }
2352                 else
2353                 {
2354                     poOutputFile->close();
2355                 }
2356                 return nullptr;
2357             }
2358 
2359             CPLDebug("JP2KAK", "parse_string(%s)", osOpt.c_str());
2360         }
2361     }
2362 
2363     oCodeStream.access_siz()->finalize_all();
2364 
2365     // Some JP2 specific parameters.
2366     if( bIsJP2 )
2367     {
2368         // Set dimensional information.
2369         // All redundant with the SIZ marker segment.
2370         jp2_dimensions dims = jp2_out.access_dimensions();
2371         dims.init(&oSizeParams);
2372 
2373         // Set colour space information (mandatory).
2374         jp2_colour colour = jp2_out.access_colour();
2375 
2376         if( bHaveCT || poSrcDS->GetRasterCount() == 3 )
2377         {
2378             colour.init( JP2_sRGB_SPACE );
2379         }
2380         else if( poSrcDS->GetRasterCount() >= 4
2381                  && poSrcDS->GetRasterBand(4)->GetColorInterpretation()
2382                  == GCI_AlphaBand )
2383         {
2384             colour.init(JP2_sRGB_SPACE);
2385             jp2_out.access_channels().init(3);
2386             jp2_out.access_channels().set_colour_mapping(0, 0);
2387             jp2_out.access_channels().set_colour_mapping(1, 1);
2388             jp2_out.access_channels().set_colour_mapping(2, 2);
2389             jp2_out.access_channels().set_opacity_mapping(0, 3);
2390             jp2_out.access_channels().set_opacity_mapping(1, 3);
2391             jp2_out.access_channels().set_opacity_mapping(2, 3);
2392         }
2393         else if( poSrcDS->GetRasterCount() >= 2
2394                  && poSrcDS->GetRasterBand(2)->GetColorInterpretation()
2395                  == GCI_AlphaBand )
2396         {
2397             colour.init(JP2_sLUM_SPACE);
2398             jp2_out.access_channels().init(1);
2399             jp2_out.access_channels().set_colour_mapping(0, 0);
2400             jp2_out.access_channels().set_opacity_mapping(0, 1);
2401         }
2402         else
2403         {
2404             colour.init( JP2_sLUM_SPACE );
2405         }
2406 
2407         // Resolution.
2408         if( poSrcDS->GetMetadataItem("TIFFTAG_XRESOLUTION") != nullptr
2409             && poSrcDS->GetMetadataItem("TIFFTAG_YRESOLUTION") != nullptr
2410             && poSrcDS->GetMetadataItem("TIFFTAG_RESOLUTIONUNIT") != nullptr )
2411         {
2412             jp2_resolution res = jp2_out.access_resolution();
2413             double dfXRes =
2414                 CPLAtof(poSrcDS->GetMetadataItem("TIFFTAG_XRESOLUTION"));
2415             double dfYRes =
2416                 CPLAtof(poSrcDS->GetMetadataItem("TIFFTAG_YRESOLUTION"));
2417 
2418             if( atoi(poSrcDS->GetMetadataItem("TIFFTAG_RESOLUTIONUNIT")) == 2 )
2419             {
2420                 // Convert pixels per inch to pixels per cm.
2421                 // TODO(schwehr): Change this to 1.0 / 2.54 if correct.
2422                 const double dfInchToCm = 39.37 / 100.0;
2423                 dfXRes *= dfInchToCm;
2424                 dfYRes *= dfInchToCm;
2425             }
2426 
2427             // Convert to pixels per meter.
2428             dfXRes *= 100.0;
2429             dfYRes *= 100.0;
2430 
2431             if( dfXRes != 0.0 && dfYRes != 0.0 )
2432             {
2433                 if( fabs(dfXRes / dfYRes - 1.0) > 0.00001 )
2434                     res.init(static_cast<float>(dfYRes / dfXRes));
2435                 else
2436                     res.init(1.0);
2437                 res.set_resolution(static_cast<float>(dfXRes), true);
2438             }
2439         }
2440     }
2441 
2442     // Write JP2 pseudocolor table if available.
2443     if( bIsJP2 && bHaveCT )
2444     {
2445         GDALColorTable *const poCT = poPrototypeBand->GetColorTable();
2446         const int nCount = poCT->GetColorEntryCount();
2447         kdu_int32 *panLUT =
2448             static_cast<kdu_int32 *>(CPLMalloc(sizeof(kdu_int32) * nCount * 3));
2449 
2450         jp2_palette oJP2Palette = jp2_out.access_palette();
2451         oJP2Palette.init(3, nCount);
2452 
2453         for( int iColor = 0; iColor < nCount; iColor++ )
2454         {
2455             GDALColorEntry sEntry = { 0, 0, 0, 0 };
2456 
2457             poCT->GetColorEntryAsRGB(iColor, &sEntry);
2458             panLUT[iColor + nCount * 0] = sEntry.c1;
2459             panLUT[iColor + nCount * 1] = sEntry.c2;
2460             panLUT[iColor + nCount * 2] = sEntry.c3;
2461         }
2462 
2463         oJP2Palette.set_lut(0, panLUT + nCount * 0, 8, false);
2464         oJP2Palette.set_lut(1, panLUT + nCount * 1, 8, false);
2465         oJP2Palette.set_lut(2, panLUT + nCount * 2, 8, false);
2466 
2467         CPLFree(panLUT);
2468 
2469         jp2_channels oJP2Channels = jp2_out.access_channels();
2470 
2471         oJP2Channels.init(3);
2472         oJP2Channels.set_colour_mapping(0, 0, 0);
2473         oJP2Channels.set_colour_mapping(1, 0, 1);
2474         oJP2Channels.set_colour_mapping(2, 0, 2);
2475     }
2476 
2477     if( bIsJP2 )
2478     {
2479         try
2480         {
2481             jp2_out.write_header();
2482         }
2483         catch( ... )
2484         {
2485             CPLDebug("JP2KAK", "jp2_out.write_header() - caught exception.");
2486             oCodeStream.destroy();
2487             return nullptr;
2488         }
2489     }
2490 
2491     // Set the GeoTIFF and GML boxes if georeferencing is available,
2492     // and this is a JP2 file.
2493     double adfGeoTransform[6] = { 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 };
2494     if( bIsJP2
2495         && ((poSrcDS->GetGeoTransform(adfGeoTransform) == CE_None
2496              && (adfGeoTransform[0] != 0.0
2497                  || adfGeoTransform[1] != 1.0
2498                  || adfGeoTransform[2] != 0.0
2499                  || adfGeoTransform[3] != 0.0
2500                  || adfGeoTransform[4] != 0.0
2501                  || std::abs(adfGeoTransform[5]) != 1.0))
2502             || poSrcDS->GetGCPCount() > 0
2503             || poSrcDS->GetMetadata("RPC") != nullptr) )
2504     {
2505         GDALJP2Metadata oJP2MD;
2506 
2507         if( poSrcDS->GetGCPCount() > 0 )
2508         {
2509             oJP2MD.SetProjection(poSrcDS->GetGCPProjection());
2510             oJP2MD.SetGCPs(poSrcDS->GetGCPCount(), poSrcDS->GetGCPs());
2511         }
2512         else
2513         {
2514             oJP2MD.SetProjection(poSrcDS->GetProjectionRef());
2515             oJP2MD.SetGeoTransform(adfGeoTransform);
2516         }
2517 
2518         oJP2MD.SetRPCMD(poSrcDS->GetMetadata("RPC"));
2519 
2520         const char *const pszAreaOrPoint =
2521             poSrcDS->GetMetadataItem(GDALMD_AREA_OR_POINT);
2522         oJP2MD.bPixelIsPoint =
2523             pszAreaOrPoint != nullptr && EQUAL(pszAreaOrPoint, GDALMD_AOP_POINT);
2524 
2525         if( CPLFetchBool(papszOptions, "GMLJP2", true) )
2526         {
2527             const char *pszGMLJP2V2Def =
2528                 CSLFetchNameValue(papszOptions, "GMLJP2V2_DEF");
2529             GDALJP2Box* poBox;
2530             if( pszGMLJP2V2Def != nullptr ) {
2531                 poBox = oJP2MD.CreateGMLJP2V2(
2532                     nXSize,nYSize,pszGMLJP2V2Def,poSrcDS);
2533             } else {
2534                 poBox = oJP2MD.CreateGMLJP2(nXSize, nYSize);
2535             }
2536             try
2537             {
2538                 JP2KAKWriteBox(&family, poBox);
2539             }
2540             catch( ... )
2541             {
2542                 CPLDebug("JP2KAK", "JP2KAKWriteBox) - caught exception.");
2543                 oCodeStream.destroy();
2544                 delete poBox;
2545                 return nullptr;
2546             }
2547         }
2548         if( CPLFetchBool(papszOptions, "GeoJP2", true) ) {
2549             GDALJP2Box* poBox = oJP2MD.CreateJP2GeoTIFF();
2550             try
2551             {
2552                 JP2KAKWriteBox(&family, poBox);
2553             }
2554             catch( ... )
2555             {
2556                 CPLDebug("JP2KAK", "JP2KAKWriteBox) - caught exception.");
2557                 oCodeStream.destroy();
2558                 delete poBox;
2559                 return nullptr;
2560             }
2561         }
2562     }
2563 
2564     // Do we have any XML boxes we want to preserve?
2565 
2566     for( int iBox = 0; true; iBox++ )
2567     {
2568         CPLString oName;
2569         oName.Printf("xml:BOX_%d", iBox);
2570         char **papszMD = poSrcDS->GetMetadata(oName);
2571 
2572         if( papszMD == nullptr || CSLCount(papszMD) != 1 )
2573             break;
2574 
2575         GDALJP2Box *poXMLBox = new GDALJP2Box();
2576 
2577         poXMLBox->SetType("xml ");
2578         poXMLBox->SetWritableData(static_cast<int>(strlen(papszMD[0]) + 1),
2579                                   reinterpret_cast<GByte *>(papszMD[0]));
2580         JP2KAKWriteBox(&family, poXMLBox);
2581     }
2582 
2583     // Open codestream box.
2584     if( bIsJP2 )
2585         jp2_out.open_codestream();
2586 
2587     // Create one big tile, and a compressing engine, and line
2588     // buffer for each component.
2589     double dfPixelsDone = 0.0;
2590     const double dfPixelsTotal = nXSize * static_cast<double>(nYSize);
2591     const bool bFlushEnabled = CPLFetchBool(papszOptions, "FLUSH", true);
2592 
2593     for( int iTileYOff = 0; iTileYOff < nYSize; iTileYOff += nTileYSize )
2594     {
2595         for( int iTileXOff = 0; iTileXOff < nXSize; iTileXOff += nTileXSize )
2596         {
2597             kdu_tile oTile = oCodeStream.open_tile(
2598                 kdu_coords(iTileXOff / nTileXSize, iTileYOff / nTileYSize));
2599 
2600             // Is this a partial tile on the right or bottom?
2601             const int nThisTileXSize =
2602                 iTileXOff + nTileXSize < nXSize
2603                 ? nTileXSize
2604                 : nXSize - iTileXOff;
2605 
2606             const int nThisTileYSize =
2607                 iTileYOff + nTileYSize < nYSize
2608                 ? nTileYSize
2609                 : nYSize - iTileYOff;
2610 
2611             // Setup scaled progress monitor.
2612 
2613             const double dfPixelsDoneAfter =
2614                 dfPixelsDone + static_cast<double>(nThisTileXSize) * nThisTileYSize;
2615 
2616             void *pScaledProgressData = GDALCreateScaledProgress(
2617                 dfPixelsDone / dfPixelsTotal, dfPixelsDoneAfter / dfPixelsTotal,
2618                 pfnProgress, pProgressData);
2619 
2620             if( !JP2KAKCreateCopy_WriteTile(poSrcDS, oTile, poROIImage,
2621                                             iTileXOff, iTileYOff,
2622                                             nThisTileXSize, nThisTileYSize,
2623                                             bReversible, nBits, eType,
2624                                             oCodeStream, bFlushEnabled,
2625                                             layer_bytes.data(), layer_count,
2626                                             GDALScaledProgress,
2627                                             pScaledProgressData, bComseg) )
2628             {
2629                 GDALDestroyScaledProgress(pScaledProgressData);
2630 
2631                 oCodeStream.destroy();
2632                 poOutputFile->close();
2633                 VSIUnlink(pszFilename);
2634                 return nullptr;
2635             }
2636 
2637             GDALDestroyScaledProgress(pScaledProgressData);
2638             dfPixelsDone = dfPixelsDoneAfter;
2639 
2640             oTile.close();
2641         }
2642     }
2643 
2644     // Finish flushing out results.
2645     oCodeStream.flush(layer_bytes.data(), layer_count, nullptr, true, bComseg);
2646     oCodeStream.destroy();
2647 
2648     if( bIsJP2 )
2649     {
2650         jp2_out.close();
2651         family.close();
2652     }
2653     else
2654     {
2655         poOutputFile->close();
2656     }
2657 
2658     oVSILTarget.close();
2659 
2660     if( !pfnProgress(1.0, nullptr, pProgressData) )
2661         return nullptr;
2662 
2663     // Re-open dataset, and copy any auxiliary pam information.
2664     GDALOpenInfo oOpenInfo(pszFilename, GA_ReadOnly);
2665     GDALPamDataset *poDS =
2666         static_cast<GDALPamDataset *>(JP2KAKDataset::Open(&oOpenInfo));
2667 
2668     if( poDS )
2669         poDS->CloneInfo(poSrcDS, GCIF_PAM_DEFAULT);
2670 
2671     return poDS;
2672 }
2673 
2674 /************************************************************************/
2675 /*                        GDALRegister_JP2KAK()                         */
2676 /************************************************************************/
2677 
2678 void GDALRegister_JP2KAK()
2679 
2680 {
2681     if( !GDAL_CHECK_VERSION("JP2KAK driver") )
2682         return;
2683 
2684     if( GDALGetDriverByName("JP2KAK") != nullptr )
2685         return;
2686 
2687     GDALDriver *poDriver = new GDALDriver();
2688 
2689     poDriver->SetDescription("JP2KAK");
2690     poDriver->SetMetadataItem(GDAL_DCAP_RASTER, "YES");
2691     poDriver->SetMetadataItem(GDAL_DCAP_VECTOR, "YES");
2692     poDriver->SetMetadataItem(
2693         GDAL_DMD_LONGNAME, "JPEG-2000 (based on Kakadu " KDU_CORE_VERSION ")");
2694     poDriver->SetMetadataItem(GDAL_DMD_HELPTOPIC, "drivers/raster/jp2kak.html");
2695     poDriver->SetMetadataItem(GDAL_DMD_CREATIONDATATYPES, "Byte Int16 UInt16");
2696     poDriver->SetMetadataItem(GDAL_DMD_MIMETYPE, "image/jp2");
2697     poDriver->SetMetadataItem(GDAL_DMD_EXTENSION, "jp2 j2k");
2698     poDriver->SetMetadataItem(GDAL_DCAP_VIRTUALIO, "YES");
2699 
2700     poDriver->SetMetadataItem(GDAL_DMD_OPENOPTIONLIST,
2701 "<OpenOptionList>"
2702 "   <Option name='1BIT_ALPHA_PROMOTION' type='boolean' description="
2703 "'Whether a 1-bit alpha channel should be promoted to 8-bit' default='YES'/>"
2704 "   <Option name='OPEN_REMOTE_GML' type='boolean' description="
2705 "'Whether to load remote vector layers referenced by "
2706 "a link in a GMLJP2 v2 box' default='NO'/>"
2707 "   <Option name='GEOREF_SOURCES' type='string' description="
2708 "'Comma separated list made with values "
2709 "INTERNAL/GMLJP2/GEOJP2/WORLDFILE/PAM/NONE that describe the priority order "
2710 "for georeferencing' default='PAM,GEOJP2,GMLJP2,WORLDFILE'/>"
2711 "</OpenOptionList>" );
2712 
2713     poDriver->SetMetadataItem(GDAL_DMD_CREATIONOPTIONLIST,
2714 "<CreationOptionList>"
2715 "   <Option name='CODEC' type='string-select' "
2716     "default='according to file extension. If unknown, default to JP2'>"
2717 "       <Value>JP2</Value>"
2718 "       <Value>J2K</Value>"
2719 "   </Option>"
2720 "   <Option name='QUALITY' type='integer' description="
2721 "'0.01-100, 100 is lossless'/>"
2722 "   <Option name='BLOCKXSIZE' type='int' description='Tile Width'/>"
2723 "   <Option name='BLOCKYSIZE' type='int' description='Tile Height'/>"
2724 "   <Option name='GeoJP2' type='boolean' description='defaults to ON'/>"
2725 "   <Option name='GMLJP2' type='boolean' description='defaults to ON'/>"
2726 "   <Option name='GMLJP2V2_DEF' type='string' description="
2727 "'Definition file to describe how a GMLJP2 v2 box should be generated. "
2728 "If set to YES, a minimal instance will be created'/>"
2729 "   <Option name='LAYERS' type='integer'/>"
2730 "   <Option name='ROI' type='string'/>"
2731 "   <Option name='COMSEG' type='boolean' />"
2732 "   <Option name='FLUSH' type='boolean' />"
2733 "   <Option name='NBITS' type='int' description="
2734 "'BITS (precision) for sub-byte files (1-7), sub-uint16 (9-15)'/>"
2735 "   <Option name='Corder' type='string'/>"
2736 "   <Option name='Cprecincts' type='string'/>"
2737 "   <Option name='Cmodes' type='string'/>"
2738 "   <Option name='Clevels' type='string'/>"
2739 "   <Option name='ORGgen_plt' type='string'/>"
2740 "   <Option name='ORGgen_tlm' type='string'/>"
2741 "   <Option name='ORGtparts' type='string'/>"
2742 "   <Option name='Qguard' type='integer'/>"
2743 "   <Option name='Sprofile' type='string'/>"
2744 "   <Option name='Rshift' type='string'/>"
2745 "   <Option name='Rlevels' type='string'/>"
2746 "   <Option name='Rweight' type='string'/>"
2747 "</CreationOptionList>" );
2748 
2749     poDriver->pfnOpen = JP2KAKDataset::Open;
2750     poDriver->pfnIdentify = JP2KAKDataset::Identify;
2751     poDriver->pfnCreateCopy = JP2KAKCreateCopy;
2752 
2753     GetGDALDriverManager()->RegisterDriver(poDriver);
2754 }
2755