//=========================================================================
//  Copyright (c) Kitware, Inc.
//  All rights reserved.
//  See LICENSE.txt 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 "vtkRawDEMReader.h"

#include "vtkAppendPolyData.h"
#include "vtkCellArray.h"
#include "vtkExecutive.h"
#include "vtkFloatArray.h"
#include "vtkGeoSphereTransform.h"
#include "vtkImageData.h"
#include "vtkInformation.h"
#include "vtkInformationVector.h"
#include "vtkMath.h"
#include "vtkNew.h"
#include "vtkObjectFactory.h"
#include "vtkPointData.h"
#include "vtkPolyData.h"
#include "vtkSmartPointer.h"
#include "vtkStreamingDemandDrivenPipeline.h"
#include "vtkStringArray.h"
#include "vtkTransform.h"
#include "vtkTriangle.h"

#include <algorithm>
#include <sys/stat.h>
#include <sys/types.h>
#include <vtksys/Glob.hxx>
#include <vtksys/SystemTools.hxx>

enum FileReadingStatus
{
  READ_OK = 0,
  READ_ERROR,
  READ_ABORT
};

vtkStandardNewMacro(vtkRawDEMReader);

vtkCxxSetObjectMacro(vtkRawDEMReader, Transform, vtkTransform);

struct RowCol
{
  RowCol(vtkIdType row, vtkIdType column)
  {
    this->Row = row;
    this->Column = column;
  }
  vtkIdType Row;
  vtkIdType Column;
};

class vtkRawDEMReaderInternals
{
public:
  int GetFilesInfo(const char* fileName, bool readSet);
  int ReadFileInfo(const char* fileName, RawDEMReaderFileInfo& fileInfo);

  std::vector<RawDEMReaderFileInfo> FileInfo;
  double LongitudeMidPoint;
  double LatitudeMidPoint;
  vtkIdType TotalNumberOfPoints;
  vtkIdType NumberOfColumns;
  vtkIdType NumberOfRows;
  vtkIdType OutputNumberOfColumns;
  vtkIdType OutputNumberOfRows;
  vtkIdType Extents[4];
  vtkIdType OnRatio;
  std::vector<RowCol> NoDataValueInstances;
};

bool SortPredicate(RawDEMReaderFileInfo elem1, RawDEMReaderFileInfo elem2)
{
  if (elem1.LatitudeOrigin > elem2.LatitudeOrigin)
  {
    return true;
  }
  else if (elem1.LatitudeOrigin == elem2.LatitudeOrigin &&
    elem1.LongitudeOrigin < elem2.LongitudeOrigin)
  {
    return true;
  }
  return false;
}

