SEGYImport.cpp

/****************************************************************************
** Copyright 2019 The Open Group
** Copyright 2019 Bluware, Inc.
**
** Licensed under the Apache License, Version 2.0 (the "License");
** you may not use this file except in compliance with the License.
** You may obtain a copy of the License at
**
**   http://www.apache.org/licenses/LICENSE-2.0
**
** Unless required by applicable law or agreed to in writing, software
** distributed under the License is distributed on an "AS IS" BASIS,
** WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
** See the License for the specific language governing permissions and
** limitations under the License.
****************************************************************************/

#define _CRT_SECURE_NO_WARNINGS 1
#include <SEGYUtils/SEGYFileInfo.h>
#include "IO/File.h"
#include "VDS/Hash.h"
#include <SEGYUtils/DataProvider.h>
#include <SEGYUtils/TraceDataManager.h>

#include <OpenVDS/OpenVDS.h>
#include <OpenVDS/MetadataContainer.h>
#include <OpenVDS/VolumeDataLayoutDescriptor.h>
#include <OpenVDS/VolumeDataAxisDescriptor.h>
#include <OpenVDS/VolumeDataChannelDescriptor.h>
#include <OpenVDS/VolumeDataAccess.h>
#include <OpenVDS/Range.h>
#include <OpenVDS/VolumeDataLayout.h>
#include <OpenVDS/KnownMetadata.h>
#include <OpenVDS/GlobalMetadataCommon.h>

#include <mutex>
#include <cstdlib>
#include <climits>
#include <cassert>
#include <algorithm>

#include "cxxopts.hpp"
#include <json/json.h>
#include <fmt/format.h>

#include <chrono>
#include <numeric>

#if defined(WIN32)
#undef WIN32_LEAN_AND_MEAN // avoid warnings if defined on command line
#define WIN32_LEAN_AND_MEAN 1
#define NOMINMAX 1
#include <io.h>
#include <windows.h>

int64_t GetTotalSystemMemory()
{
    MEMORYSTATUSEX status;
    status.dwLength = sizeof(status);
    GlobalMemoryStatusEx(&status);
    return int64_t(status.ullTotalPhys);
}

#else
#include <unistd.h>

int64_t GetTotalSystemMemory()
{
    long pages = sysconf(_SC_PHYS_PAGES);
    long page_size = sysconf(_SC_PAGE_SIZE);
    return int64_t(pages) * int64_t(page_size);
}
#endif

static DataProvider CreateDataProviderFromFile(const std::string &filename, OpenVDS::Error &error)
{
  std::unique_ptr<OpenVDS::File> file(new OpenVDS::File());
  if (!file->Open(filename, false, false, false, error))
    return DataProvider((OpenVDS::File *)nullptr);
  return DataProvider(file.release());
}

static DataProvider CreateDataProviderFromOpenOptions(const std::string &url, const std::string &connectionString, const std::string &objectId, OpenVDS::Error &error)
{
  std::unique_ptr<OpenVDS::IOManager> ioManager(OpenVDS::IOManager::CreateIOManager(url, connectionString, error));
  if (error.code)
    return DataProvider((OpenVDS::IOManager *)nullptr, "", error);
  return DataProvider(ioManager.release(), objectId, error);
}

Json::Value
SerializeSEGYBinInfo(SEGYBinInfo const& binInfo)
{
  Json::Value
    jsonBinInfo;

  jsonBinInfo["inlineNumber"] = binInfo.m_inlineNumber;
  jsonBinInfo["crosslineNumber"] = binInfo.m_crosslineNumber;

  jsonBinInfo["ensembleXCoordinate"] = binInfo.m_ensembleXCoordinate;
  jsonBinInfo["ensembleYCoordinate"] = binInfo.m_ensembleYCoordinate;

  return jsonBinInfo;
}

Json::Value
SerializeSEGYSegmentInfo(SEGYSegmentInfo const& segmentInfo)
{
  Json::Value
    jsonSegmentInfo;

  jsonSegmentInfo["primaryKey"] = segmentInfo.m_primaryKey;
  jsonSegmentInfo["traceStart"] = segmentInfo.m_traceStart;
  jsonSegmentInfo["traceStop"] = segmentInfo.m_traceStop;

  jsonSegmentInfo["binInfoStart"] = SerializeSEGYBinInfo(segmentInfo.m_binInfoStart);
  jsonSegmentInfo["binInfoStop"] = SerializeSEGYBinInfo(segmentInfo.m_binInfoStop);

  return jsonSegmentInfo;
}

std::string
ToString(SEGY::Endianness endiannness)
{
  switch (endiannness)
  {
  case SEGY::Endianness::BigEndian:  return "BigEndian";
  case SEGY::Endianness::LittleEndian: return "LittleEndian";
  default:
    assert(0); return "";
  }
}

std::string
ToString(SEGY::FieldWidth fieldWidth)
{
  switch (fieldWidth)
  {
  case SEGY::FieldWidth::TwoByte:  return "TwoByte";
  case SEGY::FieldWidth::FourByte: return "FourByte";
  default:
    assert(0); return "";
  }
}

Json::Value
SerializeSEGYHeaderField(SEGY::HeaderField const& headerField)
{
  Json::Value
    jsonHeaderField(Json::ValueType::arrayValue);

  jsonHeaderField.append(headerField.byteLocation);
  jsonHeaderField.append(ToString(headerField.fieldWidth));

  return jsonHeaderField;
}

