otbGDALImageIO.cxx 69 KB
Newer Older
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19
/*
 * Copyright (C) 2005-2017 Centre National d'Etudes Spatiales (CNES)
 *
 * This file is part of Orfeo Toolbox
 *
 *     https://www.orfeo-toolbox.org/
 *
 * Licensed under the Apache License, Version 2.0 (the "License");
 * you may not use this file except in compliance with the License.
 * You may obtain a copy of the License at
 *
 *     http://www.apache.org/licenses/LICENSE-2.0
 *
 * Unless required by applicable law or agreed to in writing, software
 * distributed under the License is distributed on an "AS IS" BASIS,
 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 * See the License for the specific language governing permissions and
 * limitations under the License.
 */
Patrick Imbo's avatar
nomsg  
Patrick Imbo committed
20

21 22
#include <iostream>
#include <fstream>
23
#include <vector>
Thomas Feuvrier's avatar
nomsg  
Thomas Feuvrier committed
24 25

#include "otbGDALImageIO.h"
Thomas Feuvrier's avatar
nomsg  
Thomas Feuvrier committed
26
#include "otbMacro.h"
27
#include "otbSystem.h"
28
#include "otbStopwatch.h"
29
#include "itksys/SystemTools.hxx"
30
#include "otbImage.h"
31
#include "otb_tinyxml.h"
32
#include "otbImageKeywordlist.h"
33

Caroline Ruffel's avatar
nomsg  
Caroline Ruffel committed
34
#include "itkMetaDataObject.h"
35
#include "otbMetaDataKey.h"
36

37 38 39
#include "itkRGBPixel.h"
#include "itkRGBAPixel.h"

40 41 42 43
#include "cpl_conv.h"
#include "ogr_spatialref.h"
#include "ogr_srs_api.h"

44 45
#include "otbGDALDriverManagerWrapper.h"

46 47
#include "otb_boost_string_header.h"

48 49
#include "otbOGRHelpers.h"

50 51
#include "stdint.h" //needed for uintptr_t

52 53 54 55
inline unsigned int uint_ceildivpow2(unsigned int a, unsigned int b) {
  return (a + (1 << b) - 1) >> b;
}

