From 51ecd17eaf1df10a01942edfcba656e1c27c1d77 Mon Sep 17 00:00:00 2001
From: Dan Lipsa <dan.lipsa@kitware.com>
Date: Wed, 21 Jun 2017 16:38:37 -0400
Subject: [PATCH 1/2] BUG: Various fixes to the 2D reader

1. Use FormatCode to correctly identify and convert IBM floating point data
2. Convert to vtkStructuredGrid (curvilinear grid). This results in decrease in running time from
   about 50 sec to less than a second.
3. Use the CoordinateMultiplier value to scale X, Y coordinates
4. Use SampleInterval to scale the Z coordinate
5. Fix memory leak
---
 src/BinaryHeaderBytesPositions.h |  12 ++-
 src/IOUtil.cpp                   | 112 +++++++++++++++++++++-
 src/IOUtil.h                     |  11 +--
 src/SegyReader.cpp               | 156 +++++++++++--------------------
 src/SegyReader.h                 |   9 +-
 src/SegyTraceReader.cpp          |  50 ++++++++--
 src/Trace.h                      |   2 +
 src/TraceHeaderBytesPositions.h  |  12 ++-
 src/main.cpp                     |  75 ++++-----------
 src/vtkSegy2DReader.cpp          |  17 +++-
 src/vtkSegy2DReader.h            |   6 +-
 11 files changed, 269 insertions(+), 193 deletions(-)

diff --git a/src/BinaryHeaderBytesPositions.h b/src/BinaryHeaderBytesPositions.h
index e11cc59..f102858 100644
--- a/src/BinaryHeaderBytesPositions.h
+++ b/src/BinaryHeaderBytesPositions.h
@@ -28,9 +28,11 @@ class BinaryHeaderBytesPositions
     int NumSamplesPerTrace;
     int NumSamplesPerTraceOriginal;
     int FormatCode;
-    int NumberExtendedHeaders;
     int EnsembleType;
-    int Version;
+    int MeasurementSystem;
+    int ByteOrderingDetection;  
+    int MajorVersion;
+    int MinorVersion;
     int FixedLengthFlag;
 public:
     BinaryHeaderBytesPositions()
@@ -51,9 +53,11 @@ class BinaryHeaderBytesPositions
         NumSamplesPerTrace = 3220;
         NumSamplesPerTraceOriginal = 3222;
         FormatCode = 3224;
-        NumberExtendedHeaders = 3502;
         EnsembleType = 3228;
-        Version = 3500;
+        MeasurementSystem = 3254;
+        ByteOrderingDetection = 3296;        
+        MajorVersion = 3500;
+        MinorVersion = 3501;
         FixedLengthFlag = 3502;
     }
 };
diff --git a/src/IOUtil.cpp b/src/IOUtil.cpp
index d2d088d..497ff27 100644
--- a/src/IOUtil.cpp
+++ b/src/IOUtil.cpp
@@ -17,13 +17,87 @@
 
 #include <cstring>
 
-IOUtil* IOUtil::instance = NULL;
+/*
+ * Downloaded from:
+ * https://www.thecodingforums.com/threads/c-code-for-converting-ibm-370-floating-point-to-ieee-754.438469/
+ */
+/* ibm2ieee - Converts a number from IBM 370 single precision floating
+point format to IEEE 754 single precision format. For normalized
+numbers, the IBM format has greater range but less precision than the
+IEEE format. Numbers within the overlapping range are converted
+exactly. Numbers which are too large are converted to IEEE Infinity
+with the correct sign. Numbers which are too small are converted to
+IEEE denormalized numbers with a potential loss of precision (including
+complete loss of precision which results in zero with the correct
+sign). When precision is lost, rounding is toward zero (because it's
+fast and easy -- if someone really wants round to nearest it shouldn't
+be TOO difficult). */
+
+#include <sys/types.h>
+#include <netinet/in.h>
+
+void ibm2ieee(void *to, const void *from, int len)
+{
+  register unsigned fr; /* fraction */
+  register int exp; /* exponent */
+  register int sgn; /* sign */
+
+  for (; len-- > 0; to = (char *)to + 4, from = (char *)from + 4) {
+/* split into sign, exponent, and fraction */
+    fr = *(long *)from; /* pick up value */
+    sgn = fr >> 31; /* save sign */
+    fr <<= 1; /* shift sign out */
+    exp = fr >> 25; /* save exponent */
+    fr <<= 7; /* shift exponent out */
+
+    if (fr == 0) { /* short-circuit for zero */
+      exp = 0;
+      goto done;
+    }
+
+/* adjust exponent from base 16 offset 64 radix point before first digit
+   to base 2 offset 127 radix point after first digit */
+/* (exp - 64) * 4 + 127 - 1 == exp * 4 - 256 + 126 == (exp << 2) - 130 */
+    exp = (exp << 2) - 130;
+
+/* (re)normalize */
+    while (fr < 0x80000000) { /* 3 times max for normalized input */
+      --exp;
+      fr <<= 1;
+    }
+
+    if (exp <= 0) { /* underflow */
+      if (exp < -24) /* complete underflow - return properly signed zero */
+        fr = 0;
+      else /* partial underflow - return denormalized number */
+        fr >>= -exp;
+      exp = 0;
+    } else if (exp >= 255) { /* overflow - return infinity */
+      fr = 0;
+      exp = 255;
+    } else { /* just a plain old number - remove the assumed high bit */
+      fr <<= 1;
+    }
+
+  done:
+/* put the pieces back together and return it */
+    *(unsigned *)to = (fr >> 9) | (exp << 23) | (sgn << 31);
+  }
+}
+
 
 IOUtil::IOUtil()
 {
     isBigEndian = checkIfBigEndian();
 }
 