int vtkRawDEMReaderInternals::GetFilesInfo(const char* fileName, bool readSet)
{
  std::vector<std::string> files;
  if (readSet)
  {
    vtksys::Glob glob;
    std::string globString = vtksys::SystemTools::GetFilenamePath(fileName) + "/" +
      vtksys::SystemTools::GetFilenameWithoutLastExtension(fileName);
    globString.erase(globString.end() - 2, globString.end());
    globString += "*.hdr";
    glob.FindFiles(globString.c_str());
    files = glob.GetFiles();
    if (files.size() == 0)
    {
      return 0;
    }
  }
  else
  {
    files.push_back(fileName);
  }

  this->TotalNumberOfPoints = 0;
  this->FileInfo.clear();
  std::vector<std::string>::const_iterator file = files.begin();
  RawDEMReaderFileInfo fileInfo;
  double longitudeRange[2] = { 500, -500 }, latitudeRange[2] = { 500, -500 };
  for (file = files.begin(); file != files.end(); file++)
  {
    if (this->ReadFileInfo(file->c_str(), fileInfo) == READ_ERROR)
    {
      std::cerr << "Unable to open file: " << *file << endl;
      return READ_ERROR;
    }
    else
    {
      this->FileInfo.push_back(fileInfo);
    }

    this->TotalNumberOfPoints += fileInfo.NumberOfColumns * fileInfo.NumberOfRows;

    double longitude[2], latitude[2];
    longitude[0] = fileInfo.LongitudeOrigin;
    longitude[1] =
      fileInfo.LongitudeOrigin + (fileInfo.NumberOfColumns - 1) * fileInfo.PointSpacing;
    latitude[1] = fileInfo.LatitudeOrigin;
    latitude[0] = fileInfo.LatitudeOrigin - (fileInfo.NumberOfRows - 1) * fileInfo.PointSpacing;

    if (latitude[0] < latitudeRange[0])
    {
      latitudeRange[0] = latitude[0];
    }
    if (latitude[1] > latitudeRange[1])
    {
      latitudeRange[1] = latitude[1];
    }
    if (longitude[0] < longitudeRange[0])
    {
      longitudeRange[0] = longitude[0];
    }
    if (longitude[1] > longitudeRange[1])
    {
      longitudeRange[1] = longitude[1];
    }
  }

  std::sort(this->FileInfo.begin(), this->FileInfo.end(), SortPredicate);

  // figure out how many (total) rows & columns; ASSUMING that every row has same number
  // of columns
  std::vector<RawDEMReaderFileInfo>::iterator fileInfoIter = this->FileInfo.begin();
  double startLatitude = fileInfoIter->LatitudeOrigin;
  double previousLatitude = startLatitude;
  double startLongitude = fileInfoIter->LongitudeOrigin;
  this->NumberOfColumns = 0;
  this->NumberOfRows = 0;
  vtkIdType rowOffset = 0, columnOffset = 0;
  for (; fileInfoIter != this->FileInfo.end(); fileInfoIter++)
  {
    if (fileInfoIter->LatitudeOrigin != previousLatitude)
    {
      previousLatitude = fileInfoIter->LatitudeOrigin;
      rowOffset += (fileInfoIter - 1)->NumberOfRows;
      columnOffset = 0;
    }
    fileInfoIter->Offset[1] = rowOffset;
    fileInfoIter->Offset[0] = columnOffset;
    columnOffset += fileInfoIter->NumberOfColumns;

    // come up with full column and row counts
    if (fileInfoIter->LatitudeOrigin == startLatitude)
    {
      this->NumberOfColumns += fileInfoIter->NumberOfColumns;
    }
    if (fileInfoIter->LongitudeOrigin == startLongitude)
    {
      this->NumberOfRows += fileInfoIter->NumberOfRows;
    }
  }

  // one more time... adjust RowOffset to be from LL instead of UL
  fileInfoIter = this->FileInfo.begin();
  for (; fileInfoIter != this->FileInfo.end(); fileInfoIter++)
  {
    fileInfoIter->Offset[1] =
      this->NumberOfRows - fileInfoIter->NumberOfRows - fileInfoIter->Offset[1];
  }

  this->LatitudeMidPoint = 0.5 * (latitudeRange[0] + latitudeRange[1]);
  this->LongitudeMidPoint = 0.5 * (longitudeRange[0] + longitudeRange[1]);

  return READ_OK;
}

int vtkRawDEMReaderInternals::ReadFileInfo(const char* fileName, RawDEMReaderFileInfo& fileInfo)
{
  // read the "hdr" file
  std::string fileNameStr = vtksys::SystemTools::GetFilenamePath(fileName) + "/" +
    vtksys::SystemTools::GetFilenameWithoutLastExtension(fileName) + ".hdr";
  ifstream fin;

  fin.open(fileNameStr.c_str(), ios::binary);
  if (!fin)
  {
    return READ_ERROR;
  }

  fileInfo.FileName = fileName;

  char buffer[64];
  fin >> buffer >> fileInfo.NumberOfColumns;
  fin >> buffer >> fileInfo.NumberOfRows;
  fin >> buffer >> fileInfo.LongitudeOrigin;
  fin >> buffer >> fileInfo.LatitudeOrigin;
  fin >> buffer >> fileInfo.PointSpacing;
  fin >> buffer >> fileInfo.NoDataValue;
  fin >> buffer >> buffer;
  fin.close();

  // LatitudeOrigin is specified as lower edge, but the data is read in from top
  // edge, so adjust the origin to be "upper-left" corner
  fileInfo.LatitudeOrigin += (fileInfo.NumberOfRows - 1) * fileInfo.PointSpacing;

  return READ_OK;
}

