// -*- c++ -*-
/*=========================================================================

  Program:   Visualization Toolkit
  Module:    vtkNrrdReader.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.

=========================================================================*/
/*----------------------------------------------------------------------------
 Copyright (c) Sandia Corporation
 See Copyright.txt or http://www.paraview.org/HTML/Copyright.html for details.
----------------------------------------------------------------------------*/

#include "vtkNrrdReader.h"

#include "vtkCharArray.h"
#include "vtkImageData.h"
#include "vtkInformationVector.h"
#include "vtkObjectFactory.h"
#include "vtkStringArray.h"

#include <string>
#include <algorithm>
#include <vector>
#include <vtksys/SystemTools.hxx>
#include <istream>
#include <sstream>

#include <math.h>
#include <string.h>
#include <ctype.h>

#include "vtkSmartPointer.h"

#define VTK_CREATE(type, name) \
  vtkSmartPointer<type> name = vtkSmartPointer<type>::New()

//=============================================================================
static std::string trim(std::string s)
{
  size_t start = 0;
  while ((start < s.length()) && (isspace(s[start])))
    {
    start++;
    }
  size_t end = s.length();
  while ((end > start) && (isspace(s[end-1])))
    {
    end--;
    }
  return s.substr(start, end-start);
}

//-----------------------------------------------------------------------------
static std::vector<std::string> split(std::string s)
{
  std::vector<std::string> result;
  size_t startValue = 0;
  while (true)
    {
    while ((startValue < s.length()) && isspace(s[startValue]))
      {
      startValue++;
      }
    if (startValue >= s.length()) return result;
    size_t endValue = startValue;
    while ((endValue < s.length()) && !isspace(s[endValue]))
      {
      endValue++;
      }
    result.push_back(s.substr(startValue, endValue-startValue));
    startValue = endValue;
    }
}

//-----------------------------------------------------------------------------
static void GetVector(std::string s, std::vector<int> &dest)
{
  std::vector<std::string> strlist = split(s);
  for (size_t i = 0; i < dest.size(); i++)
    {
    if (i < strlist.size())
      {
      dest[i] = atoi(strlist[i].c_str());
      }
    else
      {
      dest[i] = 0;
      }
    }
}

//-----------------------------------------------------------------------------
static void GetVector(std::string s, std::vector<double> &dest)
{
  std::vector<std::string> strlist = split(s);
  for (size_t i = 0; i < dest.size(); i++)
    {
    if (i < strlist.size())
      {
      dest[i] = atof(strlist[i].c_str());
      }
    else
      {
      dest[i] = 0.0;
      }
    }
}

//-----------------------------------------------------------------------------
static std::vector<double> ParseVector(std::string s)
{
  std::vector<double> result;

  s = trim(s);
  if ((s[0] != '(') || (s[s.length()-1] != ')')) return result;
  s = s.substr(1, s.length()-2);
  while (true)
    {
    size_t i = s.find(',');
    std::string value = s.substr(0, i);
    result.push_back(atof(value.c_str()));
    if (i == s.npos) break;
    s = s.substr(i+1);
    }

  return result;
}

