Rebase vtkGDALRasterReader to Paraview 5.1.0 branch

IO/GDAL/Testing/Cxx/CMakeLists.txt

+set(TestGDALRasterNoDataValue_ARGS
+  -D DATA{../Data/Input/TestGDALRasterNoDataValue.tif}
+)
+set(TestGDALRasterPalette_ARGS
+  -D DATA{../Data/Input/TestGDALRasterPalette.tif}
+)
+
 vtk_add_test_cxx(${vtk-module}CxxTests tests
   TestGDALVectorReader.cxx
   TestGDALRasterReader.cxx
+  TestGDALRasterNoDataValue.cxx,NO_VALID,NO_OUTPUT
+  TestGDALRasterPalette.cxx
   )
 vtk_test_cxx_executable(${vtk-module}CxxTests tests)

IO/GDAL/Testing/Cxx/TestGDALRasterNoDataValue.cxx

+/*=========================================================================
+
+  Program:   Visualization Toolkit
+  Module:    TestGDALRasterReader.cxx
+
+  Copyright (c) Ken Martin, Will Schroeder, Bill Lorensen
+  All rights reserved.
+  See Copyright.txt or http://www.kitware.com/Copyright.htm for details.
+
+     This software is distributed WITHOUT ANY WARRANTY; without even
+     the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
+     PURPOSE.  See the above copyright notice for more information.
+
+=========================================================================*/
+#include <vtkGDALRasterReader.h>
+#include <vtkNew.h>
+#include <vtkUniformGrid.h>
+
+#include <iostream>
+
+// Main program
+int TestGDALRasterNoDataValue(int argc, char** argv)
+{
+  if (argc < 3)
+    {
+    std::cerr << "Expected TestName -D InputFile.tif" << std::endl;
+    return -1;
+    }
+
+  std::string inputFileName(argv[2]);
+
+  // Create reader to read shape file.
+  vtkNew<vtkGDALRasterReader> reader;
+  reader->SetFileName(inputFileName.c_str());
+  reader->Update();
+
+  vtkUniformGrid *rasterImage = vtkUniformGrid::SafeDownCast(
+    reader->GetOutput());
+
+  int numErrors = 0;
+
+  if (!rasterImage->HasAnyBlankPoints())
+    {
+    std::cerr << "Error image has no blank points" << std::endl;
+    ++numErrors;
+    }
+
+  if (!rasterImage->HasAnyBlankCells())
+    {
+    std::cerr << "Error image has no blank cells" << std::endl;
+    ++numErrors;
+    }
+
+  double *scalarRange = rasterImage->GetScalarRange();
+
+  if ((scalarRange[0] < -888.5) || (scalarRange[0]) > -887.5)
+    {
+    std::cerr << "Error scalarRange[0] should be -888.0, not "
+              << scalarRange[0] << std::endl;
+    ++numErrors;
+    }
+
+  if ((scalarRange[1] < 9998.5) || (scalarRange[1] > 9999.5))
+    {
+    std::cerr << "Error scalarRange[1] should be 9999.0, not "
+              << scalarRange[1] << std::endl;
+    ++numErrors;
+    }
+
+  //std::cout << "numErrors: " << numErrors << std::endl;
+  return numErrors;
+}

IO/GDAL/Testing/Cxx/TestGDALRasterPalette.cxx