Json::Value
SerializeSEGYFileInfo(SEGYFileInfo const& fileInfo, const int fileIndex)
{
  Json::Value
    jsonFileInfo;

  jsonFileInfo["persistentID"] = fmt::format("{:X}", fileInfo.m_persistentID);
  jsonFileInfo["headerEndianness"] = ToString(fileInfo.m_headerEndianness);
  jsonFileInfo["dataSampleFormatCode"] = (int)fileInfo.m_dataSampleFormatCode;
  jsonFileInfo["sampleCount"] = fileInfo.m_sampleCount;
  jsonFileInfo["startTime"] = fileInfo.m_startTimeMilliseconds;
  jsonFileInfo["sampleInterval"] = fileInfo.m_sampleIntervalMilliseconds;
  jsonFileInfo["traceCount"] = fileInfo.m_traceCounts[fileIndex];
  jsonFileInfo["primaryKey"] = SerializeSEGYHeaderField(fileInfo.m_primaryKey);
  jsonFileInfo["secondaryKey"] = SerializeSEGYHeaderField(fileInfo.m_secondaryKey);

  Json::Value
    jsonSegmentInfoArray(Json::ValueType::arrayValue);

  for (auto const& segmentInfo : fileInfo.m_segmentInfoLists[fileIndex])
  {
    jsonSegmentInfoArray.append(SerializeSEGYSegmentInfo(segmentInfo));
  }

  jsonFileInfo["segmentInfo"] = jsonSegmentInfoArray;

  return jsonFileInfo;
}

std::map<std::string, SEGY::HeaderField>
g_traceHeaderFields =
{
 { "TraceSequenceNumber"     , SEGY::TraceHeader::TraceSequenceNumberHeaderField },
 { "TraceSequenceNumberWithinFile", SEGY::TraceHeader::TraceSequenceNumberWithinFileHeaderField },
 { "EnergySourcePointNumber"   , SEGY::TraceHeader::EnergySourcePointNumberHeaderField },
 { "EnsembleNumber"        , SEGY::TraceHeader::EnsembleNumberHeaderField },
 { "TraceNumberWithinEnsemble"  , SEGY::TraceHeader::TraceNumberWithinEnsembleHeaderField },
 { "TraceIdentificationCode"   , SEGY::TraceHeader::TraceIdentificationCodeHeaderField },
 { "CoordinateScale"       , SEGY::TraceHeader::CoordinateScaleHeaderField },
 { "SourceXCoordinate"      , SEGY::TraceHeader::SourceXCoordinateHeaderField },
 { "SourceYCoordinate"      , SEGY::TraceHeader::SourceYCoordinateHeaderField },
 { "GroupXCoordinate"       , SEGY::TraceHeader::GroupXCoordinateHeaderField },
 { "GroupYCoordinate"       , SEGY::TraceHeader::GroupYCoordinateHeaderField },
 { "CoordinateUnits"       , SEGY::TraceHeader::CoordinateUnitsHeaderField },
 { "StartTime"          , SEGY::TraceHeader::StartTimeHeaderField },
 { "NumSamples"          , SEGY::TraceHeader::NumSamplesHeaderField },
 { "SampleInterval"        , SEGY::TraceHeader::SampleIntervalHeaderField },
 { "EnsembleXCoordinate"     , SEGY::TraceHeader::EnsembleXCoordinateHeaderField },
 { "EnsembleYCoordinate"     , SEGY::TraceHeader::EnsembleYCoordinateHeaderField },
 { "InlineNumber"         , SEGY::TraceHeader::InlineNumberHeaderField },
 { "CrosslineNumber"       , SEGY::TraceHeader::CrosslineNumberHeaderField },
 { "Receiver"              , SEGY::TraceHeader::ReceiverHeaderField },
 { "Offset"                , SEGY::TraceHeader::OffsetHeaderField }
};

std::map<std::string, std::string>
g_aliases =
{
 { "Inline",       "InlineNumber" },
 { "InLine",       "InlineNumber" },
 { "InLineNumber",    "InlineNumber" },
 { "Crossline",      "CrosslineNumber" },
 { "CrossLine",      "CrosslineNumber" },
 { "CrossLineNumber",   "CrosslineNumber" },
 { "Easting",       "EnsembleXCoordinate" },
 { "Northing",      "EnsembleYCoordinate" },
 { "CDPXCoordinate",   "EnsembleXCoordinate" },
 { "CDPYCoordinate",   "EnsembleYCoordinate" },
 { "CDP-X",        "EnsembleXCoordinate" },
 { "CDP-Y",        "EnsembleYCoordinate" },
 { "Source-X",      "SourceXCoordinate" },
 { "Source-Y",      "SourceYCoordinate" },
 { "Group-X",       "GroupXCoordinate" },
 { "Group-Y",       "GroupYCoordinate" },
 { "ReceiverXCoordinate", "GroupXCoordinate" },
 { "ReceiverYCoordinate", "GroupYCoordinate" },
 { "Receiver-X",     "GroupXCoordinate" },
 { "Receiver-Y",     "GroupYCoordinate" },
 { "Scalar",       "CoordinateScale" }
};