//-----------------------------------------------------------------------------
static int NrrdType2VTKType(std::string nrrdType)
{
  nrrdType = trim(nrrdType);
  if (   (nrrdType == "signed char")
      || (nrrdType == "int8") || (nrrdType == "int8_t") )
    {
    return VTK_CHAR;
    }
  else if (   (nrrdType == "uchar") || (nrrdType == "unsigned char")
           || (nrrdType == "uint8") || (nrrdType == "uint8_t") )
    {
    return VTK_UNSIGNED_CHAR;
    }
  else if (   (nrrdType == "short") || (nrrdType == "short int")
           || (nrrdType == "signed short") || (nrrdType == "signed short int")
           || (nrrdType == "int16") || (nrrdType == "int16_t") )
    {
    return VTK_SHORT;
    }
  else if (   (nrrdType == "ushort") || (nrrdType == "unsigned short")
           || (nrrdType == "unsigned short int")
           || (nrrdType == "uint16") || (nrrdType == "uint16_t") )
    {
    return VTK_UNSIGNED_SHORT;
    }
  else if (   (nrrdType == "int") || (nrrdType == "signed int")
           || (nrrdType == "int32") || (nrrdType == "int32_t") )
    {
    return VTK_INT;
    }
  else if (   (nrrdType == "uint") || (nrrdType == "unsigned int")
           || (nrrdType == "uint32") || (nrrdType == "uint32_t") )
    {
    return VTK_UNSIGNED_INT;
    }
  else if (   (nrrdType == "longlong") || (nrrdType == "long long")
           || (nrrdType == "long long int") || (nrrdType == "signed long long")
           || (nrrdType == "signed long long int") || (nrrdType == "int64")
           || (nrrdType == "int64_t") )
    {
    return VTK_TYPE_INT64;
    }
  else if (   (nrrdType == "ulonglong") || (nrrdType == "unsigned long long")
           || (nrrdType == "unsigned long long int")
           || (nrrdType == "uint64") || (nrrdType == "uint64_t") )
    {
    return VTK_TYPE_UINT64;
    }
  else if (nrrdType == "float")
    {
    return VTK_FLOAT;
    }
  else if (nrrdType == "double")
    {
    return VTK_DOUBLE;
    }
  else if (nrrdType == "block")
    {
    vtkGenericWarningMacro(<< "Reading blocks not supported.");
    return VTK_VOID;
    }
  else
    {
    vtkGenericWarningMacro(<< "Unknown type: '" << nrrdType << "'");
    return VTK_VOID;
    }
}

//-----------------------------------------------------------------------------

vtkObjectFactoryNewMacro(vtkNrrdReader);

//-----------------------------------------------------------------------------
vtkNrrdReader::vtkNrrdReader()
{
  this->DataFiles = vtkStringArray::New();
  this->Encoding = ENCODING_RAW;
}

vtkNrrdReader::~vtkNrrdReader()
{
  this->DataFiles->Delete();
  this->DataFiles = NULL;
}

void vtkNrrdReader::PrintSelf(ostream &os, vtkIndent indent)
{
  this->Superclass::PrintSelf(os, indent);
  os << indent << "Encoding: ";
  switch (this->Encoding)
    {
    case ENCODING_RAW: os << "raw" << endl; break;
    case ENCODING_ASCII: os << "ascii" << endl; break;
    default: os << "UNKNOWN!" << endl; break;
    }
}





//-----------------------------------------------------------------------------
int vtkNrrdReader::CanReadFile(const char *filename)
{
  ifstream file(filename, ios::in | ios::binary);
  std::string firstLine;
  getline(file, firstLine);
  if (firstLine.substr(0, 4) == "NRRD")
    {
    return 2;
    }
  else
    {
    return 0;
    }
}

//-----------------------------------------------------------------------------
int vtkNrrdReader::RequestInformation(vtkInformation *request,
                                      vtkInformationVector **inputVector,
                                      vtkInformationVector *outputVector)
{
  if (!this->ReadHeader()) return 0;

  return this->Superclass::RequestInformation(request,
                                              inputVector, outputVector);
}

//-----------------------------------------------------------------------------
int vtkNrrdReader::ReadHeaderInternal(vtkCharArray* headerBuffer)
{
  if (!this->FileName)
    {
    vtkErrorMacro(<< "No filename set.");
    return 0;
    }

  ifstream file(this->FileName, ios::in | ios::binary);
  // Read in 4 MB.  Assuming that the header will be smaller than that.
  headerBuffer->SetNumberOfTuples(0x400000);
  file.read(headerBuffer->GetPointer(0), 0x400000-1);
  vtkIdType buffersize = file.gcount();
  headerBuffer->SetValue(buffersize, '\0');
  headerBuffer->SetNumberOfTuples(buffersize+1);

  // Find a blank line (which signals the end of the header).  Be careful
  // because line endings can be "\n", "\r\n", or both.  It is possible that
  // the entire file is the header (happens with "detached headers").
  char *bufferStart = headerBuffer->GetPointer(0);
  char *s = bufferStart;
  while ((s = strchr(s+1, '\n')) != NULL)
    {
    // If you find a double line ending, terminate the string and shorten the
    // buffer.
    if (s[1] == '\n')
      {
      s[2] = '\0';
      headerBuffer->SetNumberOfTuples(
                                 static_cast<vtkIdType>(s + 3 - bufferStart));
      break;
      }
    if ((s[1] == '\r') && (s[2] == '\n'))
      {
      s[3] = '\0';
      headerBuffer->SetNumberOfTuples(
                                 static_cast<vtkIdType>(s + 4 - bufferStart));
      break;
      }
    }

  return 1;
}