+/*=========================================================================
+
+  Program:   Visualization Toolkit
+  Module:    TestGDALRasterReader.cxx
+
+  Copyright (c) Ken Martin, Will Schroeder, Bill Lorensen
+  All rights reserved.
+  See Copyright.txt or http://www.kitware.com/Copyright.htm for details.
+
+     This software is distributed WITHOUT ANY WARRANTY; without even
+     the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
+     PURPOSE.  See the above copyright notice for more information.
+
+=========================================================================*/
+#include <vtkDataArray.h>
+#include <vtkGDALRasterReader.h>
+#include <vtkImageActor.h>
+#include <vtkImageProperty.h>
+#include <vtkLookupTable.h>
+#include <vtkNew.h>
+#include <vtkPointData.h>
+#include <vtkRegressionTestImage.h>
+#include <vtkRenderWindow.h>
+#include <vtkRenderer.h>
+#include <vtkRenderWindow.h>
+#include <vtkRenderWindowInteractor.h>
+#include <vtkUniformGrid.h>
+
+#include <iostream>
+
+// Main program
+int TestGDALRasterPalette(int argc, char** argv)
+{
+  if (argc < 3)
+    {
+    std::cerr << "Expected TestName -D InputFile.tif" << std::endl;
+    return -1;
+    }
+
+  std::string inputFileName(argv[2]);
+
+  // Create reader to read shape file.
+  vtkNew<vtkGDALRasterReader> reader;
+  reader->SetFileName(inputFileName.c_str());
+  reader->Update();
+  vtkUniformGrid *image = vtkUniformGrid::SafeDownCast(reader->GetOutput());
+
+  // Check that reader generated point scalars
+  if (image->GetPointData()->GetNumberOfArrays() < 1)
+    {
+    std::cerr << "ERROR: Missing point data scalars" << std::endl;
+    return 1;
+    }
+  if (image->GetPointData()->GetScalars()->GetSize() == 0)
+    {
+    std::cerr << "ERROR: Point data scalars empty" << std::endl;
+    return 1;
+    }
+  //image->GetPointData()->GetScalars()->Print(std::cout);
+
+  // Check that reader generated color table
+  vtkLookupTable *colorTable =
+    image->GetPointData()->GetScalars()->GetLookupTable();
+  if (!colorTable)
+    {
+    std::cerr << "ERROR: Missing color table" << std::endl;
+    return 1;
+    }
+  if (colorTable->GetNumberOfAvailableColors() != 256)
+    {
+    std::cerr << "ERROR: Color table does not have 256 colors."
+              << " Instead has " <<  colorTable->GetNumberOfAvailableColors()
+              << std::endl;
+    return 1;
+    }
+  //colorTable->Print(std::cout);
+
+  // Create a renderer and actor
+  vtkNew<vtkRenderer> renderer;
+  vtkNew<vtkImageActor> actor;
+  actor->SetInputData(reader->GetOutput());
+  actor->InterpolateOff();
+  //actor->GetProperty()->SetInterpolationTypeToNearest();
+  actor->GetProperty()->SetLookupTable(colorTable);
+  actor->GetProperty()->UseLookupTableScalarRangeOn();
+  renderer->AddActor(actor.GetPointer());
+
+  // Create a render window, and an interactor
+  vtkNew<vtkRenderWindow> renderWindow;
+  vtkNew<vtkRenderWindowInteractor> renderWindowInteractor;
+  renderWindow->AddRenderer(renderer.GetPointer());
+  renderWindowInteractor->SetRenderWindow(renderWindow.GetPointer());
+
+  //Add the actor to the scene
+  renderer->SetBackground(1.0, 1.0, 1.0);
+  renderWindow->SetSize(400, 400);
+  renderWindow->Render();
+  renderer->ResetCamera();
+  renderWindow->Render();
+
+  int retVal = vtkRegressionTestImage(renderWindow.GetPointer());
+  if (retVal == vtkRegressionTester::DO_INTERACTOR)
+    {
+    renderWindowInteractor->Start();
+    }
+
+  return !retVal;
+}

IO/GDAL/Testing/Data/Baseline/TestGDALRasterPalette.png.md5

+70aff17ba9e4a0d5c6cffb720f381d26

IO/GDAL/Testing/Data/Input/TestGDALRasterNoDataValue.tif.md5

+f73cb0b2de8941d9e491190d193f10a2

IO/GDAL/Testing/Data/Input/TestGDALRasterPalette.tif.md5

+cdffbf2b4f05e98cf0c0f10d34eb302a

IO/GDAL/vtkGDALRasterReader.cxx

@@ -24,6 +24,7 @@
 #include "vtkInformationVector.h"
 #include "vtkInformation.h"
 #include "vtkIntArray.h"
+#include "vtkLookupTable.h"
 #include "vtkMath.h"
 #include "vtkObjectFactory.h"
 #include "vtkPointData.h"
@@ -44,6 +45,7 @@
 // C/C++ includes
 #include <cassert>
 #include <iostream>
+#include <sstream>
 #include <vector>
 
 vtkStandardNewMacro(vtkGDALRasterReader);
@@ -84,6 +86,9 @@ public:
 
   const double* GetGeoCornerPoints();
 