vtkRawDEMReader::vtkRawDEMReader()
{
  this->FileName = NULL;
  this->FileInfoMode = true;
  this->OnRatio = 239;
  this->SetNumberOfInputPorts(0);
  this->RealNumberOfOutputPoints = 0;
  this->LimitReadToBounds = false;
  this->ReadBounds[0] = this->ReadBounds[2] = this->ReadBounds[4] = VTK_DOUBLE_MAX;
  this->ReadBounds[1] = this->ReadBounds[3] = this->ReadBounds[5] = VTK_DOUBLE_MIN;
  this->Transform = 0;
  this->TransformOutputData = false;

  this->ConvertFromLatLongToXYZ = true;
  this->TransformForZUp = true;
  this->RemoveCurvature = true;

  this->LatLongTransform1 = vtkSmartPointer<vtkGeoSphereTransform>::New();
  this->LatLongTransform2 = vtkSmartPointer<vtkTransform>::New();

  this->MaxNumberOfPoints = 1000000;
  this->LimitToMaxNumberOfPoints = false;
  this->TransformForZOrigin[0] = this->TransformForZOrigin[1] = -500;
  this->OutputImageData = true;
  this->ReadSetOfFiles = false;
  this->ZRotationAngle = 0.0;

  this->RowReadExtents[0] = this->RowReadExtents[1] = -1;
  this->ColumnReadExtents[0] = this->ColumnReadExtents[1] = -1;

  this->OutputDataTypeIsDouble = false;
  this->Internals = new vtkRawDEMReaderInternals;
  this->Origin[0] = this->Origin[1] = this->Origin[2] = 0;
}

vtkRawDEMReader::~vtkRawDEMReader()
{
  delete this->Internals;
  this->SetFileName(0);
  this->SetTransform(static_cast<vtkTransform*>(0));
}

// vtkSetStringMacro except we clear some variables if we update the value
void vtkRawDEMReader::SetFileName(const char* filename)
{
  vtkDebugMacro(<< this->GetClassName() << " (" << this << "): setting FileName to " << filename);
  if (this->FileName == NULL && filename == NULL)
  {
    return;
  }
  if (this->FileName && filename && !strcmp(this->FileName, filename))
  {
    return;
  }
  if (this->FileName)
  {
    delete[] this->FileName;
  }
  if (filename)
  {
    size_t n = strlen(filename) + 1;
    char* cp1 = new char[n];
    const char* cp2 = (filename);
    this->FileName = cp1;
    do
    {
      *cp1++ = *cp2++;
    } while (--n);
  }
  else
  {
    this->FileName = NULL;
  }

  // only want to clear these values if FileName changes!
  this->DataBounds[0] = this->DataBounds[2] = this->DataBounds[4] = VTK_DOUBLE_MAX;
  this->DataBounds[1] = this->DataBounds[3] = this->DataBounds[5] = VTK_DOUBLE_MIN;
  this->Modified();
}

void vtkRawDEMReader::SetTransform(double elements[16])
{
  vtkTransform* tmpTransform = vtkTransform::New();
  tmpTransform->SetMatrix(elements);
  this->SetTransform(tmpTransform);
  tmpTransform->Delete();
}

void vtkRawDEMReader::SetupReadExtents()
{
  vtkIdType readExtents[4] = { 0, this->Internals->NumberOfColumns - 1, 0,
    this->Internals->NumberOfRows - 1 };
  if (this->RowReadExtents[0] >= 0 && this->RowReadExtents[0] < this->Internals->NumberOfRows)
  {
    readExtents[2] = this->RowReadExtents[0];
  }
  if (this->RowReadExtents[1] >= readExtents[2] &&
    this->RowReadExtents[1] < this->Internals->NumberOfRows)
  {
    readExtents[3] = this->RowReadExtents[1];
  }
  if (this->ColumnReadExtents[0] >= 0 &&
    this->ColumnReadExtents[0] < this->Internals->NumberOfColumns)
  {
    readExtents[0] = this->ColumnReadExtents[0];
  }
  if (this->ColumnReadExtents[1] >= readExtents[0] &&
    this->ColumnReadExtents[1] < this->Internals->NumberOfColumns)
  {
    readExtents[1] = this->ColumnReadExtents[1];
  }

  this->Internals->OnRatio = this->OnRatio;
  if (this->LimitToMaxNumberOfPoints)
  {
    if (this->MaxNumberOfPoints >= this->Internals->TotalNumberOfPoints)
    {
      this->Internals->OnRatio = 1;
    }
    else
    {
      this->Internals->OnRatio =
        sqrt(this->Internals->TotalNumberOfPoints / static_cast<double>(this->MaxNumberOfPoints));
    }
  }
  int onRatio = static_cast<int>(floor(sqrt(static_cast<double>(this->Internals->OnRatio))));

  int numberOfRangeRows = readExtents[3] - readExtents[2] + 1;
  int numberOfRangeColumns = readExtents[1] - readExtents[0] + 1;
  this->Internals->OutputNumberOfColumns = 1 + (numberOfRangeColumns - 1) / onRatio;
  this->Internals->OutputNumberOfRows = 1 + (numberOfRangeRows - 1) / onRatio;

  this->RealNumberOfOutputPoints =
    this->Internals->OutputNumberOfRows * this->Internals->OutputNumberOfColumns;

  this->Internals->Extents[0] = readExtents[0];
  this->Internals->Extents[1] =
    readExtents[0] + (this->Internals->OutputNumberOfColumns - 1) * onRatio;
  this->Internals->Extents[3] = readExtents[3];
  this->Internals->Extents[2] =
    readExtents[3] - (this->Internals->OutputNumberOfRows - 1) * onRatio;
}