void
ResolveAlias(std::string& fieldName)
{
  if (g_aliases.find(fieldName) != g_aliases.end())
  {
    fieldName = g_aliases[fieldName];
  }
}

SEGY::Endianness
EndiannessFromJson(Json::Value const& jsonEndianness)
{
  std::string
    endiannessString = jsonEndianness.asString();

  if (endiannessString == "BigEndian")
  {
    return SEGY::Endianness::BigEndian;
  }
  else if (endiannessString == "LittleEndian")
  {
    return SEGY::Endianness::LittleEndian;
  }

  throw Json::Exception("Illegal endianness");
}

SEGY::FieldWidth
FieldWidthFromJson(Json::Value const& jsonFieldWidth)
{
  std::string
    fieldWidthString = jsonFieldWidth.asString();

  if (fieldWidthString == "TwoByte")
  {
    return SEGY::FieldWidth::TwoByte;
  }
  else if (fieldWidthString == "FourByte")
  {
    return SEGY::FieldWidth::FourByte;
  }

  throw Json::Exception("Illegal field width");
}

SEGY::HeaderField
HeaderFieldFromJson(Json::Value const& jsonHeaderField)
{
  int
    bytePosition = jsonHeaderField[0].asInt();

  SEGY::FieldWidth
    fieldWidth = FieldWidthFromJson(jsonHeaderField[1]);

  if (bytePosition < 1 || bytePosition > SEGY::TraceHeaderSize - ((fieldWidth == SEGY::FieldWidth::TwoByte) ? 2 : 4))
  {
    throw Json::Exception(std::string("Illegal field definition: ") + jsonHeaderField.toStyledString());
  }

  return SEGY::HeaderField(bytePosition, fieldWidth);
}

bool
ParseHeaderFormatFile(OpenVDS::File const& file, std::map<std::string, SEGY::HeaderField>& traceHeaderFields, SEGY::Endianness& headerEndianness, OpenVDS::Error &error)
{
  int64_t fileSize = file.Size(error);

  if (error.code != 0)
  {
    return false;
  }

  if (fileSize > INT_MAX)
  {
    return false;
  }

  std::unique_ptr<char[]>
    buffer(new char[fileSize]);

  file.Read(buffer.get(), 0, (int32_t)fileSize, error);

  if (error.code != 0)
  {
    return false;
  }

  try
  {
    Json::CharReaderBuilder
      rbuilder;

    rbuilder["collectComments"] = false;

    std::string
      errs;

    std::unique_ptr<Json::CharReader>
      reader(rbuilder.newCharReader());

    Json::Value
      root;

    bool
      success = reader->parse(buffer.get(), buffer.get() + fileSize, &root, &errs);

    if (!success)
    {
      throw Json::Exception(errs);
    }

    for (std::string const& fieldName : root.getMemberNames())
    {
      std::string canonicalFieldName = fieldName;
      ResolveAlias(canonicalFieldName);

      if (fieldName == "Endianness")
      {
        headerEndianness = EndiannessFromJson(root[fieldName]);
      }
      else
      {
        traceHeaderFields[canonicalFieldName] = HeaderFieldFromJson(root[fieldName]);
      }
    }
  }
  catch (Json::Exception &e)
  {
    error.code = -1;
    error.string = e.what();
    return false;
  }

  return true;
}

SEGYBinInfo
binInfoFromJson(Json::Value const& jsonBinInfo)
{
  int inlineNumber = jsonBinInfo["inlineNumber"].asInt();
  int crosslineNumber = jsonBinInfo["crosslineNumber"].asInt();
  double ensembleXCoordinate = jsonBinInfo["ensembleXCoordinate"].asDouble();
  double ensembleYCoordinate = jsonBinInfo["ensembleYCoordinate"].asDouble();

  return SEGYBinInfo(inlineNumber, crosslineNumber, ensembleXCoordinate, ensembleYCoordinate);
}

SEGYSegmentInfo
segmentInfoFromJson(Json::Value const& jsonSegmentInfo)
{
  int primaryKey = jsonSegmentInfo["primaryKey"].asInt();
  int traceStart = jsonSegmentInfo["traceStart"].asInt();
  int traceStop = jsonSegmentInfo["traceStop"].asInt();
  SEGYBinInfo binInfoStart = binInfoFromJson(jsonSegmentInfo["binInfoStart"]);
  SEGYBinInfo binInfoStop = binInfoFromJson(jsonSegmentInfo["binInfoStop"]);

  return SEGYSegmentInfo(primaryKey, traceStart, traceStop, binInfoStart, binInfoStop);
}