+IOUtil* IOUtil::Instance()
+{
+  static IOUtil ioUtil;
+  return &ioUtil;
+}
+
+
 int IOUtil::readShortInteger(int pos, ifstream &in) {
     in.seekg(pos, in.beg);
     char buffer[2];
@@ -53,6 +127,21 @@ int IOUtil::readLongInteger(int pos, ifstream &in) {
     return num;
 }
 
+int IOUtil::readLongInteger(ifstream &in) {
+    char buffer[4];
+    in.read(buffer, sizeof(buffer));
+
+    if (!isBigEndian) {
+        swap(buffer, buffer + 3);
+        swap(buffer + 1, buffer + 2);
+    }
+
+    int num;
+    memcpy(&num, buffer, 4);
+    return num;
+}
+
+
 float IOUtil::readFloat(ifstream &in) {
     char buffer[4];
     in.read(buffer, sizeof(buffer));
@@ -67,12 +156,33 @@ float IOUtil::readFloat(ifstream &in) {
     return num;
 }
 
+float IOUtil::readIBMFloat(ifstream &in) {
+    char buffer[4];
+    in.read(buffer, sizeof(buffer));
+
+    if (!isBigEndian) {
+        swap(buffer, buffer + 3);
+        swap(buffer + 1, buffer + 2);
+    }
+
+    float num;
+    ibm2ieee(&num, buffer, 1);
+    return num;
+}
+
+
 char IOUtil::readChar(ifstream &in) {
     char buffer;
     in.read(&buffer, sizeof(buffer));
     return buffer;
 }
 
+unsigned char IOUtil::readUChar(ifstream &in) {
+    char buffer;
+    in.read(&buffer, sizeof(buffer));
+    return buffer;
+}
+
 void IOUtil::swap(char *a, char *b) {
     char temp = *a;
     *a = *b;
diff --git a/src/IOUtil.h b/src/IOUtil.h
index 294cc88..f09b4b7 100644
--- a/src/IOUtil.h
+++ b/src/IOUtil.h
@@ -23,21 +23,18 @@ class IOUtil {
 public:
     int readShortInteger(int pos, ifstream &in);
     int readLongInteger(int pos, ifstream &in);
+    int readLongInteger(ifstream &in);  
     float readFloat(ifstream &in);
+    float readIBMFloat(ifstream &in);
     char readChar(ifstream &in);
+    unsigned char readUChar(ifstream &in);  
     void swap(char* a, char* b);
     bool isBigEndian;
-    static IOUtil* Instance()
-    {
-        if(instance == NULL)
-            instance = new IOUtil();
-        return instance;
-    }
+    static IOUtil* Instance();
     int getFileSize(ifstream& in);
 
 private:
     IOUtil();
-    static IOUtil* instance;
     bool checkIfBigEndian() {
         ushort a = 0x1234;
         if (*((unsigned char *) &a) == 0x12)
diff --git a/src/SegyReader.cpp b/src/SegyReader.cpp
index 1d60aa4..c13cf30 100644
--- a/src/SegyReader.cpp
+++ b/src/SegyReader.cpp
@@ -21,12 +21,8 @@
 #include <vtkSmartPointer.h>
 #include <vtkPolygon.h>
 #include <vtkStructuredPointsReader.h>
-#include <vtkVolumeTextureMapper3D.h>
-#include <vtkColorTransferFunction.h>
 #include <vtkRenderWindow.h>
-#include <vtkAxesActor.h>
-#include <vtkImageShiftScale.h>
-#include <vtkPolyDataMapper.h>
+#include <vtkStructuredGrid.h>
 #include <vtkFloatArray.h>
 #include <vtkPointData.h>
 #include <vtkArrayData.h>
@@ -36,12 +32,11 @@
 #include <vtkInformationVector.h>
 #include <vtkInformation.h>
 #include <vtkCellData.h>
-//#include <vtkXMLPolyDataWriter.h>
 
 
 SegyReader::~SegyReader()
 {
-    for(auto trace : data)
+    for(auto trace : traces)
         delete trace;
 }
 
@@ -61,8 +56,11 @@ bool SegyReader::LoadFromFile(string path) {
     {
         Trace* pTrace = new Trace();
         if (!traceReader.readTrace(traceStartPos, in, formatCode, pTrace))
+        {
+            delete pTrace;
             break;
-        data.push_back(pTrace);
+        }
+        traces.push_back(pTrace);
     }
 
     in.close();
@@ -71,8 +69,25 @@ bool SegyReader::LoadFromFile(string path) {
 
 
 bool SegyReader::readHeader(ifstream &in) {
+    short sampleInterval = IOUtil::Instance()->readShortInteger(binaryHeaderBytesPos.SampleInterval, in);
     formatCode = IOUtil::Instance()->readShortInteger(binaryHeaderBytesPos.FormatCode, in);
+    in.seekg(binaryHeaderBytesPos.MajorVersion, in.beg);
+    unsigned char majorVersion = IOUtil::Instance()->readUChar(in);
+    unsigned char minorVersion = IOUtil::Instance()->readUChar(in);
     sampleCountPerTrace =  IOUtil::Instance()->readShortInteger(binaryHeaderBytesPos.NumSamplesPerTrace, in);
+    short tracesPerEnsemble = IOUtil::Instance()->readShortInteger(binaryHeaderBytesPos.NumberTracesPerEnsemble, in);
+    short ensembleType = IOUtil::Instance()->readShortInteger(binaryHeaderBytesPos.EnsembleType, in);
+    short measurementSystem = IOUtil::Instance()->readShortInteger(binaryHeaderBytesPos.MeasurementSystem, in);    
+    int byteOrderingDetection = IOUtil::Instance()->readLongInteger(binaryHeaderBytesPos.ByteOrderingDetection, in);
+
+    std::cout << "Segy version = " << int(majorVersion) << "." << int(minorVersion) << std::endl;    
+    std::cout << "FormatCode = " << formatCode << std::endl;    
+    std::cout << "ByteOrderingDetection = " << byteOrderingDetection << std::endl;
+    std::cout << "SampleCountPerTrace=" << sampleCountPerTrace << std::endl;
+    std::cout << "ensembleType=" << ensembleType << std::endl;
+    std::cout << "measurementSystem=" << measurementSystem << std::endl;
+    std::cout << "sampleInterval=" << sampleInterval << std::endl;
+    std::cout << "tracesPerEnsemble=" << tracesPerEnsemble << std::endl << std::endl;
     return true;
 }
 
@@ -81,7 +96,7 @@ bool SegyReader::readHeader(ifstream &in) {
 bool SegyReader::ExportData3D(vtkImageData *imageData)
 {
   set<int> crosslineNumbers, inlineNumbers;
-  for(auto trace : data)
+  for(auto trace : traces)
   {
       crosslineNumbers.insert(trace->crosslineNumber);
       inlineNumbers.insert(trace->inlineNumber);
@@ -97,7 +112,7 @@ bool SegyReader::ExportData3D(vtkImageData *imageData)
   float min_data = INT_MAX;
   float max_data = INT_MIN;
 
-  for(auto trace : data)
+  for(auto trace : traces)
   {
       int cross = trace->crosslineNumber;
       auto pair = cross_inline_map.find(cross);
@@ -163,7 +178,7 @@ bool SegyReader::GetImageData(vtkImageData *imageData)
   int minCrossLineNumber = INT_MAX;
   int maxCrossLineNumber = INT_MIN;
 
-  for(auto trace : data)
+  for(auto trace : traces)
   {
       crosslineNum = trace->crosslineNumber;
       if (crosslineNum == 0)
@@ -184,7 +199,7 @@ bool SegyReader::GetImageData(vtkImageData *imageData)
   float min_data = INT_MAX;
   float max_data = INT_MIN;
 
-  for(auto trace : data)
+  for(auto trace : traces)
   {
       for(auto m : trace->data)
       {
@@ -211,7 +226,7 @@ bool SegyReader::GetImageData(vtkImageData *imageData)
         index = 0;
         remainder = 0;
 
-        for(auto trace : data)
+        for(auto trace : traces)
         {
           dataSize = trace->data.size();
           diff = aggIndex - dataSize;
@@ -228,7 +243,7 @@ bool SegyReader::GetImageData(vtkImageData *imageData)
           }
         }
 
-        *(ptr + i * sampleCountPerTrace + k) = 256.0 * (data[index]->data[remainder] - min_data)/(max_data - min_data);
+        *(ptr + i * sampleCountPerTrace + k) = 256.0 * (traces[index]->data[remainder] - min_data)/(max_data - min_data);
       }
   }
 
@@ -236,116 +251,55 @@ bool SegyReader::GetImageData(vtkImageData *imageData)
 }
 
 
-bool SegyReader::AddScalars(vtkPolyData* polyData)
+void SegyReader::AddScalars(vtkStructuredGrid* grid)
 {
 
-  vtkSmartPointer<vtkFloatArray> cellData =
+  vtkSmartPointer<vtkFloatArray> pointData =
     vtkSmartPointer<vtkFloatArray>::New();
-  cellData->SetName("trace");
-  cellData->SetNumberOfComponents(1);
-
-  int crossLineCount = data.size();
+  pointData->SetName("trace");
+  pointData->SetNumberOfComponents(1);
 
-  cellData->Allocate(crossLineCount * sampleCountPerTrace);
+  int crossLineCount = traces.size();
 
-  float min_data = INT_MAX;
-  float max_data = INT_MIN;
+  pointData->Allocate(crossLineCount * sampleCountPerTrace);
 
-  for(auto trace : data)
+  for (int k = 0; k < sampleCountPerTrace; k++)
   {
-      for(auto m : trace->data)
-      {
-          if(m < min_data)
-              min_data = m;
-          if(m > max_data )
-              max_data = m;
-      }
+    for (int i = 0; i < traces.size(); i++)
+    {
+      pointData->InsertValue( i * sampleCountPerTrace + k, traces[i]->data[k]);
+    }
   }
 
-      for (int k = 0; k < sampleCountPerTrace; k++)
-      {
-        for (int i = 0; i < data.size(); i++)
-        {
-        cellData->InsertValue( i * sampleCountPerTrace + k, 256.0 * (data[i]->data[k] - min_data)/(max_data-min_data) );
-        }
-      }
-
-    polyData->GetCellData()->SetScalars(cellData);
-    polyData->GetCellData()->SetActiveScalars("trace");
-
-
-    return true;
+  grid->GetPointData()->SetScalars(pointData);
+  grid->GetPointData()->SetActiveScalars("trace");
 }
 
 
-bool SegyReader::ExportData2D(vtkPolyData * polyData)
+void SegyReader::ExportData2D(vtkStructuredGrid * grid)
 {
+    grid->SetDimensions(traces.size(), sampleCountPerTrace, 1);
     vtkSmartPointer<vtkPoints> points =
             vtkSmartPointer<vtkPoints>::New();
 
-    vtkSmartPointer<vtkFloatArray> textureCoordinates =
-            vtkSmartPointer<vtkFloatArray>::New();
-    textureCoordinates->SetNumberOfComponents(2);
-    textureCoordinates->SetName("TextureCoordinates");
-
     for(int k=0; k<sampleCountPerTrace; k++)
     {
-        for(int i=0;i < data.size(); i++)
+      for(int i=0;i < traces.size(); i++)    
         {
-            auto trace = data[i];
-            float x = trace->xCoordinate / 100000.0;
-            float y = trace->yCoordinate / 100000.0;
-
-            float z = k * 100.0 / sampleCountPerTrace;
+            auto trace = traces[i];
+            float coordinateMultiplier = (trace->CoordinateMultiplier < 0) ?
+              1.0 / (- trace->CoordinateMultiplier) :
+              trace->CoordinateMultiplier;
+            float x = trace->xCoordinate * coordinateMultiplier;
+            float y = trace->yCoordinate * coordinateMultiplier;
+
+            float z = k * trace->SampleInterval / (sampleCountPerTrace - 1);
             points->InsertNextPoint(x, y, z);
-            textureCoordinates->InsertNextTuple2(k * 1.0 / sampleCountPerTrace, i * 1.0 / data.size());
         }
     }
 
-// Create a cell array to store the quad in
-    vtkSmartPointer<vtkCellArray> quads =
-            vtkSmartPointer<vtkCellArray>::New();
-
-    for (int k = 1; k < sampleCountPerTrace; k++)
-    {
-        for(int i=1;i < data.size(); i++)
-        {
-            auto trace = data[i];
-
-            vtkSmartPointer<vtkPolygon> polygon =
-                    vtkSmartPointer<vtkPolygon>::New();
-            polygon->GetPointIds()->SetNumberOfIds(4); //make a quad
-
-            int id1 = k * data.size() + i;
-            int id2 = (k - 1) * data.size() + i;
-            int id3 = (k - 1) * data.size() + i - 1;
-            int id4 = k * data.size() + i - 1;
-            polygon->GetPointIds()->SetId(0, id1);
-            polygon->GetPointIds()->SetId(1, id2);
-            polygon->GetPointIds()->SetId(2, id3);
-            polygon->GetPointIds()->SetId(3, id4);
-            quads->InsertNextCell(polygon);
-        }
-    }
-
-    // vtkSmartPointer<vtkTexture> texture =
-    //        vtkSmartPointer<vtkTexture>::New();
-
-    // vtkSmartPointer<vtkImageData> imageData = vtkSmartPointer<vtkImageData>::New();
-    // ExportData3D(imageData);
-    // std::cerr << "imageData " << imageData << std::endl;
-    // texture->SetInputDataObject(imageData);
-
-    polyData->SetPoints(points);
-    polyData->SetPolys(quads);
-    this->AddScalars(polyData);
-     polyData->GetPointData()->SetTCoords(textureCoordinates);
-
-//    vtkNew<vtkXMLPolyDataWriter> writer;
-//    writer->SetInputData(polyData);
-//    writer->SetFileName("test.vtp");
-//    writer->Write();
-    return polyData;
+    grid->SetPoints(points);
+    this->AddScalars(grid);
 }
 
 
diff --git a/src/SegyReader.h b/src/SegyReader.h
index 29181f5..c5abd28 100644
--- a/src/SegyReader.h
+++ b/src/SegyReader.h
@@ -24,8 +24,6 @@
 
 #include <vtkImageData.h>
 #include <vtkActor.h>
-#include <vtkPolyData.h>
-#include <vtkPolyDataAlgorithm.h>
 
 #include <iostream>
 #include <fstream>
@@ -33,6 +31,7 @@
 #include <map>
 
 using namespace std;
+class vtkStructuredGrid;
 
 class SegyReader
 {
@@ -45,11 +44,11 @@ class SegyReader
     bool GetImageData(vtkImageData *); // export the data in segy file as 3D volume
     bool ExportData3D(vtkImageData *); // export the data in segy file as 3D volume
     bool LoadFromFile(string path);
-    bool ExportData2D(vtkPolyData *);
-    bool AddScalars(vtkPolyData *);
+    void ExportData2D(vtkStructuredGrid *);
+    void AddScalars(vtkStructuredGrid *);
 private:
     bool readHeader(ifstream& in);
-    vector<Trace*> data;
+    vector<Trace*> traces;
     int formatCode;
     BinaryHeaderBytesPositions binaryHeaderBytesPos;
     SegyTraceReader traceReader;
diff --git a/src/SegyTraceReader.cpp b/src/SegyTraceReader.cpp
index 4f88807..e7cad7c 100644
--- a/src/SegyTraceReader.cpp
+++ b/src/SegyTraceReader.cpp
@@ -23,24 +23,41 @@ void SegyTraceReader::printTraceHeader(ifstream &in, int startPos) {
 //    cout << "Position:" << startPos << endl;
 
     int traceSequenceNumberInLine = IOUtil::Instance()->readLongInteger(startPos + traceHeaderBytesPos.TraceNumber, in);
-//    cout << "Trace sequence number in line : " << traceSequenceNumberInLine << endl;
+    cout << "Trace sequence number in line : " << traceSequenceNumberInLine << endl;
+
+    int traceSequenceNumberInFile = IOUtil::Instance()->readLongInteger(in);
+    cout << "Trace sequence number in file : " << traceSequenceNumberInFile << endl;
+    
 
   // Get number_of_samples from trace header position 115-116
     int numSamples = IOUtil::Instance()->readShortInteger(startPos + traceHeaderBytesPos.NumberSamples, in);
-//    cout << "number of samples: " << numSamples << endl;
+    cout << "number of samples: " << numSamples << endl;
+
+    short sampleInterval = IOUtil::Instance()->readShortInteger(startPos + traceHeaderBytesPos.SampleInterval, in);
+    cout << "sample interval: " << sampleInterval << endl;
+
 
 // Get inline number from trace header position 189-192
     int inlineNum = IOUtil::Instance()->readLongInteger(startPos + traceHeaderBytesPos.InlineNumber, in);
-//    cout << "in-line number : " << inlineNum << endl;
+    cout << "Field record number (inline number) : " << inlineNum << endl;
 
     int crosslineNum = IOUtil::Instance()->readLongInteger(startPos + traceHeaderBytesPos.CrosslineNumber, in);
-//    cout << "cross-line number : " << crosslineNum << endl;
+    cout << "cross-line number (ensemble number) : " << crosslineNum << endl;
 
+    int traceNumberWithinEnsemble = IOUtil::Instance()->readLongInteger(startPos + traceHeaderBytesPos.TraceNumberWithinEnsemble, in);
+    cout << "trace number within ensemble : " << traceNumberWithinEnsemble << endl;
+
+    short coordinateMultiplier = IOUtil::Instance()->readShortInteger(startPos + traceHeaderBytesPos.CoordinateMultiplier, in);
+    cout << "coordinate multiplier : " << coordinateMultiplier << endl;
+    
     int xCoordinate = IOUtil::Instance()->readLongInteger(startPos + traceHeaderBytesPos.XCoordinate, in);
-//    cout << "X coordinate for ensemble position of the trace : " << xCoordinate << endl;
+    cout << "X coordinate for ensemble position of the trace : " << xCoordinate << endl;
 
     int yCoordinate = IOUtil::Instance()->readLongInteger(startPos + traceHeaderBytesPos.YCoordinate, in);
-//    cout << "Y coordinate for ensemble position of the trace : " << yCoordinate << endl;
+    cout << "Y coordinate for ensemble position of the trace : " << yCoordinate << endl;
+
+    short coordinateUnits = IOUtil::Instance()->readShortInteger(startPos + traceHeaderBytesPos.CoordinateUnits, in);
+    cout << "coordinateUnits: " << coordinateUnits << endl;
 }
 
 bool SegyTraceReader::readTrace(int &startPos, ifstream &in, int formatCode, Trace* trace) {
@@ -55,16 +72,29 @@ bool SegyTraceReader::readTrace(int &startPos, ifstream &in, int formatCode, Tra
     int numSamples = IOUtil::Instance()->readShortInteger(startPos + traceHeaderBytesPos.NumberSamples, in);
     trace->xCoordinate = IOUtil::Instance()->readLongInteger(startPos + traceHeaderBytesPos.XCoordinate, in);
     trace->yCoordinate = IOUtil::Instance()->readLongInteger(startPos + traceHeaderBytesPos.YCoordinate, in);
+    trace->CoordinateMultiplier = IOUtil::Instance()->readShortInteger(startPos + traceHeaderBytesPos.CoordinateMultiplier, in);
+    trace->SampleInterval = IOUtil::Instance()->readShortInteger(startPos + traceHeaderBytesPos.SampleInterval, in);
 
     in.seekg(startPos + 240, in.beg);
+    //std::cout << "Values: ";
     for (int i = 0; i < numSamples; i++) {
         float value;
-        if(formatCode == 8)
-            value = IOUtil::Instance()->readChar(in);
-        else
-            value = IOUtil::Instance()->readFloat(in); // TODO: readChar or readFloat according to format code
+        switch (formatCode)
+        {
+        case 1:
+          value = IOUtil::Instance()->readIBMFloat(in);
+          break;
+        case 5:
+          value = IOUtil::Instance()->readFloat(in);
+          break;
+        default:
+          std::cerr << "Data sample format code " << formatCode << " not supported." << std::endl;
+          value = 0;
+        }
         trace->data.push_back(value);
+        //std::cout << value << " ";
     }
+    cout << endl;
 
     startPos += 240 + getTraceSize(numSamples, formatCode);
     return true;
diff --git a/src/Trace.h b/src/Trace.h
index dbbead0..16ab66f 100644
--- a/src/Trace.h
+++ b/src/Trace.h
@@ -23,9 +23,11 @@ class Trace {
 public:
     float xCoordinate;
     float yCoordinate;
+    short CoordinateMultiplier;
     vector<float> data;
     int inlineNumber;
     int crosslineNumber;
+    short SampleInterval;
 };
 
 
diff --git a/src/TraceHeaderBytesPositions.h b/src/TraceHeaderBytesPositions.h
index 4c7379f..5d68837 100644
--- a/src/TraceHeaderBytesPositions.h
+++ b/src/TraceHeaderBytesPositions.h
@@ -20,10 +20,14 @@ class TraceHeaderBytesPositions
 public:
     int TraceNumber;
     int InlineNumber;
-    int NumberSamples;
     int CrosslineNumber;
+    int TraceNumberWithinEnsemble;
+    int CoordinateMultiplier;
     int XCoordinate;
     int YCoordinate;
+    int CoordinateUnits;
+    int NumberSamples;
+    int SampleInterval;
 public:
     TraceHeaderBytesPositions()
     {
@@ -34,10 +38,14 @@ class TraceHeaderBytesPositions
     {
         TraceNumber = 0;
         InlineNumber = 8;
-        NumberSamples = 114;
         CrosslineNumber = 20;
+        TraceNumberWithinEnsemble = 24;
+        CoordinateMultiplier = 70;
         XCoordinate = 72;
         YCoordinate = 76;
+        CoordinateUnits = 88;
+        NumberSamples = 114;
+        SampleInterval = 116;
     }
 };
 #endif //SEGYVISUALIZER_TRACEHEADERBYTESPOSITIONS_H
diff --git a/src/main.cpp b/src/main.cpp
index 7bb613b..a5274de 100644
--- a/src/main.cpp
+++ b/src/main.cpp
@@ -33,6 +33,8 @@ using namespace std;
 #include <vtkOrientationMarkerWidget.h>
 #include <vtkPolyDataMapper.h>
 #include <vtkPolyData.h>
+#include <vtkXMLStructuredGridWriter.h>
+#include "vtksys/SystemTools.hxx"
 
 #include <vtkPlane.h>
 #include <vtkCellArray.h>
@@ -205,66 +207,29 @@ void render(vtkImageData* id)
 
 void demo2D()
 {
-    vtkSmartPointer<vtkColorTransferFunction> colorTransferFunction =
-      vtkSmartPointer<vtkColorTransferFunction>::New();
-    colorTransferFunction->AddRGBPoint(0.0, 0.0, 0.0, 1.0);
-    colorTransferFunction->AddRGBPoint(64.0, 1.0, 0.0, 0.0);
-    colorTransferFunction->AddRGBPoint(128.0, 0.0, 0.0, 1.0);
-    colorTransferFunction->AddRGBPoint(192.0, 0.0, 1.0, 0.0);
-    colorTransferFunction->AddRGBPoint(255.0, 0.0, 0.2, 0.0);
-
-    vtkSmartPointer<vtkRenderer> renderer =
-            vtkSmartPointer<vtkRenderer>::New();
-    vtkSmartPointer<vtkRenderWindow> renderWindow =
-            vtkSmartPointer<vtkRenderWindow>::New();
-    renderWindow->AddRenderer(renderer);
-    vtkSmartPointer<vtkRenderWindowInteractor> renderWindowInteractor =
-            vtkSmartPointer<vtkRenderWindowInteractor>::New();
-    renderWindowInteractor->SetRenderWindow(renderWindow);
-
     vector<string> files;
     files.push_back("data/lineA.sgy");
-//    files.push_back("data/lineB.sgy");
-//    files.push_back("data/lineC.sgy");
-//    files.push_back("data/lineD.sgy");
-//    files.push_back("data/lineE.sgy");
+    files.push_back("data/lineB.sgy");
+    files.push_back("data/lineC.sgy");
+    files.push_back("data/lineD.sgy");
+    files.push_back("data/lineE.sgy");
 
-    auto file = files[0];
-    //for(auto file : files)
+    for (auto file: files)
     {
-        vtkSmartPointer<vtkSegy2DReader> reader =
-          vtkSmartPointer<vtkSegy2DReader>::New();
-        reader->SetFileName(file.c_str());
-        reader->Update();
-
-        vtkSmartPointer<vtkActor> actor = vtkSmartPointer<vtkActor>::New();
-
-        vtkSmartPointer<vtkPolyDataMapper> mapper =
-                vtkSmartPointer<vtkPolyDataMapper>::New();
-
-
-        mapper->SetInputConnection(reader->GetOutputPort());
-        mapper->ScalarVisibilityOff();
-        actor->SetMapper(mapper);
-
-        vtkNew<vtkImageData> imageData;
-        reader->GetImageData(imageData.GetPointer());
-        vtkNew<vtkTexture> texture;
-        texture->SetInputData(imageData.GetPointer());
-        actor->SetTexture(texture.GetPointer());
-
-        renderer->AddActor(actor);
-
-        vtkCamera* camera = renderer->GetActiveCamera();
-        camera->SetViewUp(0, 0, 1.);
-        camera->SetPosition(1, 0, 0);
-        camera->SetFocalPoint(-1, 0, 0);
-        renderer->ResetCamera();
+      std::cout << "Writing to file ..." << std::endl;
+      std::string name = vtksys::SystemTools::GetFilenameName(file);
+      std::string nameNoExt = vtksys::SystemTools::GetFilenameWithoutExtension(name);
+      
+      vtkSmartPointer<vtkSegy2DReader> reader =
+        vtkSmartPointer<vtkSegy2DReader>::New();
+      reader->SetFileName(file.c_str());
+      reader->Update();
+
+      vtkNew<vtkXMLStructuredGridWriter> writer;
+      writer->SetInputData(reader->GetOutput());
+      writer->SetFileName((nameNoExt + ".vts").c_str());
+      writer->Write();
     }
-
-
-    renderWindow->Render();
-    renderWindowInteractor->Start();
 }
 
 void demo3D()
diff --git a/src/vtkSegy2DReader.cpp b/src/vtkSegy2DReader.cpp
index 4a8e871..6f31fad 100644
--- a/src/vtkSegy2DReader.cpp
+++ b/src/vtkSegy2DReader.cpp
@@ -19,17 +19,19 @@
 #include <vtkInformationVector.h>
 #include <vtkObjectFactory.h>
 
+#include <chrono>
+
 vtkStandardNewMacro(vtkSegy2DReader);
 
 vtkSegy2DReader::vtkSegy2DReader()
 {
-    this->SetFileName(0);
+    this->FileName = nullptr;
     this->SetNumberOfInputPorts( 0 );
 }
 
 vtkSegy2DReader::~vtkSegy2DReader()
 {
-    this->SetFileName(0);
+    this->SetFileName(nullptr);
 }
 
 int vtkSegy2DReader::RequestData(vtkInformation* vtkNotUsed(request),
@@ -37,7 +39,7 @@ int vtkSegy2DReader::RequestData(vtkInformation* vtkNotUsed(request),
                                   vtkInformationVector* outputVector)
 {
     vtkInformation* outInfo = outputVector->GetInformationObject(0);
-    vtkPolyData* output = vtkPolyData::SafeDownCast(outInfo->Get(vtkDataObject::DATA_OBJECT()));
+    vtkStructuredGrid* output = vtkStructuredGrid::SafeDownCast(outInfo->Get(vtkDataObject::DATA_OBJECT()));
 
     if (!this->FileName)
     {
@@ -46,8 +48,13 @@ int vtkSegy2DReader::RequestData(vtkInformation* vtkNotUsed(request),
     }
 
     reader.LoadFromFile(FileName);
-    if(!reader.ExportData2D(output))
-        cout << "Failed to request data for vtkSegy2DReader" << endl;
+    std::chrono::time_point<std::chrono::system_clock> start, end;
+    start = std::chrono::system_clock::now();
+    std::cout << "Exporting to poly data ..." << std::endl;
+    reader.ExportData2D(output);
+    end = std::chrono::system_clock::now();
+    std::chrono::duration<double> elapsed_seconds = end-start;
+    std::cout << "Elapsed time: " << elapsed_seconds.count() << std::endl;
 }
 
 void vtkSegy2DReader::PrintSelf(ostream &os, vtkIndent indent)
diff --git a/src/vtkSegy2DReader.h b/src/vtkSegy2DReader.h
index baa8247..5aa4c44 100644
--- a/src/vtkSegy2DReader.h
+++ b/src/vtkSegy2DReader.h
@@ -18,14 +18,14 @@
 
 #include "vtkSegy2DReader.h"
 
-#include <vtkPolyDataAlgorithm.h>
+#include <vtkStructuredGridAlgorithm.h>
 #include "SegyReader.h"
 
-class vtkSegy2DReader : public vtkPolyDataAlgorithm
+class vtkSegy2DReader : public vtkStructuredGridAlgorithm
 {
 public:
     static vtkSegy2DReader* New();
-    vtkTypeMacro(vtkSegy2DReader,vtkPolyDataAlgorithm);
+    vtkTypeMacro(vtkSegy2DReader,vtkStructuredGridAlgorithm);
     void PrintSelf(ostream& os, vtkIndent indent);
 
     vtkSetStringMacro(FileName);

From a41acc9e8ed9e1975f2bcbb56f08930ef2222f0f Mon Sep 17 00:00:00 2001
From: Dan Lipsa <dan.lipsa@kitware.com>
Date: Mon, 26 Jun 2017 12:04:56 -0400
Subject: [PATCH 2/2] Fix the ordering of scalars.

---
 src/SegyReader.cpp | 5 +++--
 1 file changed, 3 insertions(+), 2 deletions(-)

diff --git a/src/SegyReader.cpp b/src/SegyReader.cpp
index c13cf30..a6c4e0b 100644
--- a/src/SegyReader.cpp
+++ b/src/SegyReader.cpp
@@ -263,11 +263,12 @@ void SegyReader::AddScalars(vtkStructuredGrid* grid)
 
   pointData->Allocate(crossLineCount * sampleCountPerTrace);
 
+  int j = 0;
   for (int k = 0; k < sampleCountPerTrace; k++)
   {
     for (int i = 0; i < traces.size(); i++)
     {
-      pointData->InsertValue( i * sampleCountPerTrace + k, traces[i]->data[k]);
+      pointData->InsertValue(j++, traces[i]->data[k]);
     }
   }
 
@@ -284,7 +285,7 @@ void SegyReader::ExportData2D(vtkStructuredGrid * grid)
 
     for(int k=0; k<sampleCountPerTrace; k++)
     {
-      for(int i=0;i < traces.size(); i++)    
+        for(int i=0;i < traces.size(); i++)
         {
             auto trace = traces[i];
             float coordinateMultiplier = (trace->CoordinateMultiplier < 0) ?