//-----------------------------------------------------------------------------
int vtkNrrdReader::ReadHeader()
{
  VTK_CREATE(vtkCharArray, headerBuffer);

  if(!this->ReadHeaderInternal(headerBuffer))
    {
    return 0;
    }

  return this->ReadHeader(headerBuffer);
}

//-----------------------------------------------------------------------------
int vtkNrrdReader::ReadHeader(vtkCharArray *headerBuffer)
{
  this->HeaderSize = headerBuffer->GetNumberOfTuples();

  std::string headerStringBuffer = headerBuffer->GetPointer(0);
  std::stringstream header(headerStringBuffer);

  std::string line;
  getline(header, line);
  if (line.substr(0, 4) != "NRRD")
    {
    vtkErrorMacro(<< this->FileName << " is not a nrrd file.");
    return 0;
    }

  this->DataFiles->Initialize();
  int numDimensions = 0;
  int subDimension = -1;
  std::vector<int> dimSizes;
  std::vector<double> dimSpacing;
  this->FileLowerLeft = 1;
  this->Encoding = ENCODING_RAW;
  while (1)
    {
    getline(header, line);
    if (line.length() < 1) break;

    if (line[0] == '#')
      {
      // Comment.  Ignore.
      continue;
      }

    size_t delm = line.find(": ");
    if (delm != line.npos)
      {
      // A field/description pair.
      std::string field = line.substr(0, delm);
      std::string description = trim(line.substr(delm+2));
      if (field == "dimension")
        {
        numDimensions = atoi(description.c_str());
        }
      else if (field == "sizes")
        {
        dimSizes.resize(numDimensions);
        GetVector(description, dimSizes);
        }
      else if (field == "spacings")
        {
        dimSpacing.resize(numDimensions);
        GetVector(description, dimSpacing);
        }
      else if (field == "type")
        {
        this->DataScalarType = NrrdType2VTKType(description);
        if (this->DataScalarType == VTK_VOID) return 0;
        }
      else if (field == "encoding")
        {
        std::transform(description.begin(), description.end(),
                       description.begin(),
                       ::tolower);
        if (description == "raw")
          {
          this->Encoding = ENCODING_RAW;
          }
        else if (description == "ascii"
                 || description == "txt"
                 || description == "text")
          {
          this->Encoding = ENCODING_ASCII;
          }
        else
          {
          vtkErrorMacro(<< "Unsupported encoding: " << description);
          return 0;
          }
        }
      else if ((field == "data file") || (field == "datafile"))
        {
        std::vector<std::string> filepatterninfo = split(description);
        if (filepatterninfo[0] == "LIST")
          {
          // After LIST there is an optional subdimension (see next case below).
          subDimension = (  (filepatterninfo.size() > 1)
                          ? atoi(filepatterninfo[1].c_str()) : numDimensions );

          // In this mode files are listed one per line to the end of the file.
          while (1)
            {
            getline(header, line);
            trim(line);
            if (line.length() < 1) break;
            this->DataFiles->InsertNextValue(line);
            }
          break;
          }
        else if (filepatterninfo.size() >=4)
          {
          // description should be "<format> <min> <max> <step> [<subdim>]"
          // where <format> is a string to be processed by sprintf and <min>,
          // <max>, and <step> form the numbers.  <subdim> defines on which
          // dimension the files are split up.
          std::string format = filepatterninfo[0];
          int min = atoi(filepatterninfo[1].c_str());
          int max = atoi(filepatterninfo[2].c_str());
          int step = atoi(filepatterninfo[3].c_str());
          subDimension = (  (filepatterninfo.size() > 4)
                          ? atoi(filepatterninfo[4].c_str()) : numDimensions );
          char *filename = new char[format.size() + 20];
          for (int i = min; i <= max; i += step)
            {
            sprintf(filename, format.c_str(), i);
            this->DataFiles->InsertNextValue(filename);
            }
          delete[] filename;
          }
        else
          {
          // description is simply a filename
          this->DataFiles->InsertNextValue(description);
          }
        }
      else if (field == "space")
        {
        // All spaces are either 3D or 3D with time.
        if (description.find("time") != description.npos)
          {
          vtkErrorMacro(<< "Time in NRRD array not supported (yet).");
          return 0;
          }
        if (   (description == "left-anterior-superior")
            || (description == "LAS")
            || (description == "3D-left-handed") )
          {
          this->FileLowerLeft = 0;
          }
        numDimensions = 3;
        }
      else if (field == "labels")
        {
        std::string dataname = description.substr(description.find("\"")+1);
        dataname = dataname.substr(0, dataname.find("\""));
        delete[] this->ScalarArrayName;
        this->ScalarArrayName = new char[dataname.size()+1];
        strcpy(this->ScalarArrayName, dataname.c_str());
        }
      else if (field == "space dimension")
        {
        numDimensions = atoi(description.c_str());
        }
      else if (field == "space origin")
        {
        std::vector<double> origin = ParseVector(description);
        for (size_t i = 0; (i < 3) && (i < origin.size()); i++)
          {
          this->DataOrigin[i] = origin[i];
          }
        }
      else if (field == "space directions")
        {
        std::vector<std::string> directions = split(description);
        dimSpacing.resize(0);
        for (size_t j = 0; j < directions.size(); j++)
          {
          if (directions[j] == "none")
            {
            dimSpacing.push_back(0);
            continue;
            }
          std::vector<double> dir = ParseVector(directions[j]);
          // We don't support orientation, but we do support spacing.
          double mag = 0.0;
          for (size_t k = 0; k < dir.size(); k++)
            {
            mag += dir[k]*dir[k];
            }
          dimSpacing.push_back(sqrt(mag));
          }
        }
      else if (field == "endian")
        {
        if (description == "little")
          {
#ifdef VTK_WORDS_BIGENDIAN
          this->SwapBytes = 1;
#else
          this->SwapBytes = 0;
#endif
          }
        else if (description == "big")
          {
#ifdef VTK_WORDS_BIGENDIAN
          this->SwapBytes = 0;
#else
          this->SwapBytes = 1;
#endif
          }
        else
          {
          vtkErrorMacro(<< "Unknown endian: '" << description << "'");
          return 0;
          }
        }
      else if ((field == "line skip") || (field == "lineskip"))
        {
        if (atoi(description.c_str()) != 0)
          {
          vtkErrorMacro(<< "line skip not supported");
          return 0;
          }
        }
      else if ((field == "byte skip") || (field == "byteskip"))
        {
        if (atoi(description.c_str()) != 0)
          {
          vtkErrorMacro(<< "byte skip not supported");
          return 0;
          }
        }
      else if (   (field == "space units")
               || (field == "sample units") || (field == "sampleunits")
               || (field == "measurement frame")
               || (field == "block size") || (field == "blocksize")
               || (field == "content")
               || (field == "thicknesses")
               || (field == "axis mins") || (field == "axismins")
               || (field == "axis maxs") || (field == "axismaxs")
               || (field == "centers") || (field == "centerings")
               || (field == "units")
               || (field == "kinds")
               || (field == "min") || (field == "max")
               || (field == "old min") || (field == "oldmin")
               || (field == "old max") || (field == "oldmax")
               || (field == "number") )
        {
        // Ignored.
        }
      else
        {
        vtkWarningMacro(<< "Unknown field: '" << field << "'");
        }
      continue;
      }

    delm = line.find(":=");
    if (delm != line.npos)
      {
      // A key/value pair.
      continue;
      }
    }

  // NRRD does not distinguish between vector entries and dimensions.  For
  // example, RGB tuples are represented by adding a dimension of size 3.
  // VTK really needs to know the difference.  Here we are going to guess.
  // If the fastest changing dimension is 9 or less we consider that a tuple.
  // We will also consider any 4th dimension as a tuple.
  if (   (dimSizes.size() > 3)
      || ((dimSizes.size() > 0) && (dimSizes[0] <= 9))
      || ((dimSpacing.size() > 0) && (dimSpacing[0] == 0)) )
    {
    this->NumberOfScalarComponents = dimSizes[0];
    dimSizes.erase(dimSizes.begin());
    if (!dimSpacing.empty()) dimSpacing.erase(dimSpacing.begin());
    subDimension--;
    }
  else
    {
    this->NumberOfScalarComponents = 1;
    }

  // Record the dimensions.
  this->FileDimensionality = static_cast<int>(dimSizes.size());
  for (unsigned int i = 0; i < 3; i++)
    {
    this->DataExtent[i*2+0] = 0;
    this->DataExtent[i*2+1] = (i < dimSizes.size()) ? dimSizes[i]-1 : 0;
    this->DataSpacing[i] = (i < dimSpacing.size()) ? dimSpacing[i] : 1;
    }

  if (this->DataFiles->GetNumberOfValues() > 0)
    {
    if (this->DataFiles->GetNumberOfValues() > 1)
      {
      this->FileDimensionality--;
      if (this->FileDimensionality != 2)
        {
        vtkErrorMacro(<< "Data split into multiple files is only supported"
                      << " when each file is 2D (+ an optional vector"
                      << " dimension).");
        }
      if (subDimension != 3)
        {
        vtkErrorMacro(<< "Data split into multiple files is only supported"
                      << " when each file is 2D (+ an optional vector"
                      << " dimension).  This means the subdim must be on"
                      << " that third (or fourth in the case of a vector)"
                      << " dimension.");
        }
      }
    std::string parentDir
      = vtksys::SystemTools::GetParentDirectory(this->FileName);
    for (vtkIdType i = 0; i < this->DataFiles->GetNumberOfValues(); i++)
      {
      std::string relativePath = this->DataFiles->GetValue(i);
      std::string fullPath
        = vtksys::SystemTools::CollapseFullPath(relativePath.c_str(),
                                                parentDir.c_str());
      this->DataFiles->SetValue(i, fullPath);
      }
    // Header file pointing to data files.  Data files should have no header.
    this->HeaderSize = 0;
    }

  return 1;
}