+  void ReadColorTable(
+    GDALRasterBand *rasterBand, vtkLookupTable *colorTable) const;
+
   int NumberOfBands;
   int NumberOfBytesPerPixel;
 
@@ -101,6 +106,7 @@ public:
   // Upper left, lower left, upper right, lower right
   double CornerPoints[8];
 
+  int HasNoDataValue;
   double NoDataValue;
   vtkIdType NumberOfPoints;
 
@@ -277,17 +283,23 @@ void vtkGDALRasterReader::vtkGDALRasterReaderInternal::GenericReadData()
   // Pixel data.
   std::vector<RAW_TYPE> rawUniformGridData;
 
+  // Color table
+  vtkSmartPointer<vtkLookupTable> colorTable =
+    vtkSmartPointer<vtkLookupTable>::New();
+
   // Possible bands
   GDALRasterBand* redBand = 0;
   GDALRasterBand* greenBand = 0;
   GDALRasterBand* blueBand = 0;
   GDALRasterBand* alphaBand = 0;
   GDALRasterBand* greyBand = 0;
+  GDALRasterBand* paletteBand = 0;
 
   for (int i = 1; i <= this->NumberOfBands; ++i)
     {
     GDALRasterBand* rasterBand = this->GDALData->GetRasterBand(i);
-    NoDataValue = rasterBand->GetNoDataValue();
+    this->HasNoDataValue = 0;
+    this->NoDataValue = rasterBand->GetNoDataValue(&this->HasNoDataValue);
     if (this->NumberOfBytesPerPixel == 0)
       {
       this->TargetDataType = rasterBand->GetRasterDataType();
@@ -328,6 +340,10 @@ void vtkGDALRasterReader::vtkGDALRasterReaderInternal::GenericReadData()
       {
       greyBand = rasterBand;
       }
+    else if (rasterBand->GetColorInterpretation() == GCI_PaletteIndex)
+      {
+      paletteBand = rasterBand;
+      }
     }
 
   const int& destWidth = this->Reader->TargetDimensions[0];
@@ -434,6 +450,20 @@ void vtkGDALRasterReader::vtkGDALRasterReaderInternal::GenericReadData()
       assert(err == CE_None);
       }
     }
+  else if (paletteBand)
+    {
+    // Read indexes
+    this->Reader->SetNumberOfScalarComponents(1);
+    rawUniformGridData.resize(destWidth * destHeight * pixelSpace);
+    err = paletteBand->RasterIO(
+      GF_Read, windowX, windowY,  windowWidth, windowHeight,
+      static_cast<void*>(reinterpret_cast<GByte*>(&rawUniformGridData[0]) +
+      0 * bandSpace), destWidth, destHeight,
+      this->TargetDataType, pixelSpace, lineSpace);
+    assert(err == CE_None);
+
+    this->ReadColorTable(paletteBand, colorTable.GetPointer());
+    }
   else
     {
     std::cerr << "Unknown raster band type \n";
@@ -451,6 +481,13 @@ void vtkGDALRasterReader::vtkGDALRasterReaderInternal::GenericReadData()
   this->UniformGridData->SetSpacing(geoSpacing[0], geoSpacing[1], geoSpacing[2]);
   this->UniformGridData->SetOrigin(d[0], d[1], 0);
   this->Convert<VTK_TYPE, RAW_TYPE>(rawUniformGridData, destWidth, destHeight);
+
+  if (paletteBand)
+    {
+    this->UniformGridData->GetPointData()->GetScalars()->SetName("Categories");
+    this->UniformGridData->GetPointData()->GetScalars()->SetLookupTable(
+      colorTable);
+    }
 }
 
 //----------------------------------------------------------------------------