int vtkRawDEMReader::RequestData(vtkInformation* vtkNotUsed(request),
  vtkInformationVector** vtkNotUsed(inputVector), vtkInformationVector* outputVector)
{
  int res = this->Internals->GetFilesInfo(this->FileName, this->ReadSetOfFiles);
  if (res == READ_ERROR)
  {
    return 0;
  }
  else if (res == READ_ABORT)
  {
    return 1;
  }

  this->SetupReadExtents();

  // setup the ReadBBox, IF we're limiting the read to specifed ReadBounds
  if (this->LimitReadToBounds)
  {
    this->ReadBBox.Reset();
    this->ReadBBox.SetMinPoint(this->ReadBounds[0], this->ReadBounds[2], this->ReadBounds[4]);
    this->ReadBBox.SetMaxPoint(this->ReadBounds[1], this->ReadBounds[3], this->ReadBounds[5]);
    // the ReadBBox is guaranteed to be "valid", regardless of the whether
    // ReadBounds is valid.  If any of the MonPoint values are greater than
    // the corresponding MaxPoint, the MinPoint component will be set to be
    // the same as the MaxPoint during the SetMaxPoint fn call.
  }

  vtkInformation* outInfo = outputVector->GetInformationObject(0);
  vtkPolyData* pdOutput = vtkPolyData::SafeDownCast(outInfo->Get(vtkDataObject::DATA_OBJECT()));
  if (pdOutput)
  {
    this->ReadPolyDataOutput(pdOutput);
    //pdOutput->GetBounds(this->DataBounds);
  }
  else
  {
    vtkImageData* imageOutput =
      vtkImageData::SafeDownCast(outInfo->Get(vtkDataObject::DATA_OBJECT()));
    this->ReadImageDataOutput(imageOutput);
    //imageOutput->GetBounds(this->DataBounds);
  }
  return 1;
}

void vtkRawDEMReader::ReadPolyDataOutput(vtkPolyData* output)
{
  if (!this->OutputImageData && this->ConvertFromLatLongToXYZ && this->TransformForZUp)
  {
    this->LatLongTransform2->Identity();
    double rotationAxis[3], zAxis[3] = { 0, 0, 1 };
    double pt[3];
    RawDEMReaderFileInfo& fileInfo = this->Internals->FileInfo.front();

    // lat/long as specified, height of 0
    if (this->TransformForZOrigin[0] > -180)
    {
      pt[0] = this->TransformForZOrigin[0];
      pt[1] = this->TransformForZOrigin[1];
    }
    else
    {
      pt[0] = this->Internals->LongitudeMidPoint;
      pt[1] = this->Internals->LatitudeMidPoint;
    }
    pt[2] = 0;
    this->LatLongTransform1->TransformPoint(pt, pt);
    double rotationAngle;
    double planeNormal[3] = { pt[0], pt[1], pt[2] };
    vtkMath::Normalize(planeNormal);
    vtkMath::Cross(planeNormal, zAxis, rotationAxis);
    rotationAngle = vtkMath::DegreesFromRadians(acos(planeNormal[2]));

    // Figure out Z-rotation
    if (this->ZRotationAngle == 0)
    {
      double pt1[3] = { fileInfo.LongitudeOrigin, fileInfo.LatitudeOrigin, 0 };
      double pt3[3] = { fileInfo.LongitudeOrigin,
        fileInfo.LatitudeOrigin - (this->Internals->NumberOfRows - 1) * fileInfo.PointSpacing, 0 };
      this->LatLongTransform1->TransformPoint(pt1, pt1);
      this->LatLongTransform1->TransformPoint(pt3, pt3);

      vtkNew<vtkTransform> tmpTransform;
      tmpTransform->PreMultiply();
      tmpTransform->RotateWXYZ(rotationAngle, rotationAxis);
      tmpTransform->Translate(-pt[0], -pt[1], -pt[2]);
      tmpTransform->TransformPoint(pt3, pt3);
      tmpTransform->TransformPoint(pt1, pt1);

      double temppt[3] = { pt1[0] - pt3[0], pt1[1] - pt3[1], pt1[2] - pt3[2] };
      this->ZRotationAngle = atan2(temppt[0], temppt[1]) * 180 / vtkMath::Pi();
    }

    this->LatLongTransform2->PreMultiply();
    //    this->LatLongTransform2->Translate(pt[0], pt[1], pt[2]);
    if (this->ZRotationAngle != 0)
    {
      this->LatLongTransform2->RotateZ(this->ZRotationAngle);
    }
    this->LatLongTransform2->RotateWXYZ(rotationAngle, rotationAxis);
    this->LatLongTransform2->Translate(-pt[0], -pt[1], -pt[2]);
  }

  vtkPoints* newPts = vtkPoints::New();
  if (this->OutputDataTypeIsDouble)
  {
    newPts->SetDataTypeToDouble();
  }
  else
  {
    newPts->SetDataTypeToFloat();
  }
  vtkCellArray* newVerts = vtkCellArray::New();
  output->SetPoints(newPts);
  output->SetVerts(newVerts);
  newPts->UnRegister(this);
  newVerts->UnRegister(this);

  vtkIdType totalNumberOfPoints =
    this->Internals->OutputNumberOfRows * this->Internals->OutputNumberOfColumns;
  newPts->Allocate(totalNumberOfPoints);
  newVerts->Allocate(2 * totalNumberOfPoints);

  std::vector<RawDEMReaderFileInfo>::iterator file = this->Internals->FileInfo.begin();
  for (; file != this->Internals->FileInfo.end(); file++)
  {
    this->ReadData(*file, 0, newPts, newVerts);
  }
  newVerts->Squeeze();
  newPts->Squeeze();
}

