NexusGeometryParser.cpp

#include "MantidNexusGeometry/NexusGeometryParser.h"
#include "MantidGeometry/Instrument.h"
#include "MantidGeometry/Objects/CSGObject.h"
#include "MantidGeometry/Objects/ShapeFactory.h"
#include "MantidGeometry/Rendering/GeometryHandler.h"
#include "MantidGeometry/Rendering/ShapeInfo.h"
#include "MantidKernel/make_unique.h"
#include "MantidNexusGeometry/InstrumentBuilder.h"
#include "MantidNexusGeometry/NexusShapeFactory.h"
#include "MantidKernel/EigenConversionHelpers.h"
#include <Eigen/Core>
#include <Eigen/Geometry>
#include <H5Cpp.h>
#include <boost/algorithm/string.hpp>
#include <type_traits>
#include <numeric>
#include <tuple>

namespace Mantid {
namespace NexusGeometry {

using namespace H5;

// Eigen typedefs
typedef Eigen::Matrix<double, 3, Eigen::Dynamic> Pixels;

// Anonymous namespace
namespace {
const H5G_obj_t GROUP_TYPE = static_cast<H5G_obj_t>(0);
const H5std_string NX_CLASS = "NX_class";
const H5std_string NX_ENTRY = "NXentry";
const H5std_string NX_INSTRUMENT = "NXinstrument";
const H5std_string NX_DETECTOR = "NXdetector";
const H5std_string NX_MONITOR = "NXmonitor";
const H5std_string DETECTOR_IDS = "detector_number";
const H5std_string DETECTOR_ID = "detector_id";
const H5std_string X_PIXEL_OFFSET = "x_pixel_offset";
const H5std_string Y_PIXEL_OFFSET = "y_pixel_offset";
const H5std_string Z_PIXEL_OFFSET = "z_pixel_offset";
const H5std_string DEPENDS_ON = "depends_on";
const H5std_string NO_DEPENDENCY = ".";
const H5std_string PIXEL_SHAPE = "pixel_shape";
const H5std_string DETECTOR_SHAPE = "detector_shape";
const H5std_string SHAPE = "shape";
// Transformation types
const H5std_string TRANSFORMATION_TYPE = "transformation_type";
const H5std_string TRANSLATION = "translation";
const H5std_string ROTATION = "rotation";
const H5std_string VECTOR = "vector";
const H5std_string UNITS = "units";
// Radians and degrees
const H5std_string DEGREES = "degrees";
const static double PI = 3.1415926535;
const static double DEGREES_IN_SEMICIRCLE = 180;
// Nexus shape types
const H5std_string NX_CYLINDER = "NXcylindrical_geometry";
const H5std_string NX_OFF = "NXoff_geometry";
const H5std_string BANK_NAME = "local_name";

using FaceV = std::vector<Eigen::Vector3d>;

struct Face {
  Eigen::Vector3d v1;
  Eigen::Vector3d v2;
  Eigen::Vector3d v3;
  Eigen::Vector3d v4;
};

template <typename T, typename R>
std::vector<R> convertVector(const std::vector<T> &toConvert) {
  std::vector<R> target(toConvert.size());
  for (size_t i = 0; i < target.size(); ++i) {
    target[i] = R(toConvert[i]);
  }
  return target;
}

template <typename ExpectedT> void validateStorageType(const DataSet &data) {

  const auto typeClass = data.getTypeClass();
  const size_t sizeOfType = data.getDataType().getSize();
  // Early check to prevent reinterpretation of underlying data.
  if (std::is_floating_point<ExpectedT>::value) {
    if (H5T_FLOAT != typeClass) {
      throw std::runtime_error("Storage type mismatch. Expecting to extract a "
                               "floating point number");
    }
    if (sizeOfType != sizeof(ExpectedT)) {
      throw std::runtime_error(
          "Storage type mismatch for floats. This operation "
          "is dangerous. Nexus stored has byte size:" +
          std::to_string(sizeOfType));
    }
  } else if (std::is_integral<ExpectedT>::value) {
    if (H5T_INTEGER != typeClass) {
      throw std::runtime_error(
          "Storage type mismatch. Expecting to extract a integer");
    }
    if (sizeOfType > sizeof(ExpectedT)) {
      // endianness not checked
      throw std::runtime_error(
          "Storage type mismatch for integer. Result "
          "would result in truncation. Nexus stored has byte size:" +
          std::to_string(sizeOfType));
    }
  }
}

template <typename ValueType>
std::vector<ValueType> extractVector(const DataSet &data) {
  validateStorageType<ValueType>(data);
  DataSpace dataSpace = data.getSpace();
  std::vector<ValueType> values;
  values.resize(dataSpace.getSelectNpoints());
  // Read data into vector
  data.read(values.data(), data.getDataType(), dataSpace);

  // Return the data vector
  return values;
}

// Function to read in a dataset into a vector
template <typename ValueType>
std::vector<ValueType> get1DDataset(const H5std_string &dataset,
                                    const H5::Group &group) {
  // Open data set
  DataSet data = group.openDataSet(dataset);
  return extractVector<ValueType>(data);
}

// Function to read in a dataset into a vector
template <typename ValueType>
std::vector<ValueType> get1DDataset(const H5File &file,
                                    const H5std_string &dataset) {
  // Open data set
  DataSet data = file.openDataSet(dataset);
  return extractVector<ValueType>(data);
}

std::string get1DStringDataset(const std::string &dataset, const Group &group) {
  // Open data set
  DataSet data = group.openDataSet(dataset);
  auto dataType = data.getDataType();
  auto nCharacters = dataType.getSize();
  std::vector<char> value(nCharacters);
  data.read(value.data(), dataType, data.getSpace());
  return std::string(value.begin(), value.end());
}

std::string instrumentName(const Group &root) {
  H5std_string instrumentPath = "raw_data_1/instrument";
  const Group instrumentGroup = root.openGroup(instrumentPath);

  return get1DStringDataset("name", instrumentGroup);
}

/// Open subgroups of parent group
std::vector<Group> openSubGroups(const Group &parentGroup,
                                 const H5std_string &CLASS_TYPE) {
  std::vector<Group> subGroups;

  // Iterate over children, and determine if a group
  for (hsize_t i = 0; i < parentGroup.getNumObjs(); ++i) {
    if (parentGroup.getObjTypeByIdx(i) == GROUP_TYPE) {
      H5std_string childPath = parentGroup.getObjnameByIdx(i);
      // Open the sub group
      Group childGroup = parentGroup.openGroup(childPath);
      // Iterate through attributes to find NX_class
      for (uint32_t attribute_index = 0;
           attribute_index < static_cast<uint32_t>(childGroup.getNumAttrs());
           ++attribute_index) {
        // Test attribute at current index for NX_class
        Attribute attribute = childGroup.openAttribute(attribute_index);
        if (attribute.getName() == NX_CLASS) {
          // Get attribute data type
          DataType dataType = attribute.getDataType();
          // Get the NX_class type
          H5std_string classType;
          attribute.read(dataType, classType);
          // If group of correct type, append to subGroup vector
          if (classType == CLASS_TYPE) {
            subGroups.push_back(childGroup);
          }
        }
      }
    }
  }
  return subGroups;
}

// Open all detector groups into a vector
std::vector<Group> openDetectorGroups(const Group &root) {
  std::vector<Group> rawDataGroupPaths = openSubGroups(root, NX_ENTRY);

  // Open all instrument groups within rawDataGroups
  std::vector<Group> instrumentGroupPaths;
  for (auto const &rawDataGroupPath : rawDataGroupPaths) {
    std::vector<Group> instrumentGroups =
        openSubGroups(rawDataGroupPath, NX_INSTRUMENT);
    instrumentGroupPaths.insert(instrumentGroupPaths.end(),
                                instrumentGroups.begin(),
                                instrumentGroups.end());
  }
  // Open all detector groups within instrumentGroups
  std::vector<Group> detectorGroupPaths;
  for (auto const &instrumentGroupPath : instrumentGroupPaths) {
    // Open sub detector groups
    std::vector<Group> detectorGroups =
        openSubGroups(instrumentGroupPath, NX_DETECTOR);
    // Append to detectorGroups vector
    detectorGroupPaths.insert(detectorGroupPaths.end(), detectorGroups.begin(),
                              detectorGroups.end());
  }
  // Return the detector groups
  return detectorGroupPaths;
}

// Function to return the (x,y,z) offsets of pixels in the chosen detectorGroup
Pixels getPixelOffsets(const Group &detectorGroup) {

  // Initialise matrix
  Pixels offsetData;
  std::vector<double> xValues, yValues, zValues;
  for (unsigned int i = 0; i < detectorGroup.getNumObjs(); i++) {
    H5std_string objName = detectorGroup.getObjnameByIdx(i);
    if (objName == X_PIXEL_OFFSET) {
      xValues = get1DDataset<double>(objName, detectorGroup);
    }
    if (objName == Y_PIXEL_OFFSET) {
      yValues = get1DDataset<double>(objName, detectorGroup);
    }
    if (objName == Z_PIXEL_OFFSET) {
      zValues = get1DDataset<double>(objName, detectorGroup);
    }
  }

  // Determine size of dataset
  int rowLength = 0;
  bool xEmpty = xValues.empty();
  bool yEmpty = yValues.empty();
  bool zEmpty = zValues.empty();

  if (!xEmpty)
    rowLength = static_cast<int>(xValues.size());
  else if (!yEmpty)
    rowLength = static_cast<int>(yValues.size());
  // Need at least 2 dimensions to define points
  else
    return offsetData;

  // Default x,y,z to zero if no data provided
  offsetData.resize(3, rowLength);
  offsetData.setZero(3, rowLength);

  if (!xEmpty) {
    for (int i = 0; i < rowLength; i++)
      offsetData(0, i) = xValues[i];
  }
  if (!yEmpty) {
    for (int i = 0; i < rowLength; i++)
      offsetData(1, i) = yValues[i];
  }
  if (!zEmpty) {
    for (int i = 0; i < rowLength; i++)
      offsetData(2, i) = zValues[i];
  }
  // Return the coordinate matrix
  return offsetData;
}

/**
 * Creates a Homogemous transfomation for nexus groups
 *
 * Walks the chain of transformations described in the file where W1 is first
 *transformation and Wn is last and assembles them as
 *
 * W = Wn x ... W2 x W1
 *
 * Each W describes a Homogenous Transformation
 *
 * R | T
 * -   -
 * 0 | 1
 *
 *
 * @param file
 * @param detectorGroup
 * @return
 */
Eigen::Transform<double, 3, Eigen::Affine>
getTransformations(const H5File &file, const Group &detectorGroup) {
  H5std_string dependency;
  // Get absolute dependency path
  auto status =
      H5Gget_objinfo(detectorGroup.getId(), DEPENDS_ON.c_str(), false, nullptr);
  if (status == 0) {
    dependency = get1DStringDataset(DEPENDS_ON, detectorGroup);
  } else {
    return Eigen::Transform<double, 3, Eigen::Affine>::Identity();
  }

  // Initialise transformation holder as zero-degree rotation
  Eigen::Transform<double, 3, Eigen::Affine> transforms;
  Eigen::Vector3d axis(1.0, 0.0, 0.0);
  transforms = Eigen::AngleAxisd(0.0, axis);

  // Breaks when no more dependencies (dependency = ".")
  // Transformations must be applied in the order of direction of discovery
  // (they are _passive_ transformations)
  while (dependency != NO_DEPENDENCY) {
    // Open the transformation data set
    DataSet transformation = file.openDataSet(dependency);

    // Get magnitude of current transformation
    double magnitude = get1DDataset<double>(file, dependency)[0];
    // Containers for transformation data
    Eigen::Vector3d transformVector(0.0, 0.0, 0.0);
    H5std_string transformType;
    H5std_string transformUnits;
    for (uint32_t i = 0;
         i < static_cast<uint32_t>(transformation.getNumAttrs()); i++) {
      // Open attribute at current index
      Attribute attribute = transformation.openAttribute(i);
      H5std_string attributeName = attribute.getName();
      // Get next dependency
      if (attributeName == DEPENDS_ON) {
        DataType dataType = attribute.getDataType();
        attribute.read(dataType, dependency);
      }
      // Get transform type
      else if (attributeName == TRANSFORMATION_TYPE) {
        DataType dataType = attribute.getDataType();
        attribute.read(dataType, transformType);
      }
      // Get unit vector for transformation
      else if (attributeName == VECTOR) {
        std::vector<double> unitVector;
        DataType dataType = attribute.getDataType();

        // Get data size
        DataSpace dataSpace = attribute.getSpace();

        // Resize vector to hold data
        unitVector.resize(dataSpace.getSelectNpoints());

        // Read the data into Eigen vector
        attribute.read(dataType, unitVector.data());
        transformVector(0) = unitVector[0];
        transformVector(1) = unitVector[1];
        transformVector(2) = unitVector[2];
      } else if (attributeName == UNITS) {
        DataType dataType = attribute.getDataType();
        attribute.read(dataType, transformUnits);
      }
    }
    // Transform_type = translation
    if (transformType == TRANSLATION) {
      // Translation = magnitude*unitVector
      transformVector *= magnitude;
      Eigen::Translation3d translation(transformVector);
      transforms = translation * transforms;
    } else if (transformType == ROTATION) {
      double angle = magnitude;
      if (transformUnits == DEGREES) {
        // Convert angle from degrees to radians
        angle *= PI / DEGREES_IN_SEMICIRCLE;
      }
      Eigen::AngleAxisd rotation(angle, transformVector);
      transforms = rotation * transforms;
    }
  }
  return transforms;
}

// Function to return the detector ids in the same order as the offsets
std::vector<int> getDetectorIds(const Group &detectorGroup) {

  std::vector<int> detIds;

  for (unsigned int i = 0; i < detectorGroup.getNumObjs(); ++i) {
    H5std_string objName = detectorGroup.getObjnameByIdx(i);
    if (objName == DETECTOR_IDS) {
      const auto data = detectorGroup.openDataSet(objName);
      if (data.getDataType().getSize() == 8) {
        // Note the narrowing here!
        detIds = convertVector<int64_t, int32_t>(extractVector<int64_t>(data));
      } else {
        detIds = extractVector<int32_t>(data);
      }
    }
  }
  return detIds;
}

// Parse cylinder nexus geometry
boost::shared_ptr<const Geometry::IObject>
parseNexusCylinder(const Group &shapeGroup) {
  H5std_string pointsToVertices = "cylinders";
  std::vector<int> cPoints = get1DDataset<int>(pointsToVertices, shapeGroup);

  H5std_string verticesData = "vertices";
  // 1D reads row first, then columns
  std::vector<double> vPoints = get1DDataset<double>(verticesData, shapeGroup);
  Eigen::Map<Eigen::Matrix<double, 3, 3>> vertices(vPoints.data());
  // Read points into matrix, sorted by cPoints ordering
  Eigen::Matrix<double, 3, 3> vSorted;
  for (int i = 0; i < 3; ++i) {
    vSorted.col(cPoints[i]) = vertices.col(i);
  }
  return NexusShapeFactory::createCylinder(vSorted);
}

// Parse OFF (mesh) nexus geometry
boost::shared_ptr<const Geometry::IObject>
parseNexusMesh(const Group &shapeGroup) {
  const std::vector<uint16_t> faceIndices = convertVector<int32_t, uint16_t>(
      get1DDataset<int32_t>("faces", shapeGroup));
  const std::vector<uint16_t> windingOrder = convertVector<int32_t, uint16_t>(
      get1DDataset<int32_t>("winding_order", shapeGroup));
  const auto vertices = get1DDataset<float>("vertices", shapeGroup);
  return NexusShapeFactory::createFromOFFMesh(faceIndices, windingOrder,
                                              vertices);
}

void extractFacesAndIDs(const std::vector<uint16_t> &detFaces,
                        const std::vector<uint16_t> &windingOrder,
                        const std::vector<float> &vertices,
                        const std::unordered_map<int, uint32_t> &detIdToIndex,
                        const size_t vertsPerFace,
                        std::vector<std::vector<Eigen::Vector3d>> &detFaceVerts,
                        std::vector<std::vector<uint16_t>> &detFaceIndices,
                        std::vector<std::vector<uint16_t>> &detWindingOrder,
                        std::vector<int32_t> &detIds) {
  const size_t vertStride = 3;
  size_t detFaceIndex = 1;
  std::fill(detFaceIndices.begin(), detFaceIndices.end(),
            std::vector<uint16_t>(1, 0));
  for (size_t i = 0; i < windingOrder.size(); i += vertsPerFace) {
    auto detFaceId = detFaces[detFaceIndex];
    // Id -> Index
    auto detIndex = detIdToIndex.at(detFaceId);
    auto &detVerts = detFaceVerts[detIndex];
    auto &detIndices = detFaceIndices[detIndex];
    auto &detWinding = detWindingOrder[detIndex];
    detVerts.reserve(vertsPerFace);
    detWinding.reserve(vertsPerFace);
    for (size_t v = 0; v < vertsPerFace; ++v) {
      auto vi = windingOrder[i + v] * vertStride;
      detVerts.emplace_back(vertices[vi], vertices[vi + 1], vertices[vi + 2]);
      detWinding.push_back(static_cast<uint16_t>(detWinding.size()));
    }
    // Index -> Id
    detIds[detIndex] = detFaceId;
    detIndices.push_back(static_cast<uint16_t>(detVerts.size()));
    // Detector faces is 2N detectors
    detFaceIndex += 2;
  }
}

void parseNexusMeshAndAddDetectors(
    const std::vector<uint16_t> &detFaces,
    const std::vector<uint16_t> &faceIndices,
    const std::vector<uint16_t> &windingOrder,
    const std::vector<float> &vertices, const size_t numDets,
    const std::unordered_map<int, uint32_t> &detIdToIndex,
    const std::string &name, InstrumentBuilder &builder) {
  auto vertsPerFace = windingOrder.size() / faceIndices.size();
  std::vector<std::vector<Eigen::Vector3d>> detFaceVerts(numDets);
  std::vector<std::vector<uint16_t>> detFaceIndices(numDets);
  std::vector<std::vector<uint16_t>> detWindingOrder(numDets);
  std::vector<int> detIds(numDets);

  extractFacesAndIDs(detFaces, windingOrder, vertices, detIdToIndex,
                     vertsPerFace, detFaceVerts, detFaceIndices,
                     detWindingOrder, detIds);

  for (size_t i = 0; i < numDets; ++i) {
    auto &detVerts = detFaceVerts[i];
    const auto &detIndices = detFaceIndices[i];
    const auto &detWinding = detWindingOrder[i];
    // Calculate polygon centre
    auto centre = std::accumulate(detVerts.begin() + 1, detVerts.end(),
                                  detVerts.front()) /
                  detVerts.size();

    // translate shape to origin for shape coordinates.
    std::for_each(detVerts.begin(), detVerts.end(),
                  [&centre](Eigen::Vector3d &val) { val -= centre; });

    auto shape =
        NexusShapeFactory::createFromOFFMesh(detIndices, detWinding, detVerts);
    builder.addDetectorToLastBank(name + "_" + std::to_string(i), detIds[i],
                                  centre, std::move(shape));
  }
}

void parseAndAddBank(const Group &shapeGroup, InstrumentBuilder &builder,
                     const std::vector<int> &detectorIds,
                     const std::string &bankName) {
  // Load mapping between detector IDs and faces, winding order of vertices for
  // faces, and face corner vertices.
  const std::vector<uint16_t> detFaces = convertVector<int32_t, uint16_t>(
      get1DDataset<int32_t>("detector_faces", shapeGroup));
  const std::vector<uint16_t> faceIndices = convertVector<int32_t, uint16_t>(
      get1DDataset<int32_t>("faces", shapeGroup));
  const std::vector<uint16_t> windingOrder = convertVector<int32_t, uint16_t>(
      get1DDataset<int32_t>("winding_order", shapeGroup));
  const auto vertices = get1DDataset<float>("vertices", shapeGroup);

  // Build a map of detector IDs to the index of occurrence in the
  // "detector_number" dataset
  std::unordered_map<int, uint32_t> detIdToIndex;
  for (uint32_t i = 0; i < detectorIds.size(); ++i) {
    detIdToIndex.emplace(detectorIds[i], i);
  }

  parseNexusMeshAndAddDetectors(detFaces, faceIndices, windingOrder, vertices,
                                detectorIds.size(), detIdToIndex, bankName,
                                builder);
}

/// Choose what shape type to parse
boost::shared_ptr<const Geometry::IObject>
parseNexusShape(const Group &detectorGroup, bool &searchTubes) {
  bool isGroup = false;
  Group shapeGroup;
  try {
    shapeGroup = detectorGroup.openGroup(PIXEL_SHAPE);
    isGroup = true;
  } catch (H5::Exception &) {
    // TODO. Current assumption. Can we have pixels without specifying a
    // shape?
    try {
      shapeGroup = detectorGroup.openGroup(SHAPE);
    } catch (H5::Exception &) {
      return boost::shared_ptr<const Geometry::IObject>(nullptr);
    }
  }

  H5std_string shapeType;
  for (uint32_t i = 0; i < static_cast<uint32_t>(shapeGroup.getNumAttrs());
       ++i) {
    Attribute attribute = shapeGroup.openAttribute(i);
    H5std_string attributeName = attribute.getName();
    if (attributeName == NX_CLASS) {
      attribute.read(attribute.getDataType(), shapeType);
    }
  }
  // Give shape group to correct shape parser
  if (shapeType == NX_CYLINDER) {
    searchTubes = isGroup;
    return parseNexusCylinder(shapeGroup);
  } else if (shapeType == NX_OFF) {
    return parseNexusMesh(shapeGroup);
  } else {
    throw std::runtime_error(
        "Shape type not recognised by NexusGeometryParser");
  }
}

// Parse source and add to instrument
void parseAndAddSource(const H5File &file, const Group &root,
                       InstrumentBuilder &builder) {
  H5std_string sourcePath = "raw_data_1/instrument/source";
  Group sourceGroup = root.openGroup(sourcePath);
  auto sourceName = get1DStringDataset("name", sourceGroup);
  auto sourceTransformations = getTransformations(file, sourceGroup);
  auto defaultPos = Eigen::Vector3d(0.0, 0.0, 0.0);
  builder.addSource(sourceName, sourceTransformations * defaultPos);
}

// Parse sample and add to instrument
void parseAndAddSample(const H5File &file, const Group &root,
                       InstrumentBuilder &builder) {
  std::string sampleName = "sample";
  H5std_string samplePath = "raw_data_1/sample";
  Group sampleGroup = root.openGroup(samplePath);
  auto sampleTransforms = getTransformations(file, sampleGroup);
  auto samplePos = sampleTransforms * Eigen::Vector3d(0.0, 0.0, 0.0);
  builder.addSample(sampleName, samplePos);
}

void parseMonitors(const H5::Group &root, InstrumentBuilder &builder) {
  std::vector<Group> rawDataGroupPaths = openSubGroups(root, NX_ENTRY);

  // Open all instrument groups within rawDataGroups
  for (auto const &rawDataGroupPath : rawDataGroupPaths) {
    std::vector<Group> instrumentGroups =
        openSubGroups(rawDataGroupPath, NX_INSTRUMENT);
    for (auto &inst : instrumentGroups) {
      std::vector<Group> monitorGroups = openSubGroups(inst, NX_MONITOR);
      for (auto &monitor : monitorGroups) {
        auto detectorId = get1DDataset<int64_t>(DETECTOR_ID, monitor)[0];
        bool proxy = false;
        boost::shared_ptr<const Geometry::IObject> monitorShape =
            parseNexusShape(monitor, proxy);
        builder.addMonitor(std::to_string(detectorId),
                           static_cast<int32_t>(detectorId),
                           Eigen::Vector3d{0, 0, 0}, monitorShape);
      }
    }
  }
}

class Tube {
private:
  Eigen::Vector3d axis;
  double height = 0;
  double radius;
  std::vector<Eigen::Vector3d> m_positions;
  std::vector<int> m_detIDs;
  Eigen::Vector3d p1;
  Eigen::Vector3d p2;

public:
  const Eigen::Vector3d &location() const { return m_positions[0]; }
  const std::vector<Eigen::Vector3d> &positions() const { return m_positions; }
  const std::vector<int> &detIDs() const { return m_detIDs; }