//=============================================================================
int vtkNrrdReader::RequestData(vtkInformation *request,
                               vtkInformationVector **inputVector,
                               vtkInformationVector *outputVector)
{
  // Get rid of superclasses FileNames.  We don't support that functionality,
  // but we exploit that of the superclass.
  if (this->FileNames != NULL)
    {
    this->FileNames->Delete();
    this->FileNames = NULL;
    }

  char *saveFileName = this->FileName;

  if (this->DataFiles->GetNumberOfValues() == 1)
    {
    this->FileName = const_cast<char *>(this->DataFiles->GetValue(0).c_str());
    }
  else if (this->DataFiles->GetNumberOfValues() > 1)
    {
    this->FileNames = this->DataFiles;
    }

  int result;

  if (this->Encoding == ENCODING_RAW)
    {
    // Superclass knows how to read raw data.  Use that.

    result = this->Superclass::RequestData(request, inputVector, outputVector);
    }
  else if (this->Encoding == ENCODING_ASCII)
    {
    vtkImageData *outputData = vtkImageData::GetData(outputVector);
    this->AllocateOutputData(outputData, outputVector->GetInformationObject(0));
    if (outputData == NULL)
      {
      vtkErrorMacro(<< "Data not created correctly?");
      return 0;
      }

    result = this->ReadDataAscii(outputData);
    }
  else
    {
    vtkErrorMacro(<< "Bad encoding set");
    result = 0;
    }

  this->FileName = saveFileName;
  this->FileNames = NULL;

  return result;
}