void vtkRawDEMReader::ReadImageDataOutput(vtkImageData* output)
{
  vtkImageData* data = this->AllocateOutputData(output);
  vtkFloatArray* scalars = vtkFloatArray::SafeDownCast(data->GetPointData()->GetScalars());

  this->Internals->NoDataValueInstances.clear();

  std::vector<RawDEMReaderFileInfo>::iterator file = this->Internals->FileInfo.begin();
  for (; file != this->Internals->FileInfo.end(); file++)
  {
    this->ReadData(*file, scalars, 0, 0);
  }

  // fixup any "NoDataValue" instances
  std::vector<RowCol>::const_iterator noDataValueIter =
    this->Internals->NoDataValueInstances.begin();
  float noDataValue = this->Internals->FileInfo.front().NoDataValue;
  for (; noDataValueIter != this->Internals->NoDataValueInstances.end(); noDataValueIter++)
  {
    double sum = 0;
    vtkIdType sumCount = 0;
    for (vtkIdType row = noDataValueIter->Row - 1; row < noDataValueIter->Row + 2; row++)
    {
      if (row < 0 || row >= this->Internals->OutputNumberOfRows)
      {
        continue;
      }
      for (vtkIdType col = noDataValueIter->Column - 1; col < noDataValueIter->Column + 2; col++)
      {
        if (col < 0 || col >= this->Internals->OutputNumberOfColumns ||
          (col == noDataValueIter->Column && row == noDataValueIter->Row))
        {
          continue;
        }
        float value = scalars->GetValue(row * this->Internals->OutputNumberOfColumns + col);
        if (value != noDataValue)
        {
          sum += value;
          sumCount++;
        }
      }
      if (sumCount)
      {
        scalars->SetValue(
          noDataValueIter->Row * this->Internals->OutputNumberOfColumns + noDataValueIter->Column,
          sum / sumCount);
      }
    }
  }
}