Thomas Feuvrier's avatar
nomsg  
Thomas Feuvrier committed
56 57
namespace otb
{
Julien Malik's avatar
Julien Malik committed
58

59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75
class GDALDataTypeWrapper
{
public:
  GDALDataTypeWrapper() : pixType(GDT_Byte) {}
  ~GDALDataTypeWrapper() {}
  GDALDataTypeWrapper(const GDALDataTypeWrapper& w)
  {
    pixType = w.pixType;
  }
  GDALDataTypeWrapper& operator=(GDALDataTypeWrapper w)
  {
    pixType = w.pixType;
    return *this;
  }
  GDALDataType pixType;
}; // end of GDALDataTypeWrapper

76

Julien Malik's avatar
Julien Malik committed
77
/*
78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125
template<class InputType>
void printOutputData(InputType *pData, int nbBands, int nbPixelToRead)
{
  for (unsigned int itPxl = 0; itPxl < (unsigned int) (nbPixelToRead * nbBands); itPxl++)
    {
    std::cout << "Buffer["<< itPxl << "] = " << *(pData + itPxl) << std::endl;
    }
};

void printDataBuffer(unsigned char *pData, GDALDataType pxlType, int nbBands, int nbPixelToRead)
{
  if (pxlType == GDT_Int16)
    {
    printOutputData( static_cast<short*>( static_cast<void*>(pData) ), nbBands, nbPixelToRead);
    }
  else if (pxlType == GDT_Int32)
    {
    printOutputData( static_cast<int*>( static_cast<void*>(pData) ), nbBands, nbPixelToRead);
    }
  else if (pxlType == GDT_Float32)
    {
    printOutputData( static_cast<float*>( static_cast<void*>(pData) ), nbBands, nbPixelToRead);
    }
  else if (pxlType == GDT_Float64)
    {
    printOutputData( static_cast<double*>( static_cast<void*>(pData) ), nbBands, nbPixelToRead);
    }
  else if (pxlType == GDT_CInt16)
    {
    printOutputData( static_cast<std::complex<short>*>( static_cast<void*>(pData) ), nbBands, nbPixelToRead);
    }
  else if (pxlType == GDT_CInt32)
    {
    printOutputData( static_cast<std::complex<int>*>( static_cast<void*>(pData) ), nbBands, nbPixelToRead);
    }
  else if (pxlType == GDT_CFloat32)
    {
    printOutputData( static_cast<std::complex<float>*>( static_cast<void*>(pData) ), nbBands, nbPixelToRead);
    }
  else if (pxlType == GDT_CFloat64)
    {
    printOutputData( static_cast<std::complex<double>*>( static_cast<void*>(pData) ), nbBands, nbPixelToRead);
    }
  else
    {
    std::cerr << "Pixel type unknown" << std::endl;
    }
};
Julien Malik's avatar
Julien Malik committed
126
*/
Thomas Feuvrier's avatar
nomsg  
Thomas Feuvrier committed
127

128 129
GDALImageIO::GDALImageIO()
{
Thomas Feuvrier's avatar
nomsg  
Thomas Feuvrier committed
130
  // By default set number of dimensions to two.
131
  this->SetNumberOfDimensions(2);
Patrick Imbo's avatar
Patrick Imbo committed
132

Thomas Feuvrier's avatar
nomsg  
Thomas Feuvrier committed
133
  // By default set pixel type to scalar.
134
  m_PixelType = SCALAR;
Patrick Imbo's avatar
Patrick Imbo committed
135

Thomas Feuvrier's avatar
nomsg  
Thomas Feuvrier committed
136
  // By default set component type to unsigned char
137 138 139
  m_ComponentType = UCHAR;
  m_UseCompression = false;
  m_CompressionLevel = 4; // Range 0-9; 0 = no file compression, 9 = maximum file compression
Thomas Feuvrier's avatar
nomsg  
Thomas Feuvrier committed
140 141

  // Set default spacing to one
142 143
  m_Spacing[0] = 1.0;
  m_Spacing[1] = 1.0;
144 145 146
  // Set default origin to half a pixel (centered pixel convention)
  m_Origin[0] = 0.5;
  m_Origin[1] = 0.5;
Thomas Feuvrier's avatar
nomsg  
Thomas Feuvrier committed
147

148
  m_IsIndexed   = false;
149
  m_DatasetNumber = 0;
150 151 152
  //m_poBands     = NULL;
  //m_hDriver     = NULL;
  //m_poDataset   = NULL;
153

154 155
  m_NbBands = 0;
  m_FlagWriteImageInformation = true;
Thomas Feuvrier's avatar
nomsg  
Thomas Feuvrier committed
156

157
  m_CanStreamWrite = false;
158
  m_IsComplex = false;
Julien Malik's avatar
Julien Malik committed
159
  m_IsVectorImage = false;
160 161

  m_PxType = new GDALDataTypeWrapper;
162 163

  m_NumberOfOverviews = 0;
164
  m_ResolutionFactor = 0;
165
  m_BytePerPixel = 0;
166
  m_WriteRPCTags = false;
167
}
Thomas Feuvrier's avatar
nomsg  
Thomas Feuvrier committed
168

169 170
GDALImageIO::~GDALImageIO()
{
171
  delete m_PxType;
172
}
Thomas Feuvrier's avatar
nomsg  
Thomas Feuvrier committed
173 174

// Tell only if the file can be read with GDAL.
175 176
bool GDALImageIO::CanReadFile(const char* file)
{
Thomas Feuvrier's avatar
nomsg  
Thomas Feuvrier committed
177
  // First check the extension
178
  if (file == ITK_NULLPTR)
OTB Bot's avatar
STYLE  
OTB Bot committed
179 180
    {
    itkDebugMacro(<< "No filename specified.");
181
    return false;
OTB Bot's avatar
STYLE  
OTB Bot committed
182
    }
183 184
  m_Dataset = GDALDriverManagerWrapper::GetInstance().Open(file);
  return m_Dataset.IsNotNull();
185 186 187 188 189 190 191
}

// Used to print information about this object
void GDALImageIO::PrintSelf(std::ostream& os, itk::Indent indent) const
{
  Superclass::PrintSelf(os, indent);
  os << indent << "Compression Level : " << m_CompressionLevel << "\n";
192
  os << indent << "IsComplex (otb side) : " << m_IsComplex << "\n";
193
  os << indent << "Byte per pixel : " << m_BytePerPixel << "\n";
194
}
Thomas Feuvrier's avatar
nomsg  
Thomas Feuvrier committed
195 196

// Read a 3D image (or event more bands)... not implemented yet
197 198 199
void GDALImageIO::ReadVolume(void*)
{
}
Thomas Feuvrier's avatar
nomsg  
Thomas Feuvrier committed
200 201

// Read image with GDAL
202 203
void GDALImageIO::Read(void* buffer)
{
204 205
  // Convert buffer from void * to unsigned char *
  unsigned char *p = static_cast<unsigned char *>(buffer);
206 207

  // Check if conversion succeed
208
  if (p == ITK_NULLPTR)
OTB Bot's avatar
STYLE  
OTB Bot committed
209
    {
210
    itkExceptionMacro(<< "GDAL : Bad alloc");
211
    return;
OTB Bot's avatar
STYLE  
OTB Bot committed
212
    }
Guillaume Borrut's avatar
Guillaume Borrut committed
213

214 215 216 217
  // Get the origin of the region to read
  int lFirstLineRegion   = this->GetIORegion().GetIndex()[1];
  int lFirstColumnRegion = this->GetIORegion().GetIndex()[0];

218
  // Get nb. of lines and columns of the region to read
219 220 221 222 223 224 225
  int lNbLinesRegion   = this->GetIORegion().GetSize()[1];
  int lNbColumnsRegion = this->GetIORegion().GetSize()[0];

  //std::cout << "OriginBuffer= " <<  lFirstLineRegion << " x " << lFirstColumnRegion << std::endl;
  //std::cout << "SizeBuffer= " <<  lNbLinesRegion << " x " << lNbColumnsRegion << std::endl;

  // Compute the origin of the image region to read at the initial resolution
226 227
  int lFirstLine   = lFirstLineRegion * (1 << m_ResolutionFactor);
  int lFirstColumn = lFirstColumnRegion * (1 << m_ResolutionFactor);
228 229 230 231

  //std::cout << "OriginImage= " <<  lFirstLine << " x " << lFirstColumn << std::endl;

  // Compute the size of the image region to read at the initial resolution
232 233
  int lNbLines     = lNbLinesRegion * (1 << m_ResolutionFactor);
  int lNbColumns   = lNbColumnsRegion * (1 << m_ResolutionFactor);
234 235

  // Check if the image region is correct
236 237 238 239
  if (lFirstLine + lNbLines > static_cast<int>(m_OriginalDimensions[1]))
    lNbLines = static_cast<int>(m_OriginalDimensions[1]-lFirstLine);
  if (lFirstColumn + lNbColumns > static_cast<int>(m_OriginalDimensions[0]))
    lNbColumns = static_cast<int>(m_OriginalDimensions[0]-lFirstColumn);
240 241 242

  //std::cout << "SizeImage= " <<  lNbLines << " x " << lNbColumns << std::endl;

243

244
  GDALDataset* dataset = m_Dataset->GetDataSet();
245

246
  // This special case is due to the fact the CINT/CLONG types
247
  // do not exists in ITK. In this case we only report the first band
248
  // TODO This should be fixed
249 250 251
  /*if (GDALDataTypeIsComplex(m_PxType->pixType)
      && (m_PxType->pixType != GDT_CFloat32)
      && (m_PxType->pixType != GDT_CFloat64))
OTB Bot's avatar
STYLE  
OTB Bot committed
252
    {
253 254 255 256 257 258 259 260 261 262
    int pixelOffset = m_BytePerPixel * m_NbBands;
    int lineOffset  = m_BytePerPixel * m_NbBands * lNbColumns;
    int bandOffset  = m_BytePerPixel;
    int nbBands     = m_NbBands;

    int nbPixelToRead = lNbColumns *  lNbLines;
    std::streamoff nbBytes = static_cast<std::streamoff>(m_NbBands) * static_cast<std::streamoff>(nbPixelToRead) * static_cast<std::streamoff>(m_BytePerPixel);
    unsigned char *pBufferTemp = new unsigned char[static_cast<unsigned int>(nbBytes)];

    // keep it for the moment
263
    otbMsgDevMacro(<< "Parameters RasterIO (case CInt and CShort):"
Julien Malik's avatar
Julien Malik committed
264 265 266 267
                   << "\n indX = " << lFirstColumn
                   << "\n indY = " << lFirstLine
                   << "\n sizeX = " << lNbColumns
                   << "\n sizeY = " << lNbLines
268
                   << "\n GDAL Data Type = " << GDALGetDataTypeName(m_PxType->pixType)
Julien Malik's avatar
Julien Malik committed
269 270 271
                   << "\n pixelOffset = " << pixelOffset
                   << "\n lineOffset = " << lineOffset
                   << "\n bandOffset = " << bandOffset);
272 273 274 275 276 277 278 279 280

    CPLErr lCrGdal = m_Dataset->GetDataSet()->RasterIO(GF_Read,
                                                       lFirstColumn,
                                                       lFirstLine,
                                                       lNbColumns,
                                                       lNbLines,
                                                       pBufferTemp, //p, // pData
                                                       lNbColumns,
                                                       lNbLines,
281
                                                       m_PxType->pixType,
282 283 284 285 286 287
                                                       nbBands,
                                                       // We want to read all bands
                                                       NULL,
                                                       pixelOffset,
                                                       lineOffset,
                                                       bandOffset);
288
    // Check for gdal error
289
    if (lCrGdal == CE_Failure)
OTB Bot's avatar
STYLE  
OTB Bot committed
290
      {
291
      itkExceptionMacro(<< "Error while reading image (GDAL format) " << m_FileName );
292
      delete[] pBufferTemp;
293
      return;
OTB Bot's avatar
STYLE  
OTB Bot committed
294
      }
295
    //std::cout << "RAW BUFFER:" <<std::endl;
296
    //printDataBuffer(pBufferTemp, m_PxType->pixType, m_NbBands, lNbColumns*lNbLines);
297

298
    // Convert the buffer to GDT_Float64 type
299
    typedef std::complex<float>           RealType;
300 301
    typedef double                         ScalarRealType;

302
    if (m_PxType->pixType == GDT_CInt32)
303
      {
304
      //std::cout << "Convert input File from GDT_CInt32 to GDT_CFloat32" << std::endl;
305 306 307 308
      typedef std::complex<int>              ComplexIntType;

      for (unsigned int itPxl = 0; itPxl < (unsigned int) (nbPixelToRead * m_NbBands); itPxl++)
        {
309
        ComplexIntType pxlValue = *(static_cast<ComplexIntType*>( static_cast<void*>(pBufferTemp)) + itPxl );
310

311
        RealType    pxlValueReal( static_cast<ScalarRealType>(pxlValue.real()), static_cast<ScalarRealType>(pxlValue.imag()) );
312

313
        memcpy((void*) (&(p[itPxl*sizeof(RealType)])), (const void*) (&(pxlValueReal)), (size_t) (sizeof(RealType)));
314 315
        }
      }
316
    else if (m_PxType->pixType == GDT_CInt16)
317
      {
318
      //std::cout << "Convert input File from GDT_CInt16 to GDT_CFloat32" << std::endl;
319
      typedef std::complex<short>            ComplexShortType;
320

321 322
      for (unsigned int itPxl = 0; itPxl < (unsigned int) (nbPixelToRead * m_NbBands); itPxl++)
        {
323
        ComplexShortType pxlValue = *(static_cast<ComplexShortType*>( static_cast<void*>(pBufferTemp)) + itPxl );
324

325
        RealType    pxlValueReal( static_cast<ScalarRealType>(pxlValue.real()), static_cast<ScalarRealType>(pxlValue.imag()) );
326

327
        memcpy((void*) (&(p[itPxl*sizeof(RealType)])), (const void*) (&(pxlValueReal)), (size_t) (sizeof(RealType)));
328 329
        }
      }
330 331
    //std::cout << "CONVERTED BUFFER:" <<std::endl;
    //printDataBuffer(p, GDT_CFloat64, m_NbBands, lNbColumns*lNbLines);
332
    delete[] pBufferTemp;
333
    }
334

335 336
  // In the indexed case, one has to retrieve the index image and the
  // color table, and translate p to a 4 components color values buffer
337
  else*/ if (m_IsIndexed)
OTB Bot's avatar
STYLE  
OTB Bot committed
338
    {
339 340 341
    // TODO: This is a very special case and seems to be working only
    // for unsigned char pixels. There might be a gdal method to do
    // the work in a cleaner way
342 343
    std::streamoff lNbPixels = (static_cast<std::streamoff>(lNbColumnsRegion))
                             * (static_cast<std::streamoff>(lNbLinesRegion));
344 345 346 347 348
    std::streamoff lBufferSize = static_cast<std::streamoff>(m_BytePerPixel) * lNbPixels;
    itk::VariableLengthVector<unsigned char> value(lBufferSize);

   std::streamoff step = static_cast<std::streamoff>(this->GetNumberOfComponents())
                       * static_cast<std::streamoff>(m_BytePerPixel);
349

350
    CPLErr lCrGdal = dataset->GetRasterBand(1)->RasterIO(GF_Read,
OTB Bot's avatar
STYLE  
OTB Bot committed
351 352 353 354
                                     lFirstColumn,
                                     lFirstLine,
                                     lNbColumns,
                                     lNbLines,
355
                                     const_cast<unsigned char*>(value.GetDataPointer()),
356 357
                                     lNbColumnsRegion,
                                     lNbLinesRegion,
358
                                     m_PxType->pixType,
OTB Bot's avatar
STYLE  
OTB Bot committed
359 360
                                     0,
                                     0);
361
    if (lCrGdal == CE_Failure)
OTB Bot's avatar
STYLE  
OTB Bot committed
362
      {
OTB Bot's avatar
STYLE  
OTB Bot committed
363
      itkExceptionMacro(<< "Error while reading image (GDAL format) '"
364
        << m_FileName.c_str() << "' : " << CPLGetLastErrorMsg());
OTB Bot's avatar
STYLE  
OTB Bot committed
365
      }
366 367
    // Interpret index as color
    std::streamoff cpt(0);
368
    GDALColorTable* colorTable = dataset->GetRasterBand(1)->GetColorTable();
369
    for (std::streamoff i = 0; i < lBufferSize; i = i + static_cast<std::streamoff>(m_BytePerPixel))
OTB Bot's avatar
STYLE  
OTB Bot committed
370
      {
371
      GDALColorEntry color;
372
      colorTable->GetColorEntryAsRGB(value[i], &color);
373
      p[cpt] = color.c1;
OTB Bot's avatar
STYLE  
OTB Bot committed
374 375 376
      p[cpt + 1] = color.c2;
      p[cpt + 2] = color.c3;
      p[cpt + 3] = color.c4;
377
      cpt += step;
OTB Bot's avatar
STYLE  
OTB Bot committed
378 379
      }
    }
380
  else
OTB Bot's avatar
STYLE  
OTB Bot committed
381
    {
382
    /********  Nominal case ***********/
383
    int pixelOffset = m_BytePerPixel * m_NbBands;
384
    int lineOffset  = m_BytePerPixel * m_NbBands * lNbColumnsRegion;
385 386
    int bandOffset  = m_BytePerPixel;
    int nbBands     = m_NbBands;
387 388

    // In some cases, we need to change some parameters for RasterIO
389
    if(!GDALDataTypeIsComplex(m_PxType->pixType) && m_IsComplex && m_IsVectorImage && (m_NbBands > 1))
390
      {
OTB Bot's avatar
STYLE  
OTB Bot committed
391
      pixelOffset = m_BytePerPixel * 2;
392
      lineOffset  = pixelOffset * lNbColumnsRegion;
393
      bandOffset  = m_BytePerPixel;
394
      }
395 396

    // keep it for the moment
397 398
    //otbMsgDevMacro(<< "Number of bands inside input file: " << m_NbBands);
    otbMsgDevMacro(<< "Parameters RasterIO : \n"
Julien Malik's avatar
Julien Malik committed
399 400 401 402
                   << " indX = " << lFirstColumn << "\n"
                   << " indY = " << lFirstLine << "\n"
                   << " sizeX = " << lNbColumns << "\n"
                   << " sizeY = " << lNbLines << "\n"
403 404
                   << " Buffer Size X = " << lNbColumnsRegion << "\n"
                   << " Buffer Size Y = " << lNbLinesRegion << "\n"
405
                   << " GDAL Data Type = " << GDALGetDataTypeName(m_PxType->pixType) << "\n"
406
                   << " nbBands = " << nbBands << "\n"
Julien Malik's avatar
Julien Malik committed
407 408
                   << " pixelOffset = " << pixelOffset << "\n"
                   << " lineOffset = " << lineOffset << "\n"
409
                   << " bandOffset = " << bandOffset );
Julien Malik's avatar
Julien Malik committed
410

411
    otb::Stopwatch chrono = otb::Stopwatch::StartNew();
412 413 414 415 416
    CPLErr lCrGdal = m_Dataset->GetDataSet()->RasterIO(GF_Read,
                                                       lFirstColumn,
                                                       lFirstLine,
                                                       lNbColumns,
                                                       lNbLines,
417
                                                       p,
418 419
                                                       lNbColumnsRegion,
                                                       lNbLinesRegion,
420
                                                       m_PxType->pixType,
421
                                                       nbBands,
422
                                                       // We want to read all bands
423
                                                       ITK_NULLPTR,
424 425 426
                                                       pixelOffset,
                                                       lineOffset,
                                                       bandOffset);
Julien Malik's avatar
Julien Malik committed
427
    chrono.Stop();
428
    otbMsgDevMacro(<< "RasterIO Read took " << chrono.GetElapsedMilliseconds() << " ms")
429

430 431
    // Check if gdal call succeed
    if (lCrGdal == CE_Failure)
432
      {
OTB Bot's avatar
STYLE  
OTB Bot committed
433
      itkExceptionMacro(<< "Error while reading image (GDAL format) '"
434
        << m_FileName.c_str() << "' : " << CPLGetLastErrorMsg());
435
      return;
436
      }
437
    //printDataBuffer(p, m_PxType->pixType, m_NbBands, lNbColumnsRegion*lNbLinesRegion);
438
    }
439 440
}

441 442
bool GDALImageIO::GetSubDatasetInfo(std::vector<std::string> &names, std::vector<std::string> &desc)
{
443
  // Note: we assume that the subdatasets are in order : SUBDATASET_ID_NAME, SUBDATASET_ID_DESC, SUBDATASET_ID+1_NAME, SUBDATASET_ID+1_DESC
444 445 446
  char** papszMetadata;
  papszMetadata = m_Dataset->GetDataSet()->GetMetadata("SUBDATASETS");

447 448 449 450
  // Have we find some dataSet ?
  // This feature is supported only for hdf4 and hdf5 file (regards to the bug 270)
  if ( (CSLCount(papszMetadata) > 0) &&
       ( (strcmp(m_Dataset->GetDataSet()->GetDriver()->GetDescription(),"HDF4") == 0) ||
451 452
         (strcmp(m_Dataset->GetDataSet()->GetDriver()->GetDescription(),"HDF5") == 0) ||
	 (strcmp(m_Dataset->GetDataSet()->GetDriver()->GetDescription(),"SENTINEL2") == 0) ) )
453
    {
454
    for (int cpt = 0; papszMetadata[cpt] != ITK_NULLPTR; ++cpt)
455 456 457 458 459 460 461 462 463 464 465 466 467 468 469 470 471 472
      {
      std::string key, name;
      if (System::ParseHdfSubsetName(papszMetadata[cpt], key, name))
        {
        otbMsgDevMacro(<< "- key:  " << key);
        otbMsgDevMacro(<< "- name: " << name);
        // check if this is a dataset name
        if (key.find("_NAME") != std::string::npos) names.push_back(name);
        // check if this is a dataset descriptor
        if (key.find("_DESC") != std::string::npos) desc.push_back(name);
        }
      }
    }
  else
    {
    return false;
    }
  if (names.empty() || desc.empty()) return false;
473 474 475 476 477 478
  if (names.size() != desc.size())
    {
    names.clear();
    desc.clear();
    return false;
    }
479 480

  return true;
481 482
}

483 484 485 486 487
bool GDALImageIO::GDALPixelTypeIsComplex()
{
  return GDALDataTypeIsComplex(m_PxType->pixType);
}

488 489 490 491 492 493
void GDALImageIO::ReadImageInformation()
{
  //std::ifstream file;
  this->InternalReadImageInformation();
}

494
unsigned int GDALImageIO::GetOverviewsCount()
495
{
496
  GDALDataset* dataset = m_Dataset->GetDataSet();
497

498
  // JPEG2000 case : use the number of overviews actually in the dataset
499
  if (m_Dataset->IsJPEG2000())
500
    {
501 502
    // Include the full resolution in overviews count
    return dataset->GetRasterBand(1)->GetOverviewCount()+1;
503
    }
504

505
  if (dataset->GetRasterBand(1)->GetOverviewCount())
506 507
    // Include the full resolution in overviews count
    return dataset->GetRasterBand(1)->GetOverviewCount()+1;
508

509
  // default case: compute overviews until one of the dimensions is 1
510
  bool flagStop = false;
511
  unsigned int possibleOverviewCount = 0;
512 513
  while (!flagStop)
    {
514 515
    unsigned int tDimX = uint_ceildivpow2(dataset->GetRasterXSize(),possibleOverviewCount);
    unsigned int tDimY = uint_ceildivpow2(dataset->GetRasterYSize(),possibleOverviewCount);
516

517
    possibleOverviewCount++;
518 519 520 521 522
    if ( (tDimX == 1) || (tDimY == 1) )
      {
      flagStop = true;
      }
    }
523
  return possibleOverviewCount;
524 525 526
}


527 528 529
std::vector<std::string> GDALImageIO::GetOverviewsInfo()
{
  std::vector<std::string> desc;
530 531 532

  // This should never happen, according to implementation of GetOverviewCount()
  if (this->GetOverviewsCount() == 0)
533
    return desc;
534 535

  std::ostringstream oss;
536

537 538
  // If gdal exposes actual overviews
  unsigned int lOverviewsCount = m_Dataset->GetDataSet()->GetRasterBand(1)->GetOverviewCount();
539

540
  if (lOverviewsCount)
541
    {
542 543
    unsigned int x = m_OriginalDimensions[0];
    unsigned int y = m_OriginalDimensions[1];
544

545 546
    oss.str("");
    oss << "Resolution: 0 (Image [w x h]: " << x << "x" << y << ")";
547
    desc.push_back(oss.str());
548

549 550 551 552 553 554 555 556
    for( unsigned int iOverview = 0; iOverview < lOverviewsCount; iOverview++ )
      {
      x = m_Dataset->GetDataSet()->GetRasterBand(1)->GetOverview(iOverview)->GetXSize();
      y = m_Dataset->GetDataSet()->GetRasterBand(1)->GetOverview(iOverview)->GetYSize();
      oss.str("");
      oss << "Resolution: " << iOverview+1 << " (Image [w x h]: " << x << "x" << y << ")";
      desc.push_back(oss.str());
      }
557
    }
558 559 560 561
  else
    {
    // Fall back to gdal implicit overviews
    lOverviewsCount = this->GetOverviewsCount();
562

563
    unsigned int originalWidth = m_OriginalDimensions[0];
564 565
    unsigned int originalHeight = m_OriginalDimensions[1];

566 567 568 569 570 571 572 573 574 575 576
    // Get the overview sizes
    for( unsigned int iOverview = 0; iOverview < lOverviewsCount; iOverview++ )
      {
      // For each resolution we will compute the tile dim and image dim
      unsigned int w = uint_ceildivpow2( originalWidth, iOverview);
      unsigned int h = uint_ceildivpow2( originalHeight, iOverview);
      oss.str("");
      oss << "Resolution: " << iOverview << " (Image [w x h]: " << w << "x" << h << ")";
      desc.push_back(oss.str());
      }
    }
577

578
  return desc;
579 580
}

581 582
void GDALImageIO::InternalReadImageInformation()
{
583 584 585
  itk::ExposeMetaData<unsigned int>(this->GetMetaDataDictionary(),
                                    MetaDataKey::ResolutionFactor,
                                    m_ResolutionFactor);
586

587 588 589
  itk::ExposeMetaData<unsigned int>(this->GetMetaDataDictionary(),
                                    MetaDataKey::SubDatasetIndex,
                                    m_DatasetNumber);
590

Emmanuel Christophe's avatar
Emmanuel Christophe committed
591 592 593 594 595 596
  // Detecting if we are in the case of an image with subdatasets
  // example: hdf Modis data
  // in this situation, we are going to change the filename to the
  // supported gdal format using the m_DatasetNumber value
  // HDF4_SDS:UNKNOWN:"myfile.hdf":2
  // and make m_Dataset point to it.
597 598
  if (m_Dataset->GetDataSet()->GetRasterCount() == 0)
    {
Emmanuel Christophe's avatar
Emmanuel Christophe committed
599 600
    // this happen in the case of a hdf file with SUBDATASETS
    // Note: we assume that the datasets are in order
601 602
    char** papszMetadata;
    papszMetadata = m_Dataset->GetDataSet()->GetMetadata("SUBDATASETS");
Emmanuel Christophe's avatar
Emmanuel Christophe committed
603
    //TODO: we might want to keep the list of names somewhere, at least the number of datasets
604 605 606
    std::vector<std::string> names;
    if( CSLCount(papszMetadata) > 0 )
      {
607
      for( int cpt = 0; papszMetadata[cpt] != ITK_NULLPTR; ++cpt )
608 609
        {
        std::string key, name;
610
        if (System::ParseHdfSubsetName(papszMetadata[cpt], key, name))
611 612 613 614 615 616 617 618 619 620 621 622 623 624 625
          {
          otbMsgDevMacro(<< "- key:  " << key);
          otbMsgDevMacro(<< "- name: " << name);
          // check if this is a dataset name
          if (key.find("_NAME") != std::string::npos) names.push_back(name);
          }
        }
      }
    if (m_DatasetNumber < names.size())
      {
      otbMsgDevMacro(<< "Reading: " << names[m_DatasetNumber]);
      m_Dataset = GDALDriverManagerWrapper::GetInstance().Open(names[m_DatasetNumber]);
      }
    else
      {
626
      itkExceptionMacro(<< "Dataset requested does not exist (" << names.size() << " datasets)");
627 628 629
      }
    }

630
  GDALDataset* dataset = m_Dataset->GetDataSet();
631

632 633
  // Get image dimensions
  if ( dataset->GetRasterXSize() == 0 || dataset->GetRasterYSize() == 0 )
Thomas Feuvrier's avatar
nomsg  
Thomas Feuvrier committed
634
    {
OTB Bot's avatar
STYLE  
OTB Bot committed
635
    itkExceptionMacro(<< "Dimension is undefined.");
Thomas Feuvrier's avatar
nomsg  
Thomas Feuvrier committed
636
    }
637 638

  // Set image dimensions into IO
639 640 641 642 643 644 645 646
  m_Dimensions[0] = uint_ceildivpow2(dataset->GetRasterXSize(),m_ResolutionFactor);
  m_Dimensions[1] = uint_ceildivpow2(dataset->GetRasterYSize(),m_ResolutionFactor);

  // Keep the original dimension of the image
  m_OriginalDimensions.push_back(dataset->GetRasterXSize());
  m_OriginalDimensions.push_back(dataset->GetRasterYSize());

  otbMsgDevMacro(<< "Original Dimensions of the input file: " << m_OriginalDimensions[0] << " x " << m_OriginalDimensions[1]);
647

OTB Bot's avatar
STYLE  
OTB Bot committed
648
  // Get Number of Bands
649
  m_NbBands = dataset->GetRasterCount();
650

651
  // Get the number of overviews of the file (based on the first band)
652
  m_NumberOfOverviews = dataset->GetRasterBand(1)->GetOverviewCount();
653

654
  // Get the overview sizes
655 656
  for( unsigned int iOverview = 0; iOverview < m_NumberOfOverviews; iOverview++ )
  {
657 658 659 660 661 662
      std::pair <unsigned int, unsigned int> tempSize;
      tempSize.first  = GDALGetRasterBandXSize( dataset->GetRasterBand(1)->GetOverview(iOverview) );
      tempSize.second = GDALGetRasterBandYSize( dataset->GetRasterBand(1)->GetOverview(iOverview) );
      m_OverviewsSize.push_back(tempSize);

      /*std::cout << "Overviews size of input file" << m_FileName << ": "
OTB Bot's avatar
STYLE  
OTB Bot committed
663
                <<  m_OverviewsSize.back().first << " x " << m_OverviewsSize.back().second <<   std::endl; */
664 665
      otbMsgDevMacro( << "Overviews size of input file" << m_FileName << ": "
                      <<  m_OverviewsSize.back().first << " x " << m_OverviewsSize.back().second);
666 667
  }

668
  otbMsgDevMacro(<< "Number of Overviews inside input file: " << m_NumberOfOverviews);
669 670
  otbMsgDevMacro(<< "Input file dimension: " << m_Dimensions[0] << ", " << m_Dimensions[1]);
  otbMsgDevMacro(<< "Number of bands inside input file: " << m_NbBands);
671

OTB Bot's avatar
STYLE  
OTB Bot committed
672
  this->SetNumberOfComponents(m_NbBands);
Guillaume Borrut's avatar
Guillaume Borrut committed
673

OTB Bot's avatar
STYLE  
OTB Bot committed
674 675
  // Set the number of dimensions (verify for the dim )
  this->SetNumberOfDimensions(2);
Thomas Feuvrier's avatar
nomsg  
Thomas Feuvrier committed
676

677
  otbMsgDevMacro(<< "Nb of Dimensions of the input file: " << m_NumberOfDimensions);
Thomas Feuvrier's avatar
nomsg  
Thomas Feuvrier committed
678

OTB Bot's avatar
STYLE  
OTB Bot committed
679 680
  // Automatically set the Type to Binary for GDAL data
  this->SetFileTypeToBinary();
681

OTB Bot's avatar
STYLE  
OTB Bot committed
682
  // Get Data Type
683
  // Consider only the data type given by the first band
OTB Bot's avatar
STYLE  
OTB Bot committed
684
  // Maybe be could changed (to check)
685 686 687
  m_PxType->pixType = dataset->GetRasterBand(1)->GetRasterDataType();
  otbMsgDevMacro(<< "PixelType inside input file: "<< GDALGetDataTypeName(m_PxType->pixType) );
  if (m_PxType->pixType == GDT_Byte)
688
    {
OTB Bot's avatar
STYLE  
OTB Bot committed
689
    SetComponentType(UCHAR);
690
    }
691
  else if (m_PxType->pixType == GDT_UInt16)
692
    {
OTB Bot's avatar
STYLE  
OTB Bot committed
693
    SetComponentType(USHORT);
694
    }
695
  else if (m_PxType->pixType == GDT_Int16)
696
    {
OTB Bot's avatar
STYLE  
OTB Bot committed
697
    SetComponentType(SHORT);
698
    }
699
  else if (m_PxType->pixType == GDT_UInt32)
700
    {
OTB Bot's avatar
STYLE  
OTB Bot committed
701
    SetComponentType(UINT);
702
    }
703
  else if (m_PxType->pixType == GDT_Int32)
704
    {
OTB Bot's avatar
STYLE  
OTB Bot committed
705
    SetComponentType(INT);
706
    }
707
  else if (m_PxType->pixType == GDT_Float32)
708
    {
OTB Bot's avatar
STYLE  
OTB Bot committed
709
    SetComponentType(FLOAT);
710
    }
711
  else if (m_PxType->pixType == GDT_Float64)
712
    {
OTB Bot's avatar
STYLE  
OTB Bot committed
713
    SetComponentType(DOUBLE);
714
    }
715
  else if (m_PxType->pixType == GDT_CInt16)
716 717 718
    {
    SetComponentType(CSHORT);
    }
719
  else if (m_PxType->pixType == GDT_CInt32)
720 721 722
    {
    SetComponentType(CINT);
    }
723
  else if (m_PxType->pixType == GDT_CFloat32)
724 725 726
    {
    SetComponentType(CFLOAT);
    }
727
  else if (m_PxType->pixType == GDT_CFloat64)
728 729 730
    {
    SetComponentType(CDOUBLE);
    }
OTB Bot's avatar
STYLE  
OTB Bot committed
731
  else
732
    {
OTB Bot's avatar
STYLE  
OTB Bot committed
733
    itkExceptionMacro(<< "Pixel type unknown");
734
    }
Thomas Feuvrier's avatar
nomsg  
Thomas Feuvrier committed
735

OTB Bot's avatar
STYLE  
OTB Bot committed
736
  if (this->GetComponentType() == CHAR)
737
    {
738
    m_BytePerPixel = 1;
739
    }
OTB Bot's avatar
STYLE  
OTB Bot committed
740
  else if (this->GetComponentType() == UCHAR)
741
    {
742
    m_BytePerPixel = 1;
743
    }
OTB Bot's avatar
STYLE  
OTB Bot committed
744
  else if (this->GetComponentType() == USHORT)
745
    {
746
    m_BytePerPixel = 2;
747
    }
OTB Bot's avatar
STYLE  
OTB Bot committed
748
  else if (this->GetComponentType() == SHORT)
749
    {
750
    m_BytePerPixel = 2;
751
    }
OTB Bot's avatar
STYLE  
OTB Bot committed
752
  else if (this->GetComponentType() == INT)
753
    {
754
    m_BytePerPixel = 4;
755
    }
OTB Bot's avatar
STYLE  
OTB Bot committed
756
  else if (this->GetComponentType() == UINT)
757
    {
758
    m_BytePerPixel = 4;
759
    }
760 761
  else if (this->GetComponentType() == LONG)
    {
762
    m_BytePerPixel = sizeof(long);
763 764 765
    }
  else if (this->GetComponentType() == ULONG)
    {
766
    m_BytePerPixel = sizeof(unsigned long);
767
    }
OTB Bot's avatar
STYLE  
OTB Bot committed
768
  else if (this->GetComponentType() == FLOAT)
769
    {
770
    m_BytePerPixel = 4;
771
    }
OTB Bot's avatar
STYLE  
OTB Bot committed
772
  else if (this->GetComponentType() == DOUBLE)
773
    {
774
    m_BytePerPixel = 8;
775
    }
776 777 778 779 780 781 782 783
  else if (this->GetComponentType() == CSHORT)
    {
    m_BytePerPixel = sizeof(std::complex<short>);
    }
  else if (this->GetComponentType() == CINT)
    {
    m_BytePerPixel = sizeof(std::complex<int>);
    }
784
  else if (this->GetComponentType() == CFLOAT)
785
    {
786
    /*if (m_PxType->pixType == GDT_CInt16)
787
      m_BytePerPixel = sizeof(std::complex<short>);
788
    else if (m_PxType->pixType == GDT_CInt32)
789
      m_BytePerPixel = sizeof(std::complex<int>);
790
    else*/
791 792 793 794
      m_BytePerPixel = sizeof(std::complex<float>);
    }
  else if (this->GetComponentType() == CDOUBLE)
    {
795
      m_BytePerPixel = sizeof(std::complex<double>);
796
    }
OTB Bot's avatar
STYLE  
OTB Bot committed
797
  else
798
    {
OTB Bot's avatar
STYLE  
OTB Bot committed
799
    itkExceptionMacro(<< "Component type unknown");
800
    }
801

OTB Bot's avatar
STYLE  
OTB Bot committed
802
  /******************************************************************/
803
  // Set the pixel type with some special cases linked to the fact
804
  //  we read some data with complex type.
805
  if ( GDALDataTypeIsComplex(m_PxType->pixType) ) // Try to read data with complex type with GDAL
806
    {
807 808 809 810 811
    if ( !m_IsComplex && m_IsVectorImage )
      {
      // we are reading a complex data set into an image where the pixel
      // type is Vector<real>: we have to double the number of component
      // for that to work
812
      otbMsgDevMacro( << "GDALtypeIO= Complex and IFReader::InternalPixelType= Scalar and IFReader::PixelType= Vector");
813 814 815 816 817 818 819
      this->SetNumberOfComponents(m_NbBands*2);
      this->SetPixelType(VECTOR);
      }
    else
      {
      this->SetPixelType(COMPLEX);
      }
820
    }
821
  else // Try to read data with scalar type with GDAL
822
    {
OTB Bot's avatar
STYLE  
OTB Bot committed
823 824
    this->SetNumberOfComponents(m_NbBands);
    if (this->GetNumberOfComponents() == 1)
825
      {
OTB Bot's avatar
STYLE  
OTB Bot committed
826
      this->SetPixelType(SCALAR);
827
      }
OTB Bot's avatar
STYLE  
OTB Bot committed
828
    else
829
      {
OTB Bot's avatar
STYLE  
OTB Bot committed
830
      this->SetPixelType(VECTOR);
831
      }
Thomas Feuvrier's avatar
nomsg  
Thomas Feuvrier committed
832
    }
833

834
  /*** Parameters set by Internal Read function ***/
835 836 837 838
  otbMsgDevMacro( << "Pixel Type IFReader = " << GetPixelTypeAsString(this->GetPixelType()) )
  otbMsgDevMacro( << "Number of component IFReader = " << this->GetNumberOfComponents() )
  otbMsgDevMacro( << "Byte per pixel set = " << m_BytePerPixel )
  otbMsgDevMacro( << "Component Type set = " << GetComponentTypeAsString(this->GetComponentType()) );
839

840 841 842
  /*----------------------------------------------------------------------*/
  /*-------------------------- METADATA ----------------------------------*/
  /*----------------------------------------------------------------------*/
843

844
  // Now initialize the itk dictionary
OTB Bot's avatar
STYLE  
OTB Bot committed
845
  itk::MetaDataDictionary& dict = this->GetMetaDataDictionary();
846

847 848 849 850 851 852
  // Report the typical block size if possible
  if (dataset->GetRasterCount() > 0)
  {
    int blockSizeX = 0;
    int blockSizeY = 0;

OTB Bot's avatar
STYLE  
OTB Bot committed
853
    dataset->GetRasterBand(1)->GetBlockSize(&blockSizeX, &blockSizeY);
854

855 856
    if(blockSizeX > 0 && blockSizeY > 0)
      {
857 858 859
      otbMsgDevMacro(<< "Original blockSize: "<< blockSizeX << " x " << blockSizeY );

      blockSizeX = uint_ceildivpow2(blockSizeX,m_ResolutionFactor);
860
      if (m_Dataset->IsJPEG2000())
861 862 863 864 865 866 867 868 869
        {
        // Jpeg2000 case : use the real block size Y
        blockSizeY = uint_ceildivpow2(blockSizeY,m_ResolutionFactor);
        }
      else
        {
        // Try to keep the GDAL block memory constant
        blockSizeY = blockSizeY * (1 << m_ResolutionFactor);
        }
870 871 872

      otbMsgDevMacro(<< "Decimated blockSize: "<< blockSizeX << " x " << blockSizeY );

873 874 875 876 877 878
      itk::EncapsulateMetaData<unsigned int>(dict, MetaDataKey::TileHintX, blockSizeX);
      itk::EncapsulateMetaData<unsigned int>(dict, MetaDataKey::TileHintY, blockSizeY);
      }
  }


879
  /* -------------------------------------------------------------------- */
Guillaume Borrut's avatar
Guillaume Borrut committed
880
  /*  Get Spacing                */
881
  /* -------------------------------------------------------------------- */
882

883
  // Default Spacing
OTB Bot's avatar
STYLE  
OTB Bot committed
884 885
  m_Spacing[0] = 1;
  m_Spacing[1] = 1;
886

887 888 889 890
  // Reset origin to GDAL convention default
  m_Origin[0] = 0.0;
  m_Origin[1] = 0.0;

891 892
  // flag to detect images in sensor geometry
  bool isSensor = false;
893

OTB Bot's avatar
STYLE  
OTB Bot committed
894
  if (m_NumberOfDimensions == 3) m_Spacing[2] = 1;
895

896
  char** papszMetadata = dataset->GetMetadata(ITK_NULLPTR);
897

898 899 900
  /* -------------------------------------------------------------------- */
  /*      Report general info.                                            */
  /* -------------------------------------------------------------------- */
OTB Bot's avatar
STYLE  
OTB Bot committed
901
  GDALDriverH hDriver;
902

903
  hDriver = dataset->GetDriver();
904

OTB Bot's avatar
STYLE  
OTB Bot committed
905
  std::string driverShortName =  static_cast<std::string>(GDALGetDriverShortName(hDriver));
906
  std::string driverLongName  =  static_cast<std::string>(GDALGetDriverLongName(hDriver));
907

OTB Bot's avatar
STYLE  
OTB Bot committed
908
  itk::EncapsulateMetaData<std::string>(dict, MetaDataKey::DriverShortNameKey, driverShortName);
909
  itk::EncapsulateMetaData<std::string>(dict, MetaDataKey::DriverLongNameKey,  driverLongName);
910

911
  if (m_Dataset->IsJPEG2000())
912 913 914 915
    {
    // store the cache size used for Jpeg2000 files
    itk::EncapsulateMetaData<unsigned int>(dict, MetaDataKey::CacheSizeInBytes , GDALGetCacheMax64());
    }
Caroline Ruffel's avatar
nomsg  
Caroline Ruffel committed