//=============================================================================
template<typename T>
int vtkNrrdReaderReadDataAsciiTemplate(vtkNrrdReader *self,
                                       vtkImageData *output,
                                       T *outBuffer)
{
  // Get the requested extent
  int outExtent[6];
  output->GetExtent(outExtent);
  int numComponents = output->GetNumberOfScalarComponents();

  int fileDataExtent[6];
  self->GetDataExtent(fileDataExtent);
  vtkIdType fileDataIncrements[3];
  fileDataIncrements[0] = numComponents;
  fileDataIncrements[1] = fileDataIncrements[0] * fileDataExtent[1];
  fileDataIncrements[2] = fileDataIncrements[1] * fileDataExtent[3];

  vtkStringArray *filenames = self->GetFileNames();
  vtkStdString filename = self->GetFileName();

  std::ifstream file;
  if (self->GetFileDimensionality() == 3)
    {
    if (filenames != NULL)
      {
      filename = filenames->GetValue(0);
      }
    file.open(filename.c_str(), std::ios::in);
    if (file.fail())
      {
      vtkErrorWithObjectMacro(self, "Could not open file " << filename);
      return 0;
      }
    // Skip to start of first slab.
    for (vtkIdType skip = 0; skip < fileDataIncrements[2]*outExtent[4]; skip++)
      {
      T dummy;
      file >> dummy;
      }
    }

  vtkIdType bufferIndex = 0;
  int fileDataIndex[3];
  for (fileDataIndex[2] = outExtent[4];
       fileDataIndex[2] <= outExtent[5];
       fileDataIndex[2]++)
    {
    if (self->GetFileDimensionality() == 2)
      {
      if (file.is_open()) { file.close(); }
      if (filenames != NULL)
        {
        filename = filenames->GetValue(fileDataIndex[2]);
        }
      file.open(filename.c_str(), std::ios::in);
      if (file.fail())
        {
        vtkErrorWithObjectMacro(self, "Could not open file " << filename);
        return 0;
        }
      }
    // Skip unused rows.
    for (vtkIdType skip = 0; skip < fileDataIncrements[1]*outExtent[2]; skip++)
      {
      T dummy;
      file >> dummy;
      }
    for (fileDataIndex[1] = outExtent[2];
         fileDataIndex[1] <= outExtent[3];
         fileDataIndex[1]++)
      {
      // Skip unused columns.
      for (vtkIdType skip = 0; skip<fileDataIncrements[0]*outExtent[0]; skip++)
        {
        T dummy;
        file >> dummy;
        }
      for (fileDataIndex[0] = outExtent[0];
           fileDataIndex[0] <= outExtent[1];
           fileDataIndex[0]++)
        {
        file >> outBuffer[bufferIndex];
        bufferIndex++;
        }
      // Skip unused columns.
      for (vtkIdType skip = 0;
           skip < fileDataIncrements[0]*(fileDataExtent[1] - outExtent[1]);
           skip++)
        {
        T dummy;
        file >> dummy;
        }
      }
    // Skip unused rows.
    for (vtkIdType skip = 0;
         skip < fileDataIncrements[1]*(fileDataExtent[2] - outExtent[2]);
         skip++)
      {
      T dummy;
      file >> dummy;
      }
    }

  file.close();
  return 1;
}

//=============================================================================
int vtkNrrdReader::ReadDataAscii(vtkImageData *output)
{
  void *outBuffer = output->GetScalarPointer();
  switch (output->GetScalarType())
    {
    vtkTemplateMacro(vtkNrrdReaderReadDataAsciiTemplate(this,
                                                        output,
                                                        (VTK_TT *)(outBuffer)));
    default:
      vtkErrorMacro("Unkown data type");
      return 0;
    }

  return 1;
}