void vtkRawDEMReader::ReadData(
  RawDEMReaderFileInfo& fileInfo, vtkFloatArray* scalars, vtkPoints* pts, vtkCellArray* verts)
{
  ifstream fin;

  // read the binary "FLT" file
  std::string fileNameStr = vtksys::SystemTools::GetFilenamePath(fileInfo.FileName) + "/" +
    vtksys::SystemTools::GetFilenameWithoutLastExtension(fileInfo.FileName) + ".FLT";

  fin.open(fileNameStr.c_str(), ios::binary);

  // is there even anything to read from this file
  if (this->Internals->Extents[2] >= fileInfo.Offset[1] + fileInfo.NumberOfRows ||
    this->Internals->Extents[3] < fileInfo.Offset[1] ||
    this->Internals->Extents[0] >= fileInfo.Offset[0] + fileInfo.NumberOfColumns ||
    this->Internals->Extents[1] < fileInfo.Offset[0])
  {
    return;
  }

  // from onRatio and Row/Column offsets, determine if we need to skip at the
  // beginning rows / columns
  vtkIdType onRatio =
    static_cast<vtkIdType>(floor(sqrt(static_cast<double>(this->Internals->OnRatio))));
  vtkIdType topRowRelativeToExtent =
    this->Internals->Extents[3] - (fileInfo.Offset[1] + fileInfo.NumberOfRows - 1);
  vtkIdType initialRowSkip;
  if (topRowRelativeToExtent > 0)
  {
    initialRowSkip =
      topRowRelativeToExtent % onRatio == 0 ? 0 : onRatio - (topRowRelativeToExtent % onRatio);
  }
  else
  {
    initialRowSkip = (fileInfo.Offset[1] + fileInfo.NumberOfRows - 1) - this->Internals->Extents[3];
  }
  if (initialRowSkip)
  {
    fin.seekg(4 * initialRowSkip * fileInfo.NumberOfColumns, ios::cur);
  }
  vtkIdType leftColumnRelativeToExtent = fileInfo.Offset[0] - this->Internals->Extents[0];
  vtkIdType initialColumnSkip;
  if (leftColumnRelativeToExtent > 0)
  {
    initialColumnSkip = leftColumnRelativeToExtent % onRatio == 0 ? 0 : onRatio -
        (leftColumnRelativeToExtent % onRatio);
  }
  else
  {
    initialColumnSkip = this->Internals->Extents[0] - fileInfo.Offset[0];
  }

  // how many rows/ columns do we not read at the end
  vtkIdType lastPossibleRow = fileInfo.NumberOfRows - 1;
  if (this->Internals->Extents[2] > fileInfo.Offset[1])
  {
    lastPossibleRow -= this->Internals->Extents[2] - fileInfo.Offset[1];
  }
  vtkIdType numberOfRows = 1 + (lastPossibleRow - initialRowSkip) / onRatio;

  vtkIdType lastPossibleColumn = fileInfo.NumberOfColumns - 1;
  if (this->Internals->Extents[1] < fileInfo.Offset[0] + fileInfo.NumberOfColumns - 1)
  {
    lastPossibleColumn -=
      (fileInfo.Offset[0] + fileInfo.NumberOfColumns - 1) - this->Internals->Extents[1];
  }
  vtkIdType numberOfColumns = 1 + (lastPossibleColumn - initialColumnSkip) / onRatio;

  // how much do we skip at the end of each row
  vtkIdType lastActualColumn = initialColumnSkip + (numberOfColumns - 1) * onRatio;
  vtkIdType endRowSkip = (fileInfo.NumberOfColumns - 1) - lastActualColumn;
  // should be integer (no remainder) results
  vtkIdType outputRowOffset = ((fileInfo.Offset[1] + fileInfo.NumberOfRows - 1) - initialRowSkip -
                                this->Internals->Extents[2]) /
    onRatio;
  vtkIdType outputColumnOffset =
    (fileInfo.Offset[0] + initialColumnSkip - this->Internals->Extents[0]) / onRatio;

  // since data comes in from upper-left, we start with 1st row for output
  // as "top" row... thus, due to skipping, we might not output some # of
  // bottom rows, and thus need to adjust origin accordingly
  double dataOrigin[3] = { /*this->Internals->LongitudeMidPoint + */ fileInfo.LongitudeOrigin +
      static_cast<double>(this->Internals->Extents[0]) * fileInfo.PointSpacing,
    /*this->Internals->LatitudeMidPoint + */ fileInfo.LatitudeOrigin -
      static_cast<double>(this->Internals->NumberOfRows - 1 - this->Internals->Extents[2]) *
        fileInfo.PointSpacing,
    0 };

  double rowProgressDelta = static_cast<double>(numberOfColumns) /
    (static_cast<double>(this->Internals->OutputNumberOfColumns) *
                              static_cast<double>(this->Internals->OutputNumberOfRows));

  vtkBoundingBox bbox;
  float rawData;
  for (vtkIdType row = 0; row < numberOfRows; row++)
  {
    vtkIdType outputImageOffset =
      outputColumnOffset + (outputRowOffset - row) * this->Internals->OutputNumberOfColumns;
    if (initialColumnSkip)
    {
      fin.ignore(4 * initialColumnSkip);
    }
    for (vtkIdType column = 0; column < numberOfColumns; column++)
    {
      fin.read(reinterpret_cast<char*>(&rawData), 4);

      if (scalars)
      {
        // set, regardless of whether "NoDataValue", but if NoDataValue,
        // save location to allow fixing
        scalars->SetValue(outputImageOffset + column, rawData);
        if (rawData == fileInfo.NoDataValue)
        {
          RowCol tmp(outputRowOffset - row, outputColumnOffset + column);
          this->Internals->NoDataValueInstances.push_back(tmp);
        }
      }
      else
      {
        double pt[3];
        pt[0] = dataOrigin[0] +
          fileInfo.PointSpacing * static_cast<double>(column * onRatio + initialColumnSkip);
        pt[1] = dataOrigin[1] -
          fileInfo.PointSpacing * static_cast<double>(row * onRatio + initialRowSkip);
        pt[2] = rawData;

        if (pt[2] != fileInfo.NoDataValue)
        {
          if (this->ConvertFromLatLongToXYZ)
          {
            this->LatLongTransform1->TransformPoint(pt, pt);
            if (this->TransformForZUp)
            {
              this->LatLongTransform2->TransformPoint(pt, pt);
              if (this->RemoveCurvature)
              {
                pt[2] = rawData;
              }
            }
          }
          pt[0] -= Origin[0];
          pt[1] -= Origin[1];
          pt[2] -= Origin[2];
          bbox.AddPoint(pt[0], pt[1], pt[2]);

          // add the point, but 1st make sure it is in the ReadBounds (if specified);
          // consider the Transform if set (and "on")
          double transformedPt[3];
          bool addPt = true;
          if (this->Transform)
          {
            // only need the transformed pt if we're limiting read based on bounds or
            // we're transforming the output
            if (this->LimitReadToBounds || this->TransformOutputData)
            {
              this->Transform->TransformPoint(pt, transformedPt);
            }
            if (this->LimitReadToBounds &&
              !this->ReadBBox.ContainsPoint(transformedPt[0], transformedPt[1], transformedPt[2]))
            {
              addPt = false;
            }
          }
          else // not transformed, use as read in
          {
            if (this->LimitReadToBounds && !this->ReadBBox.ContainsPoint(pt[0], pt[1], pt[2]))
            {
              addPt = false;
            }
          }

          if (addPt)
          {
            vtkIdType outputIdx;
            if (this->Transform && this->TransformOutputData)
            {
              outputIdx = pts->InsertNextPoint(transformedPt);
            }
            else
            {
              //double tmppts[3]={pt[0], pt[2], pt[1]};
              outputIdx = pts->InsertNextPoint(pt);
            }
            verts->InsertNextCell(1, &outputIdx);
          }
        }
      }

      if (column == numberOfColumns - 1 && endRowSkip) // read the rest of the row
      {
        fin.ignore(4 * endRowSkip);
      }
      else if (onRatio > 1)
      {
        fin.ignore(4 * (onRatio - 1));
      }
    }
    this->UpdateProgress(this->GetProgress() + rowProgressDelta);
    //    this->UpdateProgress( (double)(row + 1) / (double)numberOfRows );
    if (this->GetAbortExecute())
    {
      break;
    }

    if (onRatio > 1)
    {
      if (row < numberOfRows - 1) // skip over appropriate # of rows
      {
        fin.seekg(4 * (onRatio - 1) * fileInfo.NumberOfColumns, ios::cur);
      }
    }
  }
  fin.close();
  this->UpdateProgress(1.0);

  //// set our DataBounds
  bbox.GetBounds(this->DataBounds);
}