std::vector<int>
getOrderedSegmentListIndices(SEGYFileInfo const& fileInfo, size_t& globalTotalSegments)
{
  // Generate a list of indices that will traverse m_segyFileInfo.m_segmentInfoLists in primary key order, which
  // may be different from the order that the files were given.

  std::vector<int>
    orderedListIndices;
  int
    longestList = 0;
  globalTotalSegments = 0;
  for (int i = 0; i < fileInfo.m_segmentInfoLists.size(); ++i)
  {
    orderedListIndices.push_back(i);
    globalTotalSegments += fileInfo.m_segmentInfoLists[i].size();
    if (fileInfo.m_segmentInfoLists[i].size() > fileInfo.m_segmentInfoLists[longestList].size())
    {
      longestList = i;
    }
  }
  const bool
    isAscending = fileInfo.m_segmentInfoLists[longestList].front().m_binInfoStart.m_inlineNumber <= fileInfo.m_segmentInfoLists[longestList].back().m_binInfoStart.m_inlineNumber;
  auto
    comparator = [&](int i1, int i2)
  {
    const auto
      & v1 = fileInfo.m_segmentInfoLists[i1],
      & v2 = fileInfo.m_segmentInfoLists[i2];
    return isAscending ? v1.front().m_binInfoStart.m_inlineNumber < v2.front().m_binInfoStart.m_inlineNumber : v2.front().m_binInfoStart.m_inlineNumber < v1.front().m_binInfoStart.m_inlineNumber;
  };
  std::sort(orderedListIndices.begin(), orderedListIndices.end(), comparator);

  return orderedListIndices;
}

SEGYSegmentInfo const&
findRepresentativeSegment(SEGYFileInfo const& fileInfo, int& primaryStep, int& bestListIndex)
{
  // Since we give more weight to segments near the center of the data we need a sorted index of segment lists so that we can
  // traverse the lists in data order, instead of the arbitrary order given by the filename ordering.
  size_t
    globalTotalSegments = 0;
  auto
    orderedListIndices = getOrderedSegmentListIndices(fileInfo, globalTotalSegments);

  primaryStep = 0;
  bestListIndex = 0;

  float bestScore = 0.0f;
  int bestIndex = 0;

  int segmentPrimaryStep = 0;

  size_t
    globalOffset = 0;

  for (const auto listIndex : orderedListIndices)
  {
    const auto&
      segmentInfoList = fileInfo.m_segmentInfoLists[listIndex];

    for (int i = 0; i < segmentInfoList.size(); i++)
    {
      int64_t
        numTraces = (segmentInfoList[i].m_traceStop - segmentInfoList[i].m_traceStart + 1);

      // index of this segment within the entirety of segments from all input files
      const auto
        globalIndex = globalOffset + i;

      float
        multiplier = 1.5f - abs(globalIndex - (float)globalTotalSegments / 2) / (float)globalTotalSegments; // give 50% more importance to a segment in the middle of the dataset

      float
        score = float(numTraces) * multiplier;

      if (score > bestScore)
      {
        bestScore = score;
        bestListIndex = listIndex;
        bestIndex = i;
      }

      // Updating the primary step with the step for the previous segment intentionally ignores the step of the last segment since it can be anomalous
      if (segmentPrimaryStep && (!primaryStep || std::abs(segmentPrimaryStep) < std::abs(primaryStep)))
      {
        primaryStep = segmentPrimaryStep;
      }

      if (i > 0)
      {
        segmentPrimaryStep = segmentInfoList[i].m_primaryKey - segmentInfoList[i - 1].m_primaryKey;
      }
    }

    globalOffset += segmentInfoList.size();
  }

  // If the primary step couldn't be determined, set it to the last step or 1
  primaryStep = primaryStep ? primaryStep : std::max(segmentPrimaryStep, 1);

  return fileInfo.m_segmentInfoLists[bestListIndex][bestIndex];
}

void
copySamples(const void* data, SEGY::BinaryHeader::DataSampleFormatCode dataSampleFormatCode, SEGY::Endianness endianness, float* target, int sampleStart, int sampleCount)
{
  if (dataSampleFormatCode == SEGY::BinaryHeader::DataSampleFormatCode::IBMFloat)
  {
    if(endianness == SEGY::Endianness::LittleEndian)
    {
      // Reverse endianness since Ibm2ieee expects big endian data
      const char * source = reinterpret_cast<const char*>((intptr_t)data + (size_t)sampleStart * 4);
      std::unique_ptr<char[]> temp(new char[sampleCount * 4]);
      for(int sample = 0; sample < sampleCount; sample++)
      {
        temp[sample * 4 + 0] = source[sample * 4 + 3];
        temp[sample * 4 + 1] = source[sample * 4 + 2];
        temp[sample * 4 + 2] = source[sample * 4 + 1];
        temp[sample * 4 + 3] = source[sample * 4 + 0];
      }
      SEGY::Ibm2ieee(target, temp.get(), sampleCount);
    }
    else
    {
      assert(endianness == SEGY::Endianness::BigEndian);
      SEGY::Ibm2ieee(target, reinterpret_cast<const uint32_t*>((intptr_t)data + (size_t)sampleStart * 4), sampleCount);
    }
  }
  else
  {
    assert(dataSampleFormatCode == SEGY::BinaryHeader::DataSampleFormatCode::IEEEFloat);
    if(endianness == SEGY::Endianness::LittleEndian)
    {
      SEGY::ConvertFromEndianness<SEGY::Endianness::LittleEndian>(target, reinterpret_cast<const char*>((intptr_t)data + (size_t)sampleStart * 4), sampleCount);
    }
    else
    {
      assert(endianness == SEGY::Endianness::BigEndian);
      SEGY::ConvertFromEndianness<SEGY::Endianness::BigEndian>(target, reinterpret_cast<const char*>((intptr_t)data + (size_t)sampleStart * 4), sampleCount);
    }
  }
}

