Orfeo Toolbox  4.0
otbGDALImageIO.cxx
Go to the documentation of this file.
1 /*=========================================================================
2 
3  Program: ORFEO Toolbox
4  Language: C++
5  Date: $Date$
6  Version: $Revision$
7 
8 
9  Copyright (c) Centre National d'Etudes Spatiales. All rights reserved.
10  See OTBCopyright.txt for details.
11 
12 
13  This software is distributed WITHOUT ANY WARRANTY; without even
14  the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
15  PURPOSE. See the above copyright notices for more information.
16 
17 =========================================================================*/
18 #include <iostream>
19 #include <fstream>
20 #include <vector>
21 
22 #include "otbGDALImageIO.h"
23 #include "otbMacro.h"
24 #include "otbSystem.h"
25 #include "itksys/SystemTools.hxx"
26 #include "otbImage.h"
28 
29 #include "itkMetaDataObject.h"
30 #include "otbMetaDataKey.h"
31 
32 #include "itkRGBPixel.h"
33 #include "itkRGBAPixel.h"
34 #include "itkTimeProbe.h"
35 
36 #include "cpl_conv.h"
37 #include "ogr_spatialref.h"
38 #include "ogr_srs_api.h"
39 
41 
42 #include <boost/algorithm/string/predicate.hpp>
43 #include "otbOGRHelpers.h"
44 
45 inline unsigned int uint_ceildivpow2(unsigned int a, unsigned int b) {
46  return (a + (1 << b) - 1) >> b;
47 }
48 
49 namespace otb
50 {
51 
53 {
54 public:
55  GDALDataTypeWrapper() : pixType(GDT_Byte) {}
58  {
59  pixType = w.pixType;
60  }
62  {
63  pixType = w.pixType;
64  return *this;
65  }
66  GDALDataType pixType;
67 }; // end of GDALDataTypeWrapper
68 
69 
70 /*
71 template<class InputType>
72 void printOutputData(InputType *pData, int nbBands, int nbPixelToRead)
73 {
74  for (unsigned int itPxl = 0; itPxl < (unsigned int) (nbPixelToRead * nbBands); itPxl++)
75  {
76  std::cout << "Buffer["<< itPxl << "] = " << *(pData + itPxl) << std::endl;
77  }
78 };
79 
80 void printDataBuffer(unsigned char *pData, GDALDataType pxlType, int nbBands, int nbPixelToRead)
81 {
82  if (pxlType == GDT_Int16)
83  {
84  printOutputData( static_cast<short*>( static_cast<void*>(pData) ), nbBands, nbPixelToRead);
85  }
86  else if (pxlType == GDT_Int32)
87  {
88  printOutputData( static_cast<int*>( static_cast<void*>(pData) ), nbBands, nbPixelToRead);
89  }
90  else if (pxlType == GDT_Float32)
91  {
92  printOutputData( static_cast<float*>( static_cast<void*>(pData) ), nbBands, nbPixelToRead);
93  }
94  else if (pxlType == GDT_Float64)
95  {
96  printOutputData( static_cast<double*>( static_cast<void*>(pData) ), nbBands, nbPixelToRead);
97  }
98  else if (pxlType == GDT_CInt16)
99  {
100  printOutputData( static_cast<std::complex<short>*>( static_cast<void*>(pData) ), nbBands, nbPixelToRead);
101  }
102  else if (pxlType == GDT_CInt32)
103  {
104  printOutputData( static_cast<std::complex<int>*>( static_cast<void*>(pData) ), nbBands, nbPixelToRead);
105  }
106  else if (pxlType == GDT_CFloat32)
107  {
108  printOutputData( static_cast<std::complex<float>*>( static_cast<void*>(pData) ), nbBands, nbPixelToRead);
109  }
110  else if (pxlType == GDT_CFloat64)
111  {
112  printOutputData( static_cast<std::complex<double>*>( static_cast<void*>(pData) ), nbBands, nbPixelToRead);
113  }
114  else
115  {
116  std::cerr << "Pixel type unknown" << std::endl;
117  }
118 };
119 */
120 
122 {
123  // By default set number of dimensions to two.
124  this->SetNumberOfDimensions(2);
125 
126  // By default set pixel type to scalar.
128 
129  // By default set component type to unsigned char
131  m_UseCompression = false;
132  m_CompressionLevel = 4; // Range 0-9; 0 = no file compression, 9 = maximum file compression
133 
134  // Set default spacing to one
135  m_Spacing[0] = 1.0;
136  m_Spacing[1] = 1.0;
137  // Set default origin to zero
138  m_Origin[0] = 0.0;
139  m_Origin[1] = 0.0;
140 
141  m_IsIndexed = false;
142  m_DatasetNumber = 0;
143  //m_poBands = NULL;
144  //m_hDriver = NULL;
145  //m_poDataset = NULL;
146 
147  m_NbBands = 0;
149 
150  m_CanStreamWrite = false;
151  m_IsComplex = false;
152  m_IsVectorImage = false;
153 
155 
157  m_ResolutionFactor = 0;
158 
159 }
160 
162 {
163  delete m_PxType;
164 }
165 
166 // Tell only if the file can be read with GDAL.
167 bool GDALImageIO::CanReadFile(const char* file)
168 {
169  // First check the extension
170  if (file == NULL)
171  {
172  itkDebugMacro(<< "No filename specified.");
173  return false;
174  }
176  return m_Dataset.IsNotNull();
177 }
178 
179 // Used to print information about this object
180 void GDALImageIO::PrintSelf(std::ostream& os, itk::Indent indent) const
181 {
182  Superclass::PrintSelf(os, indent);
183  os << indent << "Compression Level : " << m_CompressionLevel << "\n";
184  os << indent << "IsComplex (otb side) : " << m_IsComplex << "\n";
185  os << indent << "Byte per pixel : " << m_BytePerPixel << "\n";
186 }
187 
188 // Read a 3D image (or event more bands)... not implemented yet
190 {
191 }
192 
193 // Read image with GDAL
194 void GDALImageIO::Read(void* buffer)
195 {
196  // Convert buffer from void * to unsigned char *
197  unsigned char *p = static_cast<unsigned char *>(buffer);
198 
199  // Check if conversion succeed
200  if (p == NULL)
201  {
202  itkExceptionMacro(<< "GDAL : Bad alloc");
203  return;
204  }
205 
206  // Get the origin of the region to read
207  int lFirstLineRegion = this->GetIORegion().GetIndex()[1];
208  int lFirstColumnRegion = this->GetIORegion().GetIndex()[0];
209 
210  // Get nb. of lines and columns of the region to read
211  int lNbLinesRegion = this->GetIORegion().GetSize()[1];
212  int lNbColumnsRegion = this->GetIORegion().GetSize()[0];
213 
214  //std::cout << "OriginBuffer= " << lFirstLineRegion << " x " << lFirstColumnRegion << std::endl;
215  //std::cout << "SizeBuffer= " << lNbLinesRegion << " x " << lNbColumnsRegion << std::endl;
216 
217  // Compute the origin of the image region to read at the initial resolution
218  int lFirstLine = lFirstLineRegion * (1 << m_ResolutionFactor);
219  int lFirstColumn = lFirstColumnRegion * (1 << m_ResolutionFactor);
220 
221  //std::cout << "OriginImage= " << lFirstLine << " x " << lFirstColumn << std::endl;
222 
223  // Compute the size of the image region to read at the initial resolution
224  int lNbLines = lNbLinesRegion * (1 << m_ResolutionFactor);
225  int lNbColumns = lNbColumnsRegion * (1 << m_ResolutionFactor);
226 
227  // Check if the image region is correct
228  if (lFirstLine + lNbLines > static_cast<int>(m_OriginalDimensions[1]))
229  lNbLines = static_cast<int>(m_OriginalDimensions[1]-lFirstLine);
230  if (lFirstColumn + lNbColumns > static_cast<int>(m_OriginalDimensions[0]))
231  lNbColumns = static_cast<int>(m_OriginalDimensions[0]-lFirstColumn);
232 
233  //std::cout << "SizeImage= " << lNbLines << " x " << lNbColumns << std::endl;
234 
235 
236  GDALDataset* dataset = m_Dataset->GetDataSet();
237 
238  // This special case is due to the fact the CINT/CLONG types
239  // do not exists in ITK. In this case we only report the first band
240  // TODO This should be fixed
241  /*if (GDALDataTypeIsComplex(m_PxType->pixType)
242  && (m_PxType->pixType != GDT_CFloat32)
243  && (m_PxType->pixType != GDT_CFloat64))
244  {
245  int pixelOffset = m_BytePerPixel * m_NbBands;
246  int lineOffset = m_BytePerPixel * m_NbBands * lNbColumns;
247  int bandOffset = m_BytePerPixel;
248  int nbBands = m_NbBands;
249 
250  int nbPixelToRead = lNbColumns * lNbLines;
251  std::streamoff nbBytes = static_cast<std::streamoff>(m_NbBands) * static_cast<std::streamoff>(nbPixelToRead) * static_cast<std::streamoff>(m_BytePerPixel);
252  unsigned char *pBufferTemp = new unsigned char[static_cast<unsigned int>(nbBytes)];
253 
254  // keep it for the moment
255  otbMsgDevMacro(<< "Parameters RasterIO (case CInt and CShort):"
256  << "\n indX = " << lFirstColumn
257  << "\n indY = " << lFirstLine
258  << "\n sizeX = " << lNbColumns
259  << "\n sizeY = " << lNbLines
260  << "\n GDAL Data Type = " << GDALGetDataTypeName(m_PxType->pixType)
261  << "\n pixelOffset = " << pixelOffset
262  << "\n lineOffset = " << lineOffset
263  << "\n bandOffset = " << bandOffset);
264 
265  CPLErr lCrGdal = m_Dataset->GetDataSet()->RasterIO(GF_Read,
266  lFirstColumn,
267  lFirstLine,
268  lNbColumns,
269  lNbLines,
270  pBufferTemp, //p, // pData
271  lNbColumns,
272  lNbLines,
273  m_PxType->pixType,
274  nbBands,
275  // We want to read all bands
276  NULL,
277  pixelOffset,
278  lineOffset,
279  bandOffset);
280  // Check for gdal error
281  if (lCrGdal == CE_Failure)
282  {
283  itkExceptionMacro(<< "Error while reading image (GDAL format) " << m_FileName );
284  delete[] pBufferTemp;
285  return;
286  }
287  //std::cout << "RAW BUFFER:" <<std::endl;
288  //printDataBuffer(pBufferTemp, m_PxType->pixType, m_NbBands, lNbColumns*lNbLines);
289 
290  // Convert the buffer to GDT_Float64 type
291  typedef std::complex<float> RealType;
292  typedef double ScalarRealType;
293 
294  if (m_PxType->pixType == GDT_CInt32)
295  {
296  //std::cout << "Convert input File from GDT_CInt32 to GDT_CFloat32" << std::endl;
297  typedef std::complex<int> ComplexIntType;
298 
299  for (unsigned int itPxl = 0; itPxl < (unsigned int) (nbPixelToRead * m_NbBands); itPxl++)
300  {
301  ComplexIntType pxlValue = *(static_cast<ComplexIntType*>( static_cast<void*>(pBufferTemp)) + itPxl );
302 
303  RealType pxlValueReal( static_cast<ScalarRealType>(pxlValue.real()), static_cast<ScalarRealType>(pxlValue.imag()) );
304 
305  memcpy((void*) (&(p[itPxl*sizeof(RealType)])), (const void*) (&(pxlValueReal)), (size_t) (sizeof(RealType)));
306  }
307  }
308  else if (m_PxType->pixType == GDT_CInt16)
309  {
310  //std::cout << "Convert input File from GDT_CInt16 to GDT_CFloat32" << std::endl;
311  typedef std::complex<short> ComplexShortType;
312 
313  for (unsigned int itPxl = 0; itPxl < (unsigned int) (nbPixelToRead * m_NbBands); itPxl++)
314  {
315  ComplexShortType pxlValue = *(static_cast<ComplexShortType*>( static_cast<void*>(pBufferTemp)) + itPxl );
316 
317  RealType pxlValueReal( static_cast<ScalarRealType>(pxlValue.real()), static_cast<ScalarRealType>(pxlValue.imag()) );
318 
319  memcpy((void*) (&(p[itPxl*sizeof(RealType)])), (const void*) (&(pxlValueReal)), (size_t) (sizeof(RealType)));
320  }
321  }
322  //std::cout << "CONVERTED BUFFER:" <<std::endl;
323  //printDataBuffer(p, GDT_CFloat64, m_NbBands, lNbColumns*lNbLines);
324  delete[] pBufferTemp;
325  }
326 
327  // In the indexed case, one has to retrieve the index image and the
328  // color table, and translate p to a 4 components color values buffer
329  else*/ if (m_IsIndexed)
330  {
331  // TODO: This is a very special case and seems to be working only
332  // for unsigned char pixels. There might be a gdal method to do
333  // the work in a cleaner way
334  std::streamoff lNbPixels = (static_cast<std::streamoff>(lNbColumns))
335  * (static_cast<std::streamoff>(lNbLines));
336  std::streamoff lBufferSize = static_cast<std::streamoff>(m_BytePerPixel) * lNbPixels;
337  itk::VariableLengthVector<unsigned char> value(lBufferSize);
338 
339  std::streamoff step = static_cast<std::streamoff>(this->GetNumberOfComponents())
340  * static_cast<std::streamoff>(m_BytePerPixel);
341 
342  CPLErr lCrGdal = dataset->GetRasterBand(1)->RasterIO(GF_Read,
343  lFirstColumn,
344  lFirstLine,
345  lNbColumns,
346  lNbLines,
347  const_cast<unsigned char*>(value.GetDataPointer()),
348  lNbColumnsRegion,
349  lNbLinesRegion,
350  m_PxType->pixType,
351  0,
352  0);
353  if (lCrGdal == CE_Failure)
354  {
355  itkExceptionMacro(<< "Error while reading image (GDAL format) " << m_FileName.c_str() << ".");
356  }
357  // Interpret index as color
358  std::streamoff cpt(0);
359  GDALColorTable* colorTable = dataset->GetRasterBand(1)->GetColorTable();
360  for (std::streamoff i = 0; i < lBufferSize; i = i + static_cast<std::streamoff>(m_BytePerPixel))
361  {
362  GDALColorEntry color;
363  colorTable->GetColorEntryAsRGB(value[i], &color);
364  p[cpt] = color.c1;
365  p[cpt + 1] = color.c2;
366  p[cpt + 2] = color.c3;
367  p[cpt + 3] = color.c4;
368  cpt += step;
369  }
370  }
371  else
372  {
373  /******** Nominal case ***********/
374  int pixelOffset = m_BytePerPixel * m_NbBands;
375  int lineOffset = m_BytePerPixel * m_NbBands * lNbColumnsRegion;
376  int bandOffset = m_BytePerPixel;
377  int nbBands = m_NbBands;
378 
379  // In some cases, we need to change some parameters for RasterIO
380  if(!GDALDataTypeIsComplex(m_PxType->pixType) && m_IsComplex && m_IsVectorImage && (m_NbBands > 1))
381  {
382  pixelOffset = m_BytePerPixel * 2;
383  lineOffset = pixelOffset * lNbColumnsRegion;
384  bandOffset = m_BytePerPixel;
385  }
386 
387  // keep it for the moment
388  //otbMsgDevMacro(<< "Number of bands inside input file: " << m_NbBands);
389  otbMsgDevMacro(<< "Parameters RasterIO : \n"
390  << " indX = " << lFirstColumn << "\n"
391  << " indY = " << lFirstLine << "\n"
392  << " sizeX = " << lNbColumns << "\n"
393  << " sizeY = " << lNbLines << "\n"
394  << " Buffer Size X = " << lNbColumnsRegion << "\n"
395  << " Buffer Size Y = " << lNbLinesRegion << "\n"
396  << " GDAL Data Type = " << GDALGetDataTypeName(m_PxType->pixType) << "\n"
397  << " nbBands = " << nbBands << "\n"
398  << " pixelOffset = " << pixelOffset << "\n"
399  << " lineOffset = " << lineOffset << "\n"
400  << " bandOffset = " << bandOffset );
401 
402  itk::TimeProbe chrono;
403  chrono.Start();
404  CPLErr lCrGdal = m_Dataset->GetDataSet()->RasterIO(GF_Read,
405  lFirstColumn,
406  lFirstLine,
407  lNbColumns,
408  lNbLines,
409  p,
410  lNbColumnsRegion,
411  lNbLinesRegion,
412  m_PxType->pixType,
413  nbBands,
414  // We want to read all bands
415  NULL,
416  pixelOffset,
417  lineOffset,
418  bandOffset);
419  chrono.Stop();
420  otbMsgDevMacro(<< "RasterIO Read took " << chrono.GetTotal() << " sec")
421 
422  // Check if gdal call succeed
423  if (lCrGdal == CE_Failure)
424  {
425  itkExceptionMacro(<< "Error while reading image (GDAL format) " << m_FileName.c_str() << ".");
426  return;
427  }
428  //printDataBuffer(p, m_PxType->pixType, m_NbBands, lNbColumnsRegion*lNbLinesRegion);
429  }
430 }
431 
432 bool GDALImageIO::GetSubDatasetInfo(std::vector<std::string> &names, std::vector<std::string> &desc)
433 {
434  // Note: we assume that the subdatasets are in order : SUBDATASET_ID_NAME, SUBDATASET_ID_DESC, SUBDATASET_ID+1_NAME, SUBDATASET_ID+1_DESC
435  char** papszMetadata;
436  papszMetadata = m_Dataset->GetDataSet()->GetMetadata("SUBDATASETS");
437 
438  // Have we find some dataSet ?
439  // This feature is supported only for hdf4 and hdf5 file (regards to the bug 270)
440  if ( (CSLCount(papszMetadata) > 0) &&
441  ( (strcmp(m_Dataset->GetDataSet()->GetDriver()->GetDescription(),"HDF4") == 0) ||
442  (strcmp(m_Dataset->GetDataSet()->GetDriver()->GetDescription(),"HDF5") == 0) ) )
443  {
444  for (int cpt = 0; papszMetadata[cpt] != NULL; ++cpt)
445  {
446  std::string key, name;
447  if (System::ParseHdfSubsetName(papszMetadata[cpt], key, name))
448  {
449  otbMsgDevMacro(<< "- key: " << key);
450  otbMsgDevMacro(<< "- name: " << name);
451  // check if this is a dataset name
452  if (key.find("_NAME") != std::string::npos) names.push_back(name);
453  // check if this is a dataset descriptor
454  if (key.find("_DESC") != std::string::npos) desc.push_back(name);
455  }
456  }
457  }
458  else
459  {
460  return false;
461  }
462  if (names.empty() || desc.empty()) return false;
463  if (names.size() != desc.size())
464  {
465  names.clear();
466  desc.clear();
467  return false;
468  }
469 
470  return true;
471 }
472 
474 {
475  return GDALDataTypeIsComplex(m_PxType->pixType);
476 }
477 
479 {
480  //std::ifstream file;
482 }
483 
484 bool GDALImageIO::GetAvailableResolutions(std::vector<unsigned int>& res)
485 {
486  GDALDataset* dataset = m_Dataset->GetDataSet();
487 
488  bool flagStop = false;
489  unsigned int resFactor = 0;
490  while (!flagStop)
491  {
492  res.push_back(resFactor);
493 
494  unsigned int tDimX = uint_ceildivpow2(dataset->GetRasterXSize(),resFactor);
495  unsigned int tDimY = uint_ceildivpow2(dataset->GetRasterYSize(),resFactor);
496 
497  resFactor++;
498 
499  if ( (tDimX == 1) || (tDimY == 1) )
500  {
501  flagStop = true;
502  }
503  }
504 
505  return true;
506 }
507 
508 bool GDALImageIO::GetResolutionInfo(std::vector<unsigned int>& res, std::vector<std::string>& desc)
509 {
510  this->GetAvailableResolutions(res);
511 
512  if (res.empty())
513  return false;
514 
515  unsigned int originalWidth = m_OriginalDimensions[0];
516  unsigned int originalHeight = m_OriginalDimensions[1];
517 
518  for (std::vector<unsigned int>::iterator itRes = res.begin(); itRes < res.end(); itRes++)
519  {
520  // For each resolution we will compute the tile dim and image dim
521  std::ostringstream oss;
522 
523  unsigned int w = uint_ceildivpow2( originalWidth, *itRes);
524  unsigned int h = uint_ceildivpow2( originalHeight, *itRes);
525 
526  oss << "Resolution: " << *itRes << " (Image [w x h]: " << w << "x" << h << ", Tile [w x h]: " << "not defined x not defined" << ")";
527 
528  desc.push_back(oss.str());
529  }
530 
531  return false;
532 }
533 
535 {
536  itk::ExposeMetaData<unsigned int>(this->GetMetaDataDictionary(),
539 
540  // Detecting if we are in the case of an image with subdatasets
541  // example: hdf Modis data
542  // in this situation, we are going to change the filename to the
543  // supported gdal format using the m_DatasetNumber value
544  // HDF4_SDS:UNKNOWN:"myfile.hdf":2
545  // and make m_Dataset point to it.
546  if (m_Dataset->GetDataSet()->GetRasterCount() == 0)
547  {
548  // this happen in the case of a hdf file with SUBDATASETS
549  // Note: we assume that the datasets are in order
550  char** papszMetadata;
551  papszMetadata = m_Dataset->GetDataSet()->GetMetadata("SUBDATASETS");
552  //TODO: we might want to keep the list of names somewhere, at least the number of datasets
553  std::vector<std::string> names;
554  if( CSLCount(papszMetadata) > 0 )
555  {
556  for( int cpt = 0; papszMetadata[cpt] != NULL; ++cpt )
557  {
558  std::string key, name;
559  if (System::ParseHdfSubsetName(papszMetadata[cpt], key, name))
560  {
561  otbMsgDevMacro(<< "- key: " << key);
562  otbMsgDevMacro(<< "- name: " << name);
563  // check if this is a dataset name
564  if (key.find("_NAME") != std::string::npos) names.push_back(name);
565  }
566  }
567  }
568  if (m_DatasetNumber < names.size())
569  {
570  otbMsgDevMacro(<< "Reading: " << names[m_DatasetNumber]);
571  m_Dataset = GDALDriverManagerWrapper::GetInstance().Open(names[m_DatasetNumber]);
572  }
573  else
574  {
575  itkExceptionMacro(<< "Dataset requested does not exist (" << names.size() << " datasets)");
576  }
577  }
578 
579  GDALDataset* dataset = m_Dataset->GetDataSet();
580 
581  // Get image dimensions
582  if ( dataset->GetRasterXSize() == 0 || dataset->GetRasterYSize() == 0 )
583  {
584  itkExceptionMacro(<< "Dimension is undefined.");
585  }
586 
587  // Set image dimensions into IO
588  m_Dimensions[0] = uint_ceildivpow2(dataset->GetRasterXSize(),m_ResolutionFactor);
589  m_Dimensions[1] = uint_ceildivpow2(dataset->GetRasterYSize(),m_ResolutionFactor);
590 
591  // Keep the original dimension of the image
592  m_OriginalDimensions.push_back(dataset->GetRasterXSize());
593  m_OriginalDimensions.push_back(dataset->GetRasterYSize());
594 
595  otbMsgDevMacro(<< "Original Dimensions of the input file: " << m_OriginalDimensions[0] << " x " << m_OriginalDimensions[1]);
596 
597  // Get Number of Bands
598  m_NbBands = dataset->GetRasterCount();
599 
600  // Get the number of overviews of the file (based on the first band)
601  m_NumberOfOverviews = dataset->GetRasterBand(1)->GetOverviewCount();
602 
603  // Get the overview sizes
604  for( unsigned int iOverview = 0; iOverview < m_NumberOfOverviews; iOverview++ )
605  {
606  std::pair <unsigned int, unsigned int> tempSize;
607  tempSize.first = GDALGetRasterBandXSize( dataset->GetRasterBand(1)->GetOverview(iOverview) );
608  tempSize.second = GDALGetRasterBandYSize( dataset->GetRasterBand(1)->GetOverview(iOverview) );
609  m_OverviewsSize.push_back(tempSize);
610 
611  /*std::cout << "Overviews size of input file" << m_FileName << ": "
612  << m_OverviewsSize.back().first << " x " << m_OverviewsSize.back().second << std::endl; */
613  otbMsgDevMacro( << "Overviews size of input file" << m_FileName << ": "
614  << m_OverviewsSize.back().first << " x " << m_OverviewsSize.back().second);
615  }
616 
617  otbMsgDevMacro(<< "Number of Overviews inside input file: " << m_NumberOfOverviews);
618  otbMsgDevMacro(<< "Input file dimension: " << m_Dimensions[0] << ", " << m_Dimensions[1]);
619  otbMsgDevMacro(<< "Number of bands inside input file: " << m_NbBands);
620 
622 
623  // Set the number of dimensions (verify for the dim )
624  this->SetNumberOfDimensions(2);
625 
626  otbMsgDevMacro(<< "Nb of Dimensions of the input file: " << m_NumberOfDimensions);
627 
628  // Automatically set the Type to Binary for GDAL data
629  this->SetFileTypeToBinary();
630 
631  // Get Data Type
632  // Consider only the data type given by the first band
633  // Maybe be could changed (to check)
634  m_PxType->pixType = dataset->GetRasterBand(1)->GetRasterDataType();
635  otbMsgDevMacro(<< "PixelType inside input file: "<< GDALGetDataTypeName(m_PxType->pixType) );
636  if (m_PxType->pixType == GDT_Byte)
637  {
639  }
640  else if (m_PxType->pixType == GDT_UInt16)
641  {
643  }
644  else if (m_PxType->pixType == GDT_Int16)
645  {
647  }
648  else if (m_PxType->pixType == GDT_UInt32)
649  {
651  }
652  else if (m_PxType->pixType == GDT_Int32)
653  {
655  }
656  else if (m_PxType->pixType == GDT_Float32)
657  {
659  }
660  else if (m_PxType->pixType == GDT_Float64)
661  {
663  }
664  else if (m_PxType->pixType == GDT_CInt16)
665  {
667  }
668  else if (m_PxType->pixType == GDT_CInt32)
669  {
671  }
672  else if (m_PxType->pixType == GDT_CFloat32)
673  {
675  }
676  else if (m_PxType->pixType == GDT_CFloat64)
677  {
679  }
680  else
681  {
682  itkExceptionMacro(<< "Pixel type unknown");
683  }
684 
685  if (this->GetComponentType() == CHAR)
686  {
687  m_BytePerPixel = 1;
688  }
689  else if (this->GetComponentType() == UCHAR)
690  {
691  m_BytePerPixel = 1;
692  }
693  else if (this->GetComponentType() == USHORT)
694  {
695  m_BytePerPixel = 2;
696  }
697  else if (this->GetComponentType() == SHORT)
698  {
699  m_BytePerPixel = 2;
700  }
701  else if (this->GetComponentType() == INT)
702  {
703  m_BytePerPixel = 4;
704  }
705  else if (this->GetComponentType() == UINT)
706  {
707  m_BytePerPixel = 4;
708  }
709  else if (this->GetComponentType() == LONG)
710  {
711  m_BytePerPixel = sizeof(long);
712  }
713  else if (this->GetComponentType() == ULONG)
714  {
715  m_BytePerPixel = sizeof(unsigned long);
716  }
717  else if (this->GetComponentType() == FLOAT)
718  {
719  m_BytePerPixel = 4;
720  }
721  else if (this->GetComponentType() == DOUBLE)
722  {
723  m_BytePerPixel = 8;
724  }
725  else if (this->GetComponentType() == CSHORT)
726  {
727  m_BytePerPixel = sizeof(std::complex<short>);
728  }
729  else if (this->GetComponentType() == CINT)
730  {
731  m_BytePerPixel = sizeof(std::complex<int>);
732  }
733  else if (this->GetComponentType() == CFLOAT)
734  {
735  /*if (m_PxType->pixType == GDT_CInt16)
736  m_BytePerPixel = sizeof(std::complex<short>);
737  else if (m_PxType->pixType == GDT_CInt32)
738  m_BytePerPixel = sizeof(std::complex<int>);
739  else*/
740  m_BytePerPixel = sizeof(std::complex<float>);
741  }
742  else if (this->GetComponentType() == CDOUBLE)
743  {
744  m_BytePerPixel = sizeof(std::complex<double>);
745  }
746  else
747  {
748  itkExceptionMacro(<< "Component type unknown");
749  }
750 
751  /******************************************************************/
752  // Set the pixel type with some special cases linked to the fact
753  // we read some data with complex type.
754  if ( GDALDataTypeIsComplex(m_PxType->pixType) ) // Try to read data with complex type with GDAL
755  {
756  if ( !m_IsComplex && m_IsVectorImage )
757  {
758  // we are reading a complex data set into an image where the pixel
759  // type is Vector<real>: we have to double the number of component
760  // for that to work
761  otbMsgDevMacro( << "GDALtypeIO= Complex and IFReader::InternalPixelType= Scalar and IFReader::PixelType= Vector");
763  this->SetPixelType(VECTOR);
764  }
765  else
766  {
767  this->SetPixelType(COMPLEX);
768  }
769  }
770  else // Try to read data with scalar type with GDAL
771  {
773  if (this->GetNumberOfComponents() == 1)
774  {
775  this->SetPixelType(SCALAR);
776  }
777  else
778  {
779  this->SetPixelType(VECTOR);
780  }
781  }
782 
783  /*** Parameters set by Internal Read function ***/
784  otbMsgDevMacro( << "Pixel Type IFReader = " << GetPixelTypeAsString(this->GetPixelType()) )
785  otbMsgDevMacro( << "Number of component IFReader = " << this->GetNumberOfComponents() )
786  otbMsgDevMacro( << "Byte per pixel set = " << m_BytePerPixel )
787  otbMsgDevMacro( << "Component Type set = " << GetComponentTypeAsString(this->GetComponentType()) );
788 
789  /*----------------------------------------------------------------------*/
790  /*-------------------------- METADATA ----------------------------------*/
791  /*----------------------------------------------------------------------*/
792 
793  // Now initialize the itk dictionary
795 
796  // Report the typical block size if possible
797  if (dataset->GetRasterCount() > 0)
798  {
799  int blockSizeX = 0;
800  int blockSizeY = 0;
801 
802  dataset->GetRasterBand(1)->GetBlockSize(&blockSizeX, &blockSizeY);
803 
804  if(blockSizeX > 0 && blockSizeY > 0)
805  {
806  otbMsgDevMacro(<< "Original blockSize: "<< blockSizeX << " x " << blockSizeY );
807 
808  // Try to keep the GDAL block memory constant
809  blockSizeX = uint_ceildivpow2(blockSizeX,m_ResolutionFactor);
810  blockSizeY = blockSizeY * (1 << m_ResolutionFactor);
811 
812  otbMsgDevMacro(<< "Decimated blockSize: "<< blockSizeX << " x " << blockSizeY );
813 
814  itk::EncapsulateMetaData<unsigned int>(dict, MetaDataKey::TileHintX, blockSizeX);
815  itk::EncapsulateMetaData<unsigned int>(dict, MetaDataKey::TileHintY, blockSizeY);
816  }
817  }
818 
819 
820  /* -------------------------------------------------------------------- */
821  /* Get Spacing */
822  /* -------------------------------------------------------------------- */
823 
824  // Default Spacing
825  m_Spacing[0] = 1;
826  m_Spacing[1] = 1;
827 
828  if (m_NumberOfDimensions == 3) m_Spacing[2] = 1;
829 
830  char** papszMetadata;
831  papszMetadata = dataset->GetMetadata(NULL);
832 
833  /* -------------------------------------------------------------------- */
834  /* Report general info. */
835  /* -------------------------------------------------------------------- */
836  GDALDriverH hDriver;
837 
838  hDriver = dataset->GetDriver();
839 
840  std::string driverShortName = static_cast<std::string>(GDALGetDriverShortName(hDriver));
841  std::string driverLongName = static_cast<std::string>(GDALGetDriverLongName(hDriver));
842 
843  itk::EncapsulateMetaData<std::string>(dict, MetaDataKey::DriverShortNameKey, driverShortName);
844  itk::EncapsulateMetaData<std::string>(dict, MetaDataKey::DriverLongNameKey, driverLongName);
845 
846  /* -------------------------------------------------------------------- */
847  /* Get the projection coordinate system of the image : ProjectionRef */
848  /* -------------------------------------------------------------------- */
849  if (dataset->GetProjectionRef() != NULL && !std::string(dataset->GetProjectionRef()).empty())
850  {
851  OGRSpatialReferenceH pSR = OSRNewSpatialReference(NULL);
852 
853  const char * pszProjection = NULL;
854  pszProjection = dataset->GetProjectionRef();
855 
856  if (OSRImportFromWkt(pSR, (char **) (&pszProjection)) == OGRERR_NONE)
857  {
858  char * pszPrettyWkt = NULL;
859  OSRExportToPrettyWkt(pSR, &pszPrettyWkt, FALSE);
860 
861  itk::EncapsulateMetaData<std::string> (dict, MetaDataKey::ProjectionRefKey,
862  static_cast<std::string>(pszPrettyWkt));
863 
864  CPLFree(pszPrettyWkt);
865  }
866  else
867  {
868  itk::EncapsulateMetaData<std::string>(dict, MetaDataKey::ProjectionRefKey,
869  static_cast<std::string>(dataset->GetProjectionRef()));
870  }
871 
872  if (pSR != NULL)
873  {
874  OSRRelease(pSR);
875  pSR = NULL;
876  }
877  }
878 
879  /* -------------------------------------------------------------------- */
880  /* Get the GCP projection coordinates of the image : GCPProjection */
881  /* -------------------------------------------------------------------- */
882 
883  unsigned int gcpCount = 0;
884  gcpCount = dataset->GetGCPCount();
885  if (gcpCount > 0)
886  {
887  std::string gcpProjectionKey = static_cast<std::string>(dataset->GetGCPProjection());
888  itk::EncapsulateMetaData<std::string>(dict, MetaDataKey::GCPProjectionKey, gcpProjectionKey);
889 
890  if (gcpProjectionKey.empty())
891  {
892  gcpCount = 0; //fix for uninitialized gcpCount in gdal (when
893  //reading Palsar image)
894  }
895 
896  std::string key;
897 
898  itk::EncapsulateMetaData<unsigned int>(dict, MetaDataKey::GCPCountKey, gcpCount);
899 
900  for (unsigned int cpt = 0; cpt < gcpCount; ++cpt)
901  {
902 
903  const GDAL_GCP *psGCP;
904  psGCP = dataset->GetGCPs() + cpt;
905 
906  OTB_GCP pOtbGCP;
907  pOtbGCP.m_Id = std::string(psGCP->pszId);
908  pOtbGCP.m_Info = std::string(psGCP->pszInfo);
909  pOtbGCP.m_GCPRow = psGCP->dfGCPLine;
910  pOtbGCP.m_GCPCol = psGCP->dfGCPPixel;
911  pOtbGCP.m_GCPX = psGCP->dfGCPX;
912  pOtbGCP.m_GCPY = psGCP->dfGCPY;
913  pOtbGCP.m_GCPZ = psGCP->dfGCPZ;
914 
915  // Complete the key with the GCP number : GCP_i
916  std::ostringstream lStream;
917  lStream << MetaDataKey::GCPParametersKey << cpt;
918  key = lStream.str();
919 
920  itk::EncapsulateMetaData<OTB_GCP>(dict, key, pOtbGCP);
921 
922  }
923 
924  }
925 
926  /* -------------------------------------------------------------------- */
927  /* Get the six coefficients of affine geoTtransform */
928  /* -------------------------------------------------------------------- */
929 
930  double adfGeoTransform[6];
931  MetaDataKey::VectorType VadfGeoTransform;
932 
933  if (dataset->GetGeoTransform(adfGeoTransform) == CE_None)
934  {
935  for (int cpt = 0; cpt < 6; ++cpt)
936  VadfGeoTransform.push_back(adfGeoTransform[cpt]);
937 
938  itk::EncapsulateMetaData<MetaDataKey::VectorType>(dict, MetaDataKey::GeoTransformKey, VadfGeoTransform);
939 
941  m_Origin[0] = VadfGeoTransform[0];
942  m_Origin[1] = VadfGeoTransform[3];
943  m_Spacing[0] = VadfGeoTransform[1];
944  m_Spacing[1] = VadfGeoTransform[5];
945 
946  }
947 
948  // Dataset info
949  otbMsgDevMacro(<< "**** ReadImageInformation() DATASET INFO: ****" );
950  otbMsgDevMacro(<< "Projection Ref: "<< dataset->GetProjectionRef() );
951  double GT[6];
952  if (dataset->GetGeoTransform(GT) == CE_None)
953  {
954  otbMsgDevMacro( <<"Geo Transform: "<< GT[0] << ", " << GT[1] << ", "
955  << GT[2] << ", " << GT[3] << ", "
956  << GT[4] << ", " << GT[5] );
957  }
958  else
959  {
960  otbMsgDevMacro( << "No Geo Transform: ");
961  }
962  otbMsgDevMacro(<< "GCP Projection Ref: "<< dataset->GetGCPProjection() );
963  otbMsgDevMacro(<< "GCP Count: " << dataset->GetGCPCount() );
964 
965  /* -------------------------------------------------------------------- */
966  /* Report metadata. */
967  /* -------------------------------------------------------------------- */
968 
969  papszMetadata = dataset->GetMetadata(NULL);
970  if (CSLCount(papszMetadata) > 0)
971  {
972  std::string key;
973 
974  for (int cpt = 0; papszMetadata[cpt] != NULL; ++cpt)
975  {
976  std::ostringstream lStream;
977  lStream << MetaDataKey::MetadataKey << cpt;
978  key = lStream.str();
979 
980  itk::EncapsulateMetaData<std::string>(dict, key,
981  static_cast<std::string>(papszMetadata[cpt]));
982  }
983  }
984 
985  /* -------------------------------------------------------------------- */
986  /* Report subdatasets. */
987  /* -------------------------------------------------------------------- */
988 
989  papszMetadata = dataset->GetMetadata("SUBDATASETS");
990  if (CSLCount(papszMetadata) > 0)
991  {
992  std::string key;
993 
994  for (int cpt = 0; papszMetadata[cpt] != NULL; ++cpt)
995  {
996  std::ostringstream lStream;
997  lStream << MetaDataKey::SubMetadataKey << cpt;
998  key = lStream.str();
999 
1000  itk::EncapsulateMetaData<std::string>(dict, key,
1001  static_cast<std::string>(papszMetadata[cpt]));
1002  }
1003  }
1004 
1005  /* -------------------------------------------------------------------- */
1006  /* Report corners */
1007  /* -------------------------------------------------------------------- */
1008 
1009  double GeoX(0), GeoY(0);
1011 
1012  GDALInfoReportCorner("Upper Left", 0.0, 0.0, GeoX, GeoY);
1013  VGeo.push_back(GeoX);
1014  VGeo.push_back(GeoY);
1015 
1016  itk::EncapsulateMetaData<MetaDataKey::VectorType>(dict, MetaDataKey::UpperLeftCornerKey, VGeo);
1017 
1018  VGeo.clear();
1019 
1020  GDALInfoReportCorner("Upper Right", m_Dimensions[0], 0.0, GeoX, GeoY);
1021  VGeo.push_back(GeoX);
1022  VGeo.push_back(GeoY);
1023 
1024  itk::EncapsulateMetaData<MetaDataKey::VectorType>(dict, MetaDataKey::UpperRightCornerKey, VGeo);
1025 
1026  VGeo.clear();
1027 
1028  GDALInfoReportCorner("Lower Left", 0.0, m_Dimensions[1], GeoX, GeoY);
1029  VGeo.push_back(GeoX);
1030  VGeo.push_back(GeoY);
1031 
1032  itk::EncapsulateMetaData<MetaDataKey::VectorType>(dict, MetaDataKey::LowerLeftCornerKey, VGeo);
1033 
1034  VGeo.clear();
1035 
1036  GDALInfoReportCorner("Lower Right", m_Dimensions[0], m_Dimensions[1], GeoX, GeoY);
1037  VGeo.push_back(GeoX);
1038  VGeo.push_back(GeoY);
1039 
1040  itk::EncapsulateMetaData<MetaDataKey::VectorType>(dict, MetaDataKey::LowerRightCornerKey, VGeo);
1041 
1042  VGeo.clear();
1043 
1044  /* -------------------------------------------------------------------- */
1045  /* Color Table */
1046  /* -------------------------------------------------------------------- */
1047 
1048  for (int iBand = 0; iBand < dataset->GetRasterCount(); iBand++)
1049  {
1050  GDALColorTableH hTable;
1051  GDALRasterBandH hBand;
1052  hBand = GDALGetRasterBand(dataset, iBand + 1);
1053  if ((GDALGetRasterColorInterpretation(hBand) == GCI_PaletteIndex)
1054  && (hTable = GDALGetRasterColorTable(hBand)) != NULL)
1055  {
1056  m_IsIndexed = true;
1057 
1058  unsigned int ColorEntryCount = GDALGetColorEntryCount(hTable);
1059 
1060  itk::EncapsulateMetaData<std::string>(dict, MetaDataKey::ColorTableNameKey,
1061  static_cast<std::string>(GDALGetPaletteInterpretationName(
1062  GDALGetPaletteInterpretation(hTable))));
1063 
1064  itk::EncapsulateMetaData<unsigned int>(dict, MetaDataKey::ColorEntryCountKey, ColorEntryCount);
1065 
1066  for (int i = 0; i < GDALGetColorEntryCount(hTable); ++i)
1067  {
1068  GDALColorEntry sEntry;
1069  MetaDataKey::VectorType VColorEntry;
1070 
1071  GDALGetColorEntryAsRGB(hTable, i, &sEntry);
1072 
1073  VColorEntry.push_back(sEntry.c1);
1074  VColorEntry.push_back(sEntry.c2);
1075  VColorEntry.push_back(sEntry.c3);
1076  VColorEntry.push_back(sEntry.c4);
1077 
1078  itk::EncapsulateMetaData<MetaDataKey::VectorType>(dict, MetaDataKey::ColorEntryAsRGBKey, VColorEntry);
1079 
1080  }
1081  }
1082  }
1083  if (m_IsIndexed)
1084  {
1085  m_NbBands *= 4;
1087  this->SetPixelType(VECTOR);
1088  }
1089 
1090  // Adapt the spacing to the resolution read
1091  m_Spacing[0] *= vcl_pow(2.0, static_cast<double>(m_ResolutionFactor));
1092  m_Spacing[1] *= vcl_pow(2.0, static_cast<double>(m_ResolutionFactor));
1093 }
1094 
1095 bool GDALImageIO::CanWriteFile(const char* name)
1096 {
1097  m_FileName = name;
1098 
1099  // First check the filename
1100  if (name == NULL)
1101  {
1102  itkDebugMacro(<< "No filename specified.");
1103  return false;
1104  }
1105 
1106  // Get the GDAL format ID from the name
1107  std::string gdalDriverShortName = FilenameToGdalDriverShortName(name);
1108  if (gdalDriverShortName == "NOT-FOUND")
1109  {
1110  return false;
1111  }
1112 
1113  // Check the driver for support of Create or at least CreateCopy
1114  GDALDriver* driver = GDALDriverManagerWrapper::GetInstance().GetDriverByName(gdalDriverShortName);
1115  if ( GDALGetMetadataItem( driver, GDAL_DCAP_CREATE, NULL ) == NULL
1116  && GDALGetMetadataItem( driver, GDAL_DCAP_CREATECOPY, NULL ) == NULL )
1117  {
1118  itkDebugMacro(<< "The driver " << GDALGetDriverShortName(driver) << " does not support writing");
1119  return false;
1120  }
1121  return true;
1122 }
1123 
1125 {
1126  // Get the GDAL format ID from the name
1127  std::string gdalDriverShortName = FilenameToGdalDriverShortName(m_FileName);
1128  GDALDriver* driver = GDALDriverManagerWrapper::GetInstance().GetDriverByName(gdalDriverShortName);
1129 
1130  if (driver == NULL)
1131  {
1132  itkDebugMacro(<< "Unable to instantiate driver " << gdalDriverShortName);
1133  m_CanStreamWrite = false;
1134  }
1135  if ( GDALGetMetadataItem( driver, GDAL_DCAP_CREATE, NULL ) != NULL )
1136  {
1137  m_CanStreamWrite = true;
1138  }
1139  else
1140  {
1141  m_CanStreamWrite = false;
1142  }
1143 
1144  return m_CanStreamWrite;
1145 }
1146 
1147 void GDALImageIO::Write(const void* buffer)
1148 {
1149  // Check if we have to write the image information
1150  if (m_FlagWriteImageInformation == true)
1151  {
1152  this->InternalWriteImageInformation(buffer);
1154  }
1155 
1156  // Check if conversion succeed
1157  if (buffer == NULL)
1158  {
1159  itkExceptionMacro(<< "GDAL : Bad alloc");
1160  return;
1161  }
1162 
1163  // Compute offset and size
1164  unsigned int lNbLines = this->GetIORegion().GetSize()[1];
1165  unsigned int lNbColumns = this->GetIORegion().GetSize()[0];
1166  int lFirstLine = this->GetIORegion().GetIndex()[1]; // [1... ]
1167  int lFirstColumn = this->GetIORegion().GetIndex()[0]; // [1... ]
1168 
1169  // Particular case: checking that the written region is the same size
1170  // of the entire image
1171  // starting at offset 0 (when no streaming)
1172  if ((lNbLines == m_Dimensions[1]) && (lNbColumns == m_Dimensions[0]))
1173  {
1174  lFirstLine = 0;
1175  lFirstColumn = 0;
1176  }
1177 
1178  // Convert buffer from void * to unsigned char *
1179  //unsigned char *p = static_cast<unsigned char*>( const_cast<void *>(buffer));
1180  //printDataBuffer(p, m_PxType->pixType, m_NbBands, 10*2); // Buffer incorrect
1181 
1182  // If driver supports streaming
1183  if (m_CanStreamWrite)
1184  {
1185 
1186  otbMsgDevMacro(<< "RasterIO Write requested region : " << this->GetIORegion() <<
1187  "\n, lFirstColumn =" << lFirstColumn <<
1188  "\n, lFirstLine =" << lFirstLine <<
1189  "\n, lNbColumns =" << lNbColumns <<
1190  "\n, lNbLines =" << lNbLines <<
1191  "\n, m_PxType =" << GDALGetDataTypeName(m_PxType->pixType) <<
1192  "\n, m_NbBands =" << m_NbBands <<
1193  "\n, m_BytePerPixel ="<< m_BytePerPixel <<
1194  "\n, Pixel offset =" << m_BytePerPixel * m_NbBands << // is nbComp * BytePerPixel
1195  "\n, Line offset =" << m_BytePerPixel * m_NbBands * lNbColumns << // is pixelOffset * nbColumns
1196  "\n, Band offset =" << m_BytePerPixel) // is BytePerPixel
1197 
1198  itk::TimeProbe chrono;
1199  chrono.Start();
1200  CPLErr lCrGdal = m_Dataset->GetDataSet()->RasterIO(GF_Write,
1201  lFirstColumn,
1202  lFirstLine,
1203  lNbColumns,
1204  lNbLines,
1205  const_cast<void *>(buffer),
1206  lNbColumns,
1207  lNbLines,
1208  m_PxType->pixType,
1209  m_NbBands,
1210  // We want to write all bands
1211  NULL,
1212  // Pixel offset
1213  // is nbComp * BytePerPixel
1215  // Line offset
1216  // is pixelOffset * nbColumns
1217  m_BytePerPixel * m_NbBands * lNbColumns,
1218  // Band offset is BytePerPixel
1219  m_BytePerPixel);
1220  chrono.Stop();
1221  otbMsgDevMacro(<< "RasterIO Write took " << chrono.GetTotal() << " sec")
1222 
1223  // Check if writing succeed
1224  if (lCrGdal == CE_Failure)
1225  {
1226  itkExceptionMacro(<< "Error while writing image (GDAL format) " << m_FileName.c_str() << ".");
1227  }
1228  // Flush dataset cache
1229  m_Dataset->GetDataSet()->FlushCache();
1230  }
1231  else
1232  {
1233  // We only wrote data to the memory dataset
1234  // Now write it to the real file with CreateCopy()
1235  std::string gdalDriverShortName = FilenameToGdalDriverShortName(m_FileName);
1236  std::string realFileName = GetGdalWriteImageFileName(gdalDriverShortName, m_FileName);
1237 
1238  GDALDriver* driver = GDALDriverManagerWrapper::GetInstance().GetDriverByName(gdalDriverShortName);
1239  if (driver == NULL)
1240  {
1241  itkExceptionMacro(<< "Unable to instantiate driver " << gdalDriverShortName << " to write " << m_FileName);
1242  }
1243 
1244  GDALCreationOptionsType creationOptions = m_CreationOptions;
1245  GDALDataset* hOutputDS = driver->CreateCopy( realFileName.c_str(), m_Dataset->GetDataSet(), FALSE,
1246  otb::ogr::StringListConverter(creationOptions).to_ogr(),
1247  NULL, NULL );
1248  GDALClose(hOutputDS);
1249  }
1250 
1251 
1252  if (lFirstLine + lNbLines == m_Dimensions[1]
1253  && lFirstColumn + lNbColumns == m_Dimensions[0])
1254  {
1255  // Last pixel written
1256  // Reinitialize to close the file
1258  }
1259 }
1260 
1263 {
1264 }
1265 
1267 {
1268  //char ** papszOptions = NULL;
1269  std::string driverShortName;
1270  m_NbBands = this->GetNumberOfComponents();
1271 
1272  if ((m_Dimensions[0] == 0) && (m_Dimensions[1] == 0))
1273  {
1274  itkExceptionMacro(<< "Dimensions are not defined.");
1275  }
1276 
1277  if ((this->GetPixelType() == COMPLEX) /*&& (m_NbBands / 2 > 0)*/)
1278  {
1279  //m_NbBands /= 2;
1280 
1281  if (this->GetComponentType() == CSHORT)
1282  {
1283  m_BytePerPixel = 4;
1284  m_PxType->pixType = GDT_CInt16;
1285  }
1286  else if (this->GetComponentType() == CINT)
1287  {
1288  m_BytePerPixel = 8;
1289  m_PxType->pixType = GDT_CInt32;
1290  }
1291  else if (this->GetComponentType() == CFLOAT)
1292  {
1293  m_BytePerPixel = 8;
1294  m_PxType->pixType = GDT_CFloat32;
1295  }
1296  else if (this->GetComponentType() == CDOUBLE)
1297  {
1298  m_BytePerPixel = 16;
1299  m_PxType->pixType = GDT_CFloat64;
1300  }
1301  else
1302  {
1303  itkExceptionMacro(<< "This complex type is not defined :" << this->GetPixelTypeAsString(this->GetPixelType()) );
1304  }
1305  }
1306  else
1307  {
1308  if (this->GetComponentType() == CHAR)
1309  {
1310  m_BytePerPixel = 1;
1311  m_PxType->pixType = GDT_Byte;
1312  }
1313  else if (this->GetComponentType() == UCHAR)
1314  {
1315  m_BytePerPixel = 1;
1316  m_PxType->pixType = GDT_Byte;
1317  }
1318  else if (this->GetComponentType() == USHORT)
1319  {
1320  m_BytePerPixel = 2;
1321  m_PxType->pixType = GDT_UInt16;
1322  }
1323  else if (this->GetComponentType() == SHORT)
1324  {
1325  m_BytePerPixel = 2;
1326  m_PxType->pixType = GDT_Int16;
1327  }
1328  else if (this->GetComponentType() == INT)
1329  {
1330  m_BytePerPixel = 4;
1331  m_PxType->pixType = GDT_Int32;
1332  }
1333  else if (this->GetComponentType() == UINT)
1334  {
1335  m_BytePerPixel = 4;
1336  m_PxType->pixType = GDT_UInt32;
1337  }
1338  else if (this->GetComponentType() == LONG)
1339  {
1340  m_BytePerPixel = sizeof(long);
1341  if( m_BytePerPixel == 8 )
1342  {
1343  itkWarningMacro(<< "Cast a long (64 bits) image into an int (32 bits) one.")
1344  }
1345  m_PxType->pixType = GDT_Int32;
1346  }
1347  else if (this->GetComponentType() == ULONG)
1348  {
1349  m_BytePerPixel = sizeof(unsigned long);
1350  if( m_BytePerPixel == 8 )
1351  {
1352  itkWarningMacro(<< "Cast an unsigned long (64 bits) image into an unsigned int (32 bits) one.")
1353  }
1354  m_PxType->pixType = GDT_UInt32;
1355  }
1356  else if (this->GetComponentType() == FLOAT)
1357  {
1358  m_BytePerPixel = 4;
1359  m_PxType->pixType = GDT_Float32;
1360  }
1361  else if (this->GetComponentType() == DOUBLE)
1362  {
1363  m_BytePerPixel = 8;
1364  m_PxType->pixType = GDT_Float64;
1365  }
1366  else
1367  {
1368  m_BytePerPixel = 1;
1369  m_PxType->pixType = GDT_Byte;
1370  }
1371  }
1372 
1373  // Automatically set the Type to Binary for GDAL data
1374  this->SetFileTypeToBinary();
1375 
1376  driverShortName = FilenameToGdalDriverShortName(m_FileName);
1377  if (driverShortName == "NOT-FOUND")
1378  {
1379  itkExceptionMacro(
1380  << "GDAL Writing failed : the image file name '" << m_FileName.c_str() << "' is not recognized by GDAL.");
1381  }
1382 
1383  if (m_CanStreamWrite)
1384  {
1385  GDALCreationOptionsType creationOptions = m_CreationOptions;
1386 /*
1387  // Force tile mode for TIFF format if no creation option are given
1388  if( driverShortName == "GTiff" )
1389  {
1390  if ( CreationOptionContains( "TILED=YES" ) )
1391  {
1392  // User requested tiled TIFF explicitly
1393  //
1394  // Let GDAL set up the BLOCKXSIZE and BLOCKYSIZE
1395  // or suppose the user have set it also along with TILED=YES
1396  // This allows the user to have complete
1397  // control over the tiling scheme
1398  }
1399  else if ( CreationOptionContains( "BLOCKYSIZE=" ) )
1400  {
1401  // User did not set "TILED=YES" but set "BLOCKYSIZE="
1402  // -> He requested a stripped TIFF
1403  }
1404  else
1405  {
1406  // User did not specify "TILED=YES" nor "BLOCKYSIZE=?"
1407  // Switch to TILED mode automatically, and choose BLOCKXSIZE and BLOCKYSIZE for him
1408 
1409  otbMsgDevMacro(<< "Enabling TIFF Tiled mode")
1410 
1411  // Use a fixed tile size
1412  // Take as reference is a 256*256 short int 4 bands tile
1413  const unsigned int ReferenceTileSizeInBytes = 256 * 256 * 4 * 2;
1414  const unsigned int NbPixelPerTile = ReferenceTileSizeInBytes / m_BytePerPixel / m_NbBands;
1415  const unsigned int IdealTileDimension = static_cast<unsigned int>( vcl_sqrt(static_cast<float>(NbPixelPerTile)) );
1416 
1417  // Set tileDimension to the nearest power of two and aligned to
1418  // 16 pixels (needed by TIFF spec)
1419  unsigned int tileDimension = 16;
1420  while(2*tileDimension < IdealTileDimension)
1421  {
1422  tileDimension*=2;
1423  }
1424  otbMsgDevMacro(<< "Tile dimension : " << tileDimension << " * " << tileDimension)
1425 
1426  std::ostringstream tileDimensionStr;
1427  tileDimensionStr << tileDimension;
1428 
1429  creationOptions.push_back( "TILED=YES" );
1430  creationOptions.push_back( std::string("BLOCKXSIZE=") + tileDimensionStr.str() );
1431  creationOptions.push_back( std::string("BLOCKYSIZE=") + tileDimensionStr.str() );
1432  }
1433  }
1434 */
1436  driverShortName,
1437  GetGdalWriteImageFileName(driverShortName, m_FileName),
1438  m_Dimensions[0], m_Dimensions[1],
1440  otb::ogr::StringListConverter(creationOptions).to_ogr());
1441  }
1442  else
1443  {
1444  // buffer casted in unsigned long cause under Win32 the address
1445  // doesn't begin with 0x, the address in not interpreted as
1446  // hexadecimal but alpha numeric value, then the conversion to
1447  // integer make us pointing to an non allowed memory block => Crash.
1448  std::ostringstream stream;
1449  stream << "MEM:::"
1450  << "DATAPOINTER=" << (unsigned long)(buffer) << ","
1451  << "PIXELS=" << m_Dimensions[0] << ","
1452  << "LINES=" << m_Dimensions[1] << ","
1453  << "BANDS=" << m_NbBands << ","
1454  << "DATATYPE=" << GDALGetDataTypeName(m_PxType->pixType) << ","
1455  << "PIXELOFFSET=" << m_BytePerPixel * m_NbBands << ","
1456  << "LINEOFFSET=" << m_BytePerPixel * m_NbBands * m_Dimensions[0] << ","
1457  << "BANDOFFSET=" << m_BytePerPixel;
1458 
1460  }
1461 
1462  if (m_Dataset.IsNull())
1463  {
1464  itkExceptionMacro(
1465  << "GDAL Writing failed : Impossible to create the image file name '" << m_FileName << "'.");
1466  }
1467 
1468  /*----------------------------------------------------------------------*/
1469  /*-------------------------- METADATA ----------------------------------*/
1470  /*----------------------------------------------------------------------*/
1471 
1472  // Now initialize the itk dictionary
1474  std::ostringstream oss;
1475  GDALDataset* dataset = m_Dataset->GetDataSet();
1476 
1477  std::string projectionRef;
1478  itk::ExposeMetaData<std::string>(dict, MetaDataKey::ProjectionRefKey, projectionRef);
1479 
1480  /* -------------------------------------------------------------------- */
1481  /* Set the GCPs */
1482  /* -------------------------------------------------------------------- */
1483  const double Epsilon = 1E-10;
1484  if (projectionRef.empty()
1485  && (vcl_abs(m_Origin[0]) > Epsilon
1486  || vcl_abs(m_Origin[1]) > Epsilon
1487  || vcl_abs(m_Spacing[0] - 1) > Epsilon
1488  || vcl_abs(m_Spacing[1] - 1) > Epsilon) )
1489  {
1490  // See issue #303 :
1491  // If there is no ProjectionRef, and the GeoTransform is not the identity,
1492  // then saving also GCPs is undefined behavior for GDAL, and a WGS84 projection crs
1493  // is assigned arbitrarily
1494  otbMsgDevMacro(<< "Skipping GCPs saving to prevent GDAL from assigning a WGS84 projection ref to the file")
1495  }
1496  else
1497  {
1498  unsigned int gcpCount = 0;
1499  itk::ExposeMetaData<unsigned int>(dict, MetaDataKey::GCPCountKey, gcpCount);
1500 
1501  if (gcpCount > 0)
1502  {
1503 
1504  GDAL_GCP * gdalGcps = new GDAL_GCP[gcpCount];
1505 
1506  for (unsigned int gcpIndex = 0; gcpIndex < gcpCount; ++gcpIndex)
1507  {
1508  //Build the GCP string in the form of GCP_n
1509  std::ostringstream lStream;
1510  lStream << MetaDataKey::GCPParametersKey << gcpIndex;
1511  std::string key = lStream.str();
1512 
1513  OTB_GCP gcp;
1514  itk::ExposeMetaData<OTB_GCP>(dict, key, gcp);
1515 
1516  gdalGcps[gcpIndex].pszId = const_cast<char *>(gcp.m_Id.c_str());
1517  gdalGcps[gcpIndex].pszInfo = const_cast<char *>(gcp.m_Info.c_str());
1518  gdalGcps[gcpIndex].dfGCPPixel = gcp.m_GCPCol;
1519  gdalGcps[gcpIndex].dfGCPLine = gcp.m_GCPRow;
1520  gdalGcps[gcpIndex].dfGCPX = gcp.m_GCPX;
1521  gdalGcps[gcpIndex].dfGCPY = gcp.m_GCPY;
1522  gdalGcps[gcpIndex].dfGCPZ = gcp.m_GCPZ;
1523  }
1524 
1525  std::string gcpProjectionRef;
1526  itk::ExposeMetaData<std::string>(dict, MetaDataKey::GCPProjectionKey, gcpProjectionRef);
1527 
1528  dataset->SetGCPs(gcpCount, gdalGcps, gcpProjectionRef.c_str());
1529 
1530  delete[] gdalGcps;
1531  }
1532  }
1533 
1534  /* -------------------------------------------------------------------- */
1535  /* Set the projection coordinate system of the image : ProjectionRef */
1536  /* -------------------------------------------------------------------- */
1537  if (!projectionRef.empty())
1538  {
1539  dataset->SetProjection(projectionRef.c_str());
1540  }
1541 
1542  /* -------------------------------------------------------------------- */
1543  /* Set the six coefficients of affine geoTransform */
1544  /* -------------------------------------------------------------------- */
1545  itk::VariableLengthVector<double> geoTransform(6);
1547  geoTransform[0] = m_Origin[0];
1548  geoTransform[3] = m_Origin[1];
1549  geoTransform[1] = m_Spacing[0];
1550  geoTransform[5] = m_Spacing[1];
1551 
1552  // FIXME: Here component 1 and 4 should be replaced by the orientation parameters
1553  geoTransform[2] = 0.;
1554  geoTransform[4] = 0.;
1555  dataset->SetGeoTransform(const_cast<double*>(geoTransform.GetDataPointer()));
1556 
1557  /* -------------------------------------------------------------------- */
1558  /* Report metadata. */
1559  /* -------------------------------------------------------------------- */
1560 
1561  std::string svalue = "";
1562  std::vector<std::string> keys = dict.GetKeys();
1563  std::string const metadataKey = MetaDataKey::MetadataKey;
1564 
1565  for (unsigned int itkey = 0; itkey < keys.size(); ++itkey)
1566  {
1568  if (keys[itkey].compare(0, metadataKey.length(), metadataKey) == 0)
1569  {
1570  itk::ExposeMetaData<std::string>(dict, keys[itkey], svalue);
1571  unsigned int equalityPos = svalue.find_first_of('=');
1572  std::string tag = svalue.substr(0, equalityPos);
1573  std::string value = svalue.substr(equalityPos + 1);
1574  otbMsgDevMacro(<< "Metadata: " << tag << "=" << value);
1575  dataset->SetMetadataItem(tag.c_str(), value.c_str(), NULL);
1576  }
1577  }
1578  // END
1579 
1580  // Dataset info
1581  otbMsgDevMacro( << "**** WriteImageInformation() DATASET INFO: ****" );
1582  otbMsgDevMacro( << "Projection Ref: "<<dataset->GetProjectionRef() );
1583  double GT[6];
1584  if (dataset->GetGeoTransform(GT) == CE_None)
1585  {
1586  otbMsgDevMacro( <<"Geo Transform: "<< GT[0] << ", " << GT[1] << ", "
1587  << GT[2] << ", " << GT[3] << ", "
1588  << GT[4] << ", " << GT[5] );
1589  }
1590  else
1591  {
1592  otbMsgDevMacro( << "No Geo Transform: ");
1593  }
1594 
1595  otbMsgDevMacro( << "GCP Projection Ref: "<< dataset->GetGCPProjection() );
1596  otbMsgDevMacro( << "GCP Count: " << dataset->GetGCPCount() );
1597 
1598 }
1599 
1600 std::string GDALImageIO::FilenameToGdalDriverShortName(const std::string& name) const
1601 {
1602  std::string extension;
1603  std::string gdalDriverShortName;
1604 
1605  // Get extension in lowercase
1606  extension = itksys::SystemTools::LowerCase( itksys::SystemTools::GetFilenameLastExtension(name) );
1607 
1608  if ( extension == ".tif" || extension == ".tiff" )
1609  gdalDriverShortName = "GTiff";
1610  else if ( extension == ".hdr" )
1611  gdalDriverShortName = "ENVI";
1612  else if ( extension == ".img" )
1613  gdalDriverShortName = "HFA";
1614  else if ( extension == ".ntf" )
1615  gdalDriverShortName = "NITF";
1616  else if ( extension == ".png" )
1617  gdalDriverShortName="PNG";
1618  else if ( extension == ".jpg" || extension== ".jpeg" )
1619  gdalDriverShortName="JPEG";
1620  else if ( extension == ".pix" )
1621  gdalDriverShortName="PCIDSK";
1622  else if ( extension == ".lbl" || extension == ".pds" )
1623  gdalDriverShortName="ISIS2";
1624  else
1625  gdalDriverShortName = "NOT-FOUND";
1626 
1627  return gdalDriverShortName;
1628 }
1629 
1630 std::string GDALImageIO::GetGdalWriteImageFileName(const std::string& gdalDriverShortName, const std::string& filename) const
1631 {
1632  std::string gdalFileName;
1633 
1634  gdalFileName = filename;
1635  // Suppress hdr extension for ENVI format
1636  if (gdalDriverShortName == "ENVI")
1637  {
1638  gdalFileName = System::GetRootName(filename);
1639  }
1640  return gdalFileName;
1641 }
1642 
1643 bool GDALImageIO::GDALInfoReportCorner(const char * /*corner_name*/, double x, double y, double& GeoX, double& GeoY) const
1644 {
1645  double adfGeoTransform[6];
1646  bool IsTrue;
1647 
1648  /* -------------------------------------------------------------------- */
1649  /* Transform the point into georeferenced coordinates. */
1650  /* -------------------------------------------------------------------- */
1651  if (m_Dataset->GetDataSet()->GetGeoTransform(adfGeoTransform) == CE_None)
1652  {
1653  GeoX = adfGeoTransform[0] + adfGeoTransform[1] * x + adfGeoTransform[2] * y;
1654  GeoY = adfGeoTransform[3] + adfGeoTransform[4] * x + adfGeoTransform[5] * y;
1655  IsTrue = true;
1656  }
1657 
1658  else
1659  {
1660  GeoX = x;
1661  GeoY = y;
1662  IsTrue = false;
1663  }
1664 
1665  return IsTrue;
1666 }
1667 
1668 bool GDALImageIO::CreationOptionContains(std::string partialOption) const
1669 {
1670  size_t i;
1671  for (i = 0; i < m_CreationOptions.size(); ++i)
1672  {
1673  if (boost::algorithm::starts_with(m_CreationOptions[i], partialOption))
1674  {
1675  break;
1676  }
1677  }
1678  return (i != m_CreationOptions.size());
1679 }
1680 
1681 } // end namespace otb

Generated at Sat Mar 8 2014 15:56:47 for Orfeo Toolbox with doxygen 1.8.3.1