vtkIdType vtkRawDEMReader::GetTotalNumberOfPoints()
{
  if (this->Internals->GetFilesInfo(this->FileName, this->ReadSetOfFiles) == READ_ERROR)
  {
    return -1;
  }
  this->SetupReadExtents();

  return this->Internals->TotalNumberOfPoints;
}

vtkIdType vtkRawDEMReader::GetRealNumberOfOutputPoints()
{
  if (this->Internals->GetFilesInfo(this->FileName, this->ReadSetOfFiles) == READ_ERROR)
  {
    return -1;
  }
  this->SetupReadExtents();

  return this->RealNumberOfOutputPoints;
}

void vtkRawDEMReader::GatherDimensions()
{
  if (this->Internals->GetFilesInfo(this->FileName, this->ReadSetOfFiles) == READ_ERROR)
  {
    this->Dimensions[0] = -1;
    this->Dimensions[1] = -1;
    return;
  }
  this->Dimensions[0] = this->Internals->NumberOfColumns;
  this->Dimensions[1] = this->Internals->NumberOfRows;
}

void vtkRawDEMReader::PrintSelf(ostream& os, vtkIndent indent)
{
  this->Superclass::PrintSelf(os, indent);

  os << indent << "File Name: " << (this->FileName ? this->FileName : "(none)") << "\n";
  os << indent
     << "Convert From Lat/Long to xyz: " << (this->ConvertFromLatLongToXYZ ? "On" : "Off");
}