bool
analyzeSegment(DataProvider &dataProvider, SEGYFileInfo const& fileInfo, SEGYSegmentInfo const& segmentInfo, float valueRangePercentile, OpenVDS::FloatRange& valueRange, int& fold, int& secondaryStep, const SEGY::SEGYType segyType, int& offsetStart, int& offsetEnd, int& offsetStep, OpenVDS::Error& error)
{
  assert(segmentInfo.m_traceStop >= segmentInfo.m_traceStart && "A valid segment info should always have a stop trace greater or equal to the start trace");

  bool success = true;

  valueRange = OpenVDS::FloatRange(0.0f, 1.0f);
  secondaryStep = 0;
  fold = 1;
  offsetStart = 0;
  offsetEnd = 0;
  offsetStep = 0;

  const int traceByteSize = fileInfo.TraceByteSize();

  int64_t traceBufferStart = 0;
  int traceBufferSize = 0;
  std::unique_ptr<char[]> buffer;

  // Create min/max heaps for determining value range
  int heapSizeMax = int(((100.0f - valueRangePercentile) / 100.0f) * (segmentInfo.m_traceStop - segmentInfo.m_traceStart + 1) * fileInfo.m_sampleCount / 2) + 1;

  std::vector<float> minHeap, maxHeap;

  minHeap.reserve(heapSizeMax);
  maxHeap.reserve(heapSizeMax);

  // Allocate sample buffer for converting samples to float
  std::unique_ptr<float[]> sampleBuffer(new float[fileInfo.m_sampleCount]);
  float* samples = sampleBuffer.get();

  // Determine fold and secondary step
  int gatherSecondaryKey = 0, gatherFold = 0, gatherSecondaryStep = 0;

  bool
    hasPreviousGatherOffset = false;
  int
    previousGatherOffset = 0;

  for (int64_t trace = segmentInfo.m_traceStart; trace <= segmentInfo.m_traceStop; trace++)
  {
    if(trace - traceBufferStart >= traceBufferSize)
    {
      traceBufferStart = trace;
      traceBufferSize = (segmentInfo.m_traceStop - trace + 1) < 1000 ? int(segmentInfo.m_traceStop - trace + 1) : 1000;

      buffer.reset(new char[traceByteSize * traceBufferSize]);
      int64_t offset = SEGY::TextualFileHeaderSize + SEGY::BinaryFileHeaderSize + traceByteSize * traceBufferStart;
      success = dataProvider.Read(buffer.get(), offset, traceByteSize * traceBufferSize, error);

      if (!success)
      {
        break;
      }
    }

    const void *header = buffer.get() + traceByteSize * (trace - traceBufferStart);
    const void *data   = buffer.get() + traceByteSize * (trace - traceBufferStart) + SEGY::TraceHeaderSize;

    int tracePrimaryKey = SEGY::ReadFieldFromHeader(header, fileInfo.m_primaryKey, fileInfo.m_headerEndianness);
    int traceSecondaryKey = fileInfo.IsUnbinned() ? static_cast<int>(trace - segmentInfo.m_traceStart + 1) : SEGY::ReadFieldFromHeader(header, fileInfo.m_secondaryKey, fileInfo.m_headerEndianness);

    if(tracePrimaryKey != segmentInfo.m_primaryKey)
    {
      fmt::print(stderr, "Warning: trace {} has a primary key that doesn't match with the segment. This trace will be ignored.\n", segmentInfo.m_traceStart + trace);
      continue;
    }

    if(gatherFold > 0 && traceSecondaryKey == gatherSecondaryKey)
    {
      gatherFold++;
      fold = std::max(fold, gatherFold);
    }
    else
    {
      // Updating the secondary step with the step for the previous gather intentionally ignores the step of the last gather since it can be anomalous
      if(gatherSecondaryStep && (!secondaryStep || std::abs(gatherSecondaryStep) < std::abs(secondaryStep)))
      {
        secondaryStep = gatherSecondaryStep;
      }

      if(gatherFold > 0)
      {
        gatherSecondaryStep = traceSecondaryKey - gatherSecondaryKey;
      }
      gatherSecondaryKey = traceSecondaryKey;
      gatherFold = 1;
    }

    if (fileInfo.HasGatherOffset())
    {
      auto
        thisOffset = SEGY::ReadFieldFromHeader(header, g_traceHeaderFields["Offset"], fileInfo.m_headerEndianness);
      if (hasPreviousGatherOffset)
      {
        offsetStart = std::min(offsetStart, thisOffset);
        offsetEnd = std::max(offsetEnd, thisOffset);
        if (thisOffset != previousGatherOffset)
        {
          offsetStep = std::min(offsetStep, std::abs(thisOffset - previousGatherOffset));
        }
      }
      else
      {
        offsetStart = thisOffset;
        offsetEnd = thisOffset;
        offsetStep = INT32_MAX;
        hasPreviousGatherOffset = true;
      }
      previousGatherOffset = thisOffset;
    }

    // Update value range
    if (fileInfo.m_dataSampleFormatCode == SEGY::BinaryHeader::DataSampleFormatCode::IBMFloat || fileInfo.m_dataSampleFormatCode == SEGY::BinaryHeader::DataSampleFormatCode::IEEEFloat)
    {
      copySamples(data, fileInfo.m_dataSampleFormatCode, fileInfo.m_headerEndianness, samples, 0, fileInfo.m_sampleCount);

      for (int sample = 0; sample < fileInfo.m_sampleCount; sample++)
      {
        if (int(minHeap.size()) < heapSizeMax)
        {
          minHeap.push_back(samples[sample]);
          std::push_heap(minHeap.begin(), minHeap.end(), std::less<float>());
        }
        else if (samples[sample] < minHeap[0])
        {
          std::pop_heap(minHeap.begin(), minHeap.end(), std::less<float>());
          minHeap.back() = samples[sample];
          std::push_heap(minHeap.begin(), minHeap.end(), std::less<float>());
        }

        if (int(maxHeap.size()) < heapSizeMax)
        {
          maxHeap.push_back(samples[sample]);
          std::push_heap(maxHeap.begin(), maxHeap.end(), std::greater<float>());
        }
        else if (samples[sample] > maxHeap[0])
        {
          std::pop_heap(maxHeap.begin(), maxHeap.end(), std::greater<float>());
          maxHeap.back() = samples[sample];
          std::push_heap(maxHeap.begin(), maxHeap.end(), std::greater<float>());
        }
      }
    }
  }

  if (fileInfo.HasGatherOffset() && !fileInfo.IsUnbinned())
  {
    // check that offset start/end/step is consistent
    if (offsetStart + (fold - 1) * offsetStep != offsetEnd)
    {
      const auto
        msgFormat = "The detected gather offset start/end/step of '{0}/{1}/{2}' is not consistent with the detected fold of '{3}'. This usually indicates using the wrong header format for the input dataset.\n.";
      fmt::print(stderr, msgFormat, offsetStart, offsetEnd, offsetStep, fold);
      return false;

    }
  }

  // If the secondary step couldn't be determined, set it to the last step or 1
  secondaryStep = secondaryStep ? secondaryStep : std::max(gatherSecondaryStep, 1);

  // Set value range
  if (!minHeap.empty())
  {
    assert(!maxHeap.empty());

    if (minHeap[0] != maxHeap[0])
    {
      valueRange = OpenVDS::FloatRange(minHeap[0], maxHeap[0]);
    }
    else
    {
      valueRange = OpenVDS::FloatRange(minHeap[0], minHeap[0] + 1.0f);
    }
  }

  return success;
}