  boost::shared_ptr<const Mantid::Geometry::IObject> shape() {
    Eigen::Matrix<double, 3, 3> points;
    // Centre shape about (0, 0, 0)
    auto p1 = m_positions.front();
    auto p2 = m_positions.back();
    auto centre = (p1 + p2) / 2;

    points.col(0) = p1 - centre;
    points.col(2) = p2 - centre;

    // Find point on outer circle
    auto norm = axis.norm();
    auto factor = radius / norm;
    auto pointOnCircle = (axis * factor) + (p1 - centre);
    points.col(1) = pointOnCircle;

    return NexusShapeFactory::createCylinder(points);
  }

  Tube(Eigen::Vector3d axis, double radius, Eigen::Vector3d pos, int detID) {
    this->axis = axis;
    this->radius = radius;
    m_positions.push_back(pos);
    m_detIDs.push_back(detID);

    p1 = axis + pos;
    p2 = pos;
  }

  bool addIfCoLinear(const Eigen::Vector3d &pos, int detID) {
    if (checkCoLinear(pos)) {
      m_positions.push_back(pos);
      m_detIDs.push_back(detID);

      auto dist = m_positions[0] - pos;
      if (dist.norm() > height)
        height = dist.norm();
      return true;
    }

    return false;
  }

private:
  bool checkCoLinear(const Eigen::Vector3d &pos) const {
    auto temp = pos + axis; // unecessary
    auto numVec = ((p2 - p1).cross(p1 - temp));
    auto denomVec = (p2 - p1);
    auto num = numVec.norm();

    if (num == 0)
      return true;

    auto denom = denomVec.norm();
    auto d = num / denom;

    return d == 0.0;
  }
};

std::vector<Tube> findTubes(const Mantid::Geometry::IObject &shape,
                            const Pixels &positions,
                            const std::vector<int> detIDs) {
  auto axis =
      Kernel::toVector3d(shape.getGeometryHandler()->shapeInfo().points()[1]);
  auto radius = shape.getGeometryHandler()->shapeInfo().radius();
  std::vector<Tube> tubes;

  tubes.emplace_back(Tube(axis, radius, positions.col(0), detIDs[0]));
  bool newEntry;

  for (size_t i = 1; i < detIDs.size(); ++i) {
    newEntry = true;
    for (auto t = tubes.rbegin(); t != tubes.rend(); ++t) {
      auto &tube = (*t);
      if (tube.addIfCoLinear(positions.col(i), detIDs[i])) {
        newEntry = false;
        break;
      }
    }

    if (newEntry)
      tubes.push_back(Tube(axis, radius, positions.col(i), detIDs[i]));
  }

  return tubes;
}

std::unique_ptr<const Mantid::Geometry::Instrument>
extractInstrument(const H5File &file, const Group &root) {
  InstrumentBuilder builder(instrumentName(root));
  // Get path to all detector groups
  const std::vector<Group> detectorGroups = openDetectorGroups(root);
  for (auto &detectorGroup : detectorGroups) {
    // Transform in homogenous coordinates. Offsets will be rotated then bank
    // translation applied.
    Eigen::Transform<double, 3, 2> transforms =
        getTransformations(file, detectorGroup);
    // Absolute bank position
    Eigen::Vector3d bankPos = transforms * Eigen::Vector3d{0, 0, 0};
    // Absolute bank rotation
    auto bankRotation = Eigen::Quaterniond(transforms.rotation());
    auto bankName = get1DStringDataset(BANK_NAME, detectorGroup);
    builder.addBank(bankName, bankPos, bankRotation);
    // Get the pixel detIds
    auto detectorIds = getDetectorIds(detectorGroup);

    try {
      auto shapeGroup = detectorGroup.openGroup(DETECTOR_SHAPE);
      parseAndAddBank(shapeGroup, builder, detectorIds, bankName);
      continue;
    } catch (H5::Exception &) { // No detector_shape group
    }

    // Get the pixel offsets
    Pixels pixelOffsets = getPixelOffsets(detectorGroup);
    // Calculate pixel relative positions
    Pixels detectorPixels = Eigen::Affine3d::Identity() * pixelOffsets;
    bool searchTubes = false;
    // Extract shape
    auto shape = parseNexusShape(detectorGroup, searchTubes);

    std::vector<Tube> tubes;
    if (searchTubes) {
      tubes = findTubes(*shape, pixelOffsets, detectorIds);
      for (size_t i = 0; i < tubes.size(); ++i) {
        auto name = bankName + "_tube_" + std::to_string(i);
        auto &tube = tubes[i];
        builder.addObjComponentAssembly(name, tube.location(), tube.shape(),
                                        shape, tube.positions(), tube.detIDs());
      }
    }
    else {
      for (size_t i = 0; i < detectorIds.size(); ++i) {
        auto index = static_cast<int>(i);
        std::string name = bankName + "_" + std::to_string(index);

        Eigen::Vector3d relativePos = detectorPixels.col(index);
        builder.addDetectorToLastBank(name, detectorIds[index], relativePos,
          shape);
      }
    }
  }
  // Sort the detectors
  // Parse source and sample and add to instrument
  parseAndAddSample(file, root, builder);
  parseAndAddSource(file, root, builder);
  parseMonitors(root, builder);
  return builder.createInstrument();
}
} // namespace

std::unique_ptr<const Geometry::Instrument>
NexusGeometryParser::createInstrument(const std::string &fileName) {

  const H5File file(fileName, H5F_ACC_RDONLY);
  auto rootGroup = file.openGroup("/");
  return extractInstrument(file, rootGroup);
}
} // namespace NexusGeometry
} // namespace Mantid