int vtkRawDEMReader::RequestDataObject(vtkInformation*,
  vtkInformationVector** vtkNotUsed(inputVector), vtkInformationVector* outputVector)
{
  int outputType = VTK_POLY_DATA;
  if (this->OutputImageData)
  {
    outputType = VTK_IMAGE_DATA;
  }

  vtkInformation* info = outputVector->GetInformationObject(0);
  vtkDataSet* output = vtkDataSet::SafeDownCast(info->Get(vtkDataObject::DATA_OBJECT()));

  if (output && (output->GetDataObjectType() == outputType))
  {
    return 1;
  }

  if (!output || output->GetDataObjectType() != outputType)
  {
    switch (outputType)
    {
      case VTK_POLY_DATA:
        output = vtkPolyData::New();
        break;
      case VTK_IMAGE_DATA:
        output = vtkImageData::New();
        break;
      default:
        return 0;
    }

    this->GetExecutive()->SetOutputData(0, output);
    output->Delete();
  }

  return 1;
}

int vtkRawDEMReader::RequestInformation(vtkInformation* vtkNotUsed(request),
  vtkInformationVector** vtkNotUsed(inputVector), vtkInformationVector* outputVector)
{
  if (!this->FileName)
  {
    vtkErrorMacro("FileName has to be specified!");
    return 0;
  }

  // get the info objects, needed for vtkImageData output
  if (this->OutputImageData && vtksys::SystemTools::FileExists(this->FileName))
  {
    int res = this->Internals->GetFilesInfo(this->FileName, this->ReadSetOfFiles);
    if (res == READ_ERROR)
    {
      return 0;
    }
    else if (res == READ_ABORT)
    {
      return 1;
    }

    if (this->Internals->FileInfo.size() == 0)
    {
      return 0;
    }
    this->SetupReadExtents();

    RawDEMReaderFileInfo& fileInfo = this->Internals->FileInfo.front();
    // since data comes in from upper-left, we start with 1st row for output
    // as "top" row... thus, due to skipping, we might not output some # of
    // bottom rows, and thus need to adjust origin accordingly
    double dataOrigin[3] = { fileInfo.LongitudeOrigin +
        this->Internals->Extents[0] * fileInfo.PointSpacing,
      fileInfo.LatitudeOrigin -
        (this->Internals->NumberOfRows - 1 - this->Internals->Extents[2]) * fileInfo.PointSpacing,
      0 };

    vtkInformation* outInfo = outputVector->GetInformationObject(0);
    int dataExtent[6] = { 0, static_cast<int>(this->Internals->OutputNumberOfColumns - 1), 0,
      static_cast<int>(this->Internals->OutputNumberOfRows - 1), 0, 0 };
    outInfo->Set(vtkStreamingDemandDrivenPipeline::WHOLE_EXTENT(), dataExtent, 6);
    double spacing[3] = { fileInfo.PointSpacing * this->Internals->OnRatio,
      fileInfo.PointSpacing * this->Internals->OnRatio, 0 };
    outInfo->Set(vtkDataObject::SPACING(), spacing, 3);

    vtkDataObject::SetPointDataActiveScalarInfo(outInfo, VTK_FLOAT, 1);
    outInfo->Set(vtkDataObject::ORIGIN(), dataOrigin, 3);
  }

  return 1;
}

vtkImageData* vtkRawDEMReader::AllocateOutputData(vtkDataObject* output)
{
  // set the extent to be the update extent
  vtkImageData* out = vtkImageData::SafeDownCast(output);
  if (out)
  {
    // this needs to be fixed -Ken
    vtkStreamingDemandDrivenPipeline* sddp =
      vtkStreamingDemandDrivenPipeline::SafeDownCast(this->GetExecutive());
    int numInfoObj = sddp->GetNumberOfOutputPorts();
    if (sddp && numInfoObj == 1)
    {
      int extent[6];
      sddp->GetOutputInformation(0)->Get(vtkStreamingDemandDrivenPipeline::UPDATE_EXTENT(), extent);
      out->SetExtent(extent);
    }
    else
    {
      vtkWarningMacro("There are multiple output ports. You cannot use AllocateOutputData");
      return NULL;
    }
    out->AllocateScalars(VTK_FLOAT, 1);
  }
  return out;
}