bool
createSEGYMetadata(DataProvider &dataProvider, SEGYFileInfo const &fileInfo, OpenVDS::MetadataContainer& metadataContainer, SEGY::BinaryHeader::MeasurementSystem &measurementSystem, OpenVDS::Error& error)
{
  std::vector<uint8_t> textHeader(SEGY::TextualFileHeaderSize);
  std::vector<uint8_t> binaryHeader(SEGY::BinaryFileHeaderSize);

  // Read headers
  bool success = dataProvider.Read(textHeader.data(), 0, SEGY::TextualFileHeaderSize, error) &&
    dataProvider.Read(binaryHeader.data(), SEGY::TextualFileHeaderSize, SEGY::BinaryFileHeaderSize, error);

  if (!success) return false;

  // Create metadata
  metadataContainer.SetMetadataBLOB("SEGY", "TextHeader", textHeader);

  metadataContainer.SetMetadataBLOB("SEGY", "BinaryHeader", binaryHeader);

  metadataContainer.SetMetadataInt("SEGY", "Endianness", int(fileInfo.m_headerEndianness));
  metadataContainer.SetMetadataInt("SEGY", "DataSampleFormatCode", int(fileInfo.m_dataSampleFormatCode));

  measurementSystem = SEGY::BinaryHeader::MeasurementSystem(SEGY::ReadFieldFromHeader(binaryHeader.data(), SEGY::BinaryHeader::MeasurementSystemHeaderField, fileInfo.m_headerEndianness));
  return success;
}