@@ -480,6 +517,12 @@ void vtkGDALRasterReader::vtkGDALRasterReaderInternal::Convert(
   scArr->SetNumberOfComponents(this->NumberOfBands);
   scArr->SetNumberOfTuples(targetWidth * targetHeight);
 
+  RAW_TYPE TNoDataValue = 0;
+  if (this->HasNoDataValue)
+    {
+    TNoDataValue = static_cast<RAW_TYPE>(this->NoDataValue);
+    }
+
   for (int j = 0; j < targetHeight; ++j)
     {
     for (int i = 0; i < targetWidth; ++i)
@@ -491,11 +534,18 @@ void vtkGDALRasterReader::vtkGDALRasterReaderInternal::Convert(
                       j * targetWidth * NumberOfBands + bandIndex;
         sourceIndex = i + j * targetWidth +
                       bandIndex * targetWidth * targetHeight;
+
         RAW_TYPE tmp = rawUniformGridData[sourceIndex];
-        if(tmp < min) min = tmp;
-        if(tmp > max) max = tmp;
-        if(tmp == NoDataValue) this->UniformGridData->BlankPoint(targetIndex);
-        else this->NumberOfPoints++;
+        if (this->HasNoDataValue && tmp == TNoDataValue)
+          {
+          this->UniformGridData->BlankPoint(targetIndex);
+          }
+        else
+          {
+          if(tmp < min) min = tmp;
+          if(tmp > max) max = tmp;
+          this->NumberOfPoints++;
+          }
 
         scArr->InsertValue(targetIndex, rawUniformGridData[sourceIndex]);
         }
@@ -595,6 +645,53 @@ const double* vtkGDALRasterReader::vtkGDALRasterReaderInternal::GetGeoCornerPoin
   return this->CornerPoints;
 }
 
+//-----------------------------------------------------------------------------
+void vtkGDALRasterReader::vtkGDALRasterReaderInternal::ReadColorTable(
+  GDALRasterBand *rasterBand, vtkLookupTable *colorTable) const
+{
+  GDALColorTable *gdalTable = rasterBand->GetColorTable();
+  if (gdalTable->GetPaletteInterpretation() != GPI_RGB)
+    {
+    std::cerr << "Color table palette type not supported "
+              << gdalTable->GetPaletteInterpretation() << std::endl;
+    return;
+    }
+
+  char **categoryNames = rasterBand->GetCategoryNames();
+
+  colorTable->IndexedLookupOn();
+  int numEntries = gdalTable->GetColorEntryCount();
+  colorTable->SetNumberOfTableValues(numEntries);
+  std::stringstream ss;
+  for (int i=0; i< numEntries; ++i)
+    {
+    const GDALColorEntry *gdalEntry = gdalTable->GetColorEntry(i);
+    double r = static_cast<double>(gdalEntry->c1) / 255.0;
+    double g = static_cast<double>(gdalEntry->c2) / 255.0;
+    double b = static_cast<double>(gdalEntry->c3) / 255.0;
+    double a = static_cast<double>(gdalEntry->c4) / 255.0;
+    colorTable->SetTableValue(i, r, g, b, a);
+
+    // Copy category name to lookup table annotation
+    if (categoryNames)
+      {
+      // Only use non-empty names
+      if (strlen(categoryNames[i]) > 0)
+        {
+        colorTable->SetAnnotation(vtkVariant(i), categoryNames[i]);
+        }
+      }
+    else
+      {
+      // Create default annotation
+      ss.str("");
+      ss.clear();
+      ss << "Category " << i;
+      colorTable->SetAnnotation(vtkVariant(i), ss.str());
+      }
+    }
+}
+
 //-----------------------------------------------------------------------------
 void vtkGDALRasterReader::PrintSelf(std::ostream& os, vtkIndent indent)
 {

IO/GDAL/vtkGDALRasterReader.h

@@ -67,21 +67,17 @@ public:
   // Get raster width and heigth
   vtkGetVector2Macro(RasterDimensions, int);
 
-  //BTX
   // Description:
   // Return metadata as reported by GDAL
   const std::vector<std::string>& GetMetaData();
-  //ETX
 
   // Description:
   // Return the invalid value for a pixel (for blanking purposes)
   double GetInvalidValue();
 
-  //BTX
   // Description:
   // Return domain metadata
   std::vector<std::string> GetDomainMetaData(const std::string& domain);
-  //ETX
 
   // Description:
   // Return driver name which was used to read the current data
@@ -117,7 +113,7 @@ protected:
   vtkGDALRasterReaderInternal* Implementation;
 
 private:
-  vtkGDALRasterReader(const vtkGDALRasterReader&); // Not implemented.
+  vtkGDALRasterReader(const vtkGDALRasterReader&); // Not implemented
   void operator=(const vtkGDALRasterReader&); // Not implemented
 };