void
createSurveyCoordinateSystemMetadata(SEGYFileInfo const& fileInfo, SEGY::BinaryHeader::MeasurementSystem measurementSystem, std::string const &crsWkt, OpenVDS::MetadataContainer& metadataContainer)
{
  if (fileInfo.m_segmentInfoLists.empty()) return;

  double inlineSpacing[2] = { 0, 0 };
  double crosslineSpacing[2] = { 0, 0 };

  size_t
    globalTotalSegments;
  auto
    orderedListIndices = getOrderedSegmentListIndices(fileInfo, globalTotalSegments);

  // Determine crossline spacing
  int countedCrosslineSpacings = 0;

  for (size_t listIndex : orderedListIndices)
  {
    const auto
      & segmentInfoList = fileInfo.m_segmentInfoLists[listIndex];

    for (auto const& segmentInfo : segmentInfoList)
    {
      int crosslineCount = segmentInfo.m_binInfoStop.m_crosslineNumber - segmentInfo.m_binInfoStart.m_crosslineNumber;

      if (crosslineCount == 0 || segmentInfo.m_binInfoStart.m_inlineNumber != segmentInfo.m_binInfoStop.m_inlineNumber) continue;

      double segmentCrosslineSpacing[3];

      segmentCrosslineSpacing[0] = (segmentInfo.m_binInfoStop.m_ensembleXCoordinate - segmentInfo.m_binInfoStart.m_ensembleXCoordinate) / crosslineCount;
      segmentCrosslineSpacing[1] = (segmentInfo.m_binInfoStop.m_ensembleYCoordinate - segmentInfo.m_binInfoStart.m_ensembleYCoordinate) / crosslineCount;

      crosslineSpacing[0] += segmentCrosslineSpacing[0];
      crosslineSpacing[1] += segmentCrosslineSpacing[1];

      countedCrosslineSpacings++;
    }
  }

  if (countedCrosslineSpacings > 0)
  {
    crosslineSpacing[0] /= countedCrosslineSpacings;
    crosslineSpacing[1] /= countedCrosslineSpacings;
  }
  else
  {
    crosslineSpacing[0] = 0;
    crosslineSpacing[1] = 1;
  }

  // Determine inline spacing
  SEGYSegmentInfo const& firstSegmentInfo = fileInfo.m_segmentInfoLists[orderedListIndices.front()].front();
  SEGYSegmentInfo const& lastSegmentInfo = fileInfo.m_segmentInfoLists[orderedListIndices.back()].back();

  if (firstSegmentInfo.m_binInfoStart.m_inlineNumber != lastSegmentInfo.m_binInfoStart.m_inlineNumber)
  {
    int inlineNunberDelta = lastSegmentInfo.m_binInfoStart.m_inlineNumber - firstSegmentInfo.m_binInfoStart.m_inlineNumber;
    int crosslineNunberDelta = lastSegmentInfo.m_binInfoStart.m_crosslineNumber - firstSegmentInfo.m_binInfoStart.m_crosslineNumber;

    double offset[2] = { crosslineSpacing[0] * crosslineNunberDelta,
               crosslineSpacing[1] * crosslineNunberDelta };

    inlineSpacing[0] = (lastSegmentInfo.m_binInfoStart.m_ensembleXCoordinate - firstSegmentInfo.m_binInfoStart.m_ensembleXCoordinate - offset[0]) / inlineNunberDelta;
    inlineSpacing[1] = (lastSegmentInfo.m_binInfoStart.m_ensembleYCoordinate - firstSegmentInfo.m_binInfoStart.m_ensembleYCoordinate - offset[1]) / inlineNunberDelta;
  }
  else
  {
    // make square voxels
    inlineSpacing[0] = crosslineSpacing[1];
    inlineSpacing[1] = -crosslineSpacing[0];
  }

  // Determine origin
  double origin[2];

  origin[0] = firstSegmentInfo.m_binInfoStart.m_ensembleXCoordinate;
  origin[1] = firstSegmentInfo.m_binInfoStart.m_ensembleYCoordinate;

  origin[0] -= inlineSpacing[0] * firstSegmentInfo.m_binInfoStart.m_inlineNumber;
  origin[1] -= inlineSpacing[1] * firstSegmentInfo.m_binInfoStart.m_inlineNumber;

  origin[0] -= crosslineSpacing[0] * firstSegmentInfo.m_binInfoStart.m_crosslineNumber;
  origin[1] -= crosslineSpacing[1] * firstSegmentInfo.m_binInfoStart.m_crosslineNumber;

  // Set coordinate system
  metadataContainer.SetMetadataDoubleVector2(LATTICE_CATEGORY, LATTICE_ORIGIN, OpenVDS::DoubleVector2(origin[0], origin[1]));
  metadataContainer.SetMetadataDoubleVector2(LATTICE_CATEGORY, LATTICE_INLINE_SPACING, OpenVDS::DoubleVector2(inlineSpacing[0], inlineSpacing[1]));
  metadataContainer.SetMetadataDoubleVector2(LATTICE_CATEGORY, LATTICE_CROSSLINE_SPACING, OpenVDS::DoubleVector2(crosslineSpacing[0], crosslineSpacing[1]));

  if(!crsWkt.empty())
  {
    metadataContainer.SetMetadataString(LATTICE_CATEGORY, CRS_WKT, crsWkt);
  }
  if(measurementSystem == SEGY::BinaryHeader::MeasurementSystem::Meters)
  {
    metadataContainer.SetMetadataString(LATTICE_CATEGORY, LATTICE_UNIT, KNOWNMETADATA_UNIT_METER);
  }
  else if(measurementSystem == SEGY::BinaryHeader::MeasurementSystem::Feet)
  {
    metadataContainer.SetMetadataString(LATTICE_CATEGORY, LATTICE_UNIT, KNOWNMETADATA_UNIT_FOOT);
  }
}

/////////////////////////////////////////////////////////////////////////////
bool
createImportInformationMetadata(const std::vector<DataProvider> &dataProviders, OpenVDS::MetadataContainer& metadataContainer, OpenVDS::Error &error)
{
  auto now = std::chrono::system_clock::now();
  std::time_t tt = std::chrono::system_clock::to_time_t(now);
  std::tm tm = *std::gmtime(&tt);
  auto duration = now.time_since_epoch();
  int millis = std::chrono::duration_cast<std::chrono::milliseconds>(duration).count() % 1000;

  std::string importTimeStamp = fmt::format("{:04d}-{:02d}-{:02d}T{:02d}:{:02d}:{:02d}.{:03d}Z", tm.tm_year + 1900, tm.tm_mon + 1, tm.tm_mday, tm.tm_hour, tm.tm_min, tm.tm_sec, millis);

  // Join all the input names into a single comma delimited string
  const char pathSeparators[] = { '/', '\\', ':' };
  std::stringstream
    allNames;
  size_t
    nameCount = 0;
  for (auto iter = dataProviders.begin(); iter != dataProviders.end(); ++iter, ++nameCount)
  {
    // Strip the path from the file/object name
    std::string
      currentName = iter->FileOrObjectName();
    for (auto pathSeparator : pathSeparators)
    {
      size_t pos = currentName.rfind(pathSeparator);
      if (pos != std::string::npos) currentName = currentName.substr(pos + 1);
    }

    allNames << currentName;
    if (nameCount < dataProviders.size() - 1)
    {
      allNames << ",";
    }
  }

  std::string inputFileName = allNames.str();

  // In lack of a better displayName we use the file name
  std::string displayName = inputFileName;

  // Use the timestamp from the first input
  std::string inputTimeStamp = dataProviders[0].LastWriteTime(error);
  if (error.code != 0)
  {
    return false;
  }

  // Sum the sizes of all the inputs
  int64_t
    inputFileSize = 0;
  for (const auto& provider : dataProviders)
  {
    inputFileSize += provider.Size(error);
    if (error.code != 0)
    {
      return false;
    }
  }

  // Set import information
  metadataContainer.SetMetadataString(KNOWNMETADATA_CATEGORY_IMPORTINFORMATION, KNOWNMETADATA_IMPORTINFORMATION_DISPLAYNAME, displayName);
  metadataContainer.SetMetadataString(KNOWNMETADATA_CATEGORY_IMPORTINFORMATION, KNOWNMETADATA_IMPORTINFORMATION_INPUTFILENAME, inputFileName);
  metadataContainer.SetMetadataDouble(KNOWNMETADATA_CATEGORY_IMPORTINFORMATION, KNOWNMETADATA_IMPORTINFORMATION_INPUTFILESIZE, (double)inputFileSize);
  metadataContainer.SetMetadataString(KNOWNMETADATA_CATEGORY_IMPORTINFORMATION, KNOWNMETADATA_IMPORTINFORMATION_INPUTTIMESTAMP, inputTimeStamp);
  metadataContainer.SetMetadataString(KNOWNMETADATA_CATEGORY_IMPORTINFORMATION, KNOWNMETADATA_IMPORTINFORMATION_IMPORTTIMESTAMP, importTimeStamp);
  return true;
}

bool
parseSEGYFileInfoFile(OpenVDS::File const& file, SEGYFileInfo& fileInfo)
{
  OpenVDS::Error error;

  int64_t fileSize = file.Size(error);

  if (error.code != 0)
  {
    return false;
  }

  if (fileSize > INT_MAX)
  {
    return false;
  }

  std::unique_ptr<char[]>
    buffer(new char[fileSize]);

  file.Read(buffer.get(), 0, (int32_t)fileSize, error);

  if (error.code != 0)
  {
    return false;
  }

  try
  {
    Json::CharReaderBuilder
      rbuilder;

    rbuilder["collectComments"] = false;

    std::string
      errs;

    std::unique_ptr<Json::CharReader>
      reader(rbuilder.newCharReader());

    Json::Value
      jsonFileInfo;

    bool
      success = reader->parse(buffer.get(), buffer.get() + fileSize, &jsonFileInfo, &errs);

    if (!success)
    {
      throw Json::Exception(errs);
    }

    fileInfo.m_persistentID = strtoull(jsonFileInfo["persistentID"].asCString(), nullptr, 16);
    fileInfo.m_headerEndianness = EndiannessFromJson(jsonFileInfo["headerEndianness"]);
    fileInfo.m_dataSampleFormatCode = SEGY::BinaryHeader::DataSampleFormatCode(jsonFileInfo["dataSampleFormatCode"].asInt());
    fileInfo.m_sampleCount = jsonFileInfo["sampleCount"].asInt();
    fileInfo.m_startTimeMilliseconds = jsonFileInfo["startTime"].asDouble();
    fileInfo.m_sampleIntervalMilliseconds = jsonFileInfo["sampleInterval"].asDouble();
    fileInfo.m_traceCounts.clear();
    fileInfo.m_traceCounts.push_back(jsonFileInfo["traceCount"].asInt64());
    fileInfo.m_primaryKey = HeaderFieldFromJson(jsonFileInfo["primaryKey"]);
    fileInfo.m_secondaryKey = HeaderFieldFromJson(jsonFileInfo["secondaryKey"]);

    fileInfo.m_segmentInfoLists.clear();
    fileInfo.m_segmentInfoLists.emplace_back();
    auto&
      segmentInfo = fileInfo.m_segmentInfoLists.back();
    for (Json::Value jsonSegmentInfo : jsonFileInfo["segmentInfo"])
    {
      segmentInfo.push_back(segmentInfoFromJson(jsonSegmentInfo));
    }
  }
  catch (Json::Exception &e)
  {
    std::cerr << "Failed to parse JSON SEG-Y file info file: " << e.what();
    return false;
  }

  return true;
}

std::vector<OpenVDS::VolumeDataAxisDescriptor>
createAxisDescriptors(SEGYFileInfo const& fileInfo, SEGY::SampleUnits sampleUnits, int fold, int inlineStep, int crosslineStep)