Skip to content

Commit

Permalink
first pass generic geom handling #37
Browse files Browse the repository at this point in the history
  • Loading branch information
@LamarMoore
LamarMoore committed May 2, 2018
1 parent 15e3dd1 commit 7f7bd11
Show file tree
Hide file tree
Showing 8 changed files with 226 additions and 142 deletions.
11 changes: 2 additions & 9 deletions Framework/Geometry/src/Instrument/StructuredDetector.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -398,15 +398,8 @@ Detector *StructuredDetector::addDetector(CompAssembly *parent,
auto yrb = m_yvalues[(y * w) + x + 1];

// calculate midpoint of trapeziod
auto a = std::fabs(xrf - xlf);
auto b = std::fabs(xrb - xlb);
auto h = std::fabs(ylb - ylf);
auto cx = ((a + b) / 4);
auto cy = h / 2;

// store detector position before translating to origin
auto xpos = xlb + cx;
auto ypos = ylb + cy;
auto xpos = (xlb + xlf + xrf + xrb) / 4;
auto ypos = (ylb + ylf + yrf + yrb) / 4;

// Translate detector shape to origin
xlf -= xpos;
Expand Down
7 changes: 5 additions & 2 deletions Framework/Geometry/src/Objects/MeshObject.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -58,8 +58,8 @@ const Kernel::Material MeshObject::material() const { return m_material; }
*/
bool MeshObject::hasValidShape() const {
// May enclose volume if there are at
// at least 4 triangles and 4 vertices
return (numberOfTriangles() >= 4 && numberOfVertices() >= 4);
// at least 2 triangles and 4 vertices
return (numberOfTriangles() >= 2 && numberOfVertices() >= 4);
}

/**
Expand Down Expand Up @@ -446,6 +446,9 @@ const BoundingBox &MeshObject::getBoundingBox() const {
minZ = std::min(minZ, vz);
maxZ = std::max(maxZ, vz);
}
if (minZ == maxZ)
maxZ += 0.001;

// Cache bounding box, so we do not need to repeat calculation
m_boundingBox = BoundingBox(maxX, maxY, maxZ, minX, minY, minZ);
}
Expand Down
6 changes: 6 additions & 0 deletions Framework/NexusGeometry/inc/MantidNexusGeometry/NexusShapeFactory.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -4,6 +4,7 @@
#include "MantidNexusGeometry/DllConfig.h"

#include "Eigen/Core"
#include "Eigen/Geometry"
#include <memory>
#include <map>
#include <vector>
Expand Down Expand Up @@ -57,6 +58,11 @@ DLLExport std::unique_ptr<const Geometry::IObject>
createFromOFFMesh(const std::vector<uint16_t> &faceIndices,
const std::vector<uint16_t> &windingOrder,
const std::vector<float> &nexusVertices);

std::unique_ptr<const Geometry::IObject>
createFromOFFMesh(const std::vector<uint16_t> &faceIndices,
const std::vector<uint16_t> &windingOrder,
const std::vector<Eigen::Vector3d> &nexusVertices);
}
}
}
Expand Down
249 changes: 158 additions & 91 deletions Framework/NexusGeometry/src/NexusGeometryParser.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -407,93 +407,181 @@ parseNexusCylinder(const Group &shapeGroup) {

// Parse OFF (mesh) nexus geometry
boost::shared_ptr<const Geometry::IObject>
parseNexusMesh(const Group &shapeGroup) {
parseNexusMesh(const Group &detectorGroup, const Group &shapeGroup) {
const std::vector<uint16_t> detFaces = convertVector<int32_t, uint16_t>(
get1DDataset<int32_t>("detector_faces", shapeGroup));
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);

Pixels pixelOffsets = getPixelOffsets(detectorGroup);
return NexusShapeFactory::createFromOFFMesh(faceIndices, windingOrder,
vertices);
}

std::pair<boost::shared_ptr<const Geometry::IObject>, Eigen::Vector3d> parseHex(const FaceV &face) {
auto xlb = face[0].x();
auto xlf = face[1].x();
auto xrf = face[2].x();
auto xrb = face[3].x();
auto ylb = face[0].y();
auto ylf = face[1].y();
auto yrf = face[2].y();
auto yrb = face[3].y();

// calculate midpoint of trapeziod
auto a = std::fabs(xrf - xlf);
auto b = std::fabs(xrb - xlb);
auto h = std::fabs(ylb - ylf);
auto cx = ((a + b) / 4);
auto cy = h / 2;

// store detector position before translating to origin
auto xpos = xlb + cx;
auto ypos = ylb + cy;

// Translate detector shape to origin
xlf -= xpos;
xrf -= xpos;
xrb -= xpos;
xlb -= xpos;
ylf -= ypos;
yrf -= ypos;
yrb -= ypos;
ylb -= ypos;

return{ Mantid::Geometry::ShapeFactory::createHexahedralShape(
xlb, xlf, xrf, xrb, ylb, ylf, yrf, yrb), Eigen::Vector3d(xpos, ypos, 0) };
}

std::vector<std::pair<boost::shared_ptr<const Geometry::IObject>, Eigen::Vector3d>>
parseNexusMeshForBank(const Group &shapeGroup, const size_t numDets, int32_t minDetId) {
//std::pair<boost::shared_ptr<const Geometry::IObject>, Eigen::Vector3d>
//parseHex(const FaceV &face) {
// auto xlb = face[0].x();
// auto xlf = face[1].x();
// auto xrf = face[2].x();
// auto xrb = face[3].x();
// auto ylb = face[0].y();
// auto ylf = face[1].y();
// auto yrf = face[2].y();
// auto yrb = face[3].y();
//
// // calculate midpoint of trapeziod
// auto a = std::fabs(xrf - xlf);
// auto b = std::fabs(xrb - xlb);
// auto h = std::fabs(ylb - ylf);
// auto cx = ((a + b) / 4);
// auto cy = h / 2;
//
// // store detector position before translating to origin
// auto xpos = xlb + cx;
// auto ypos = ylb + cy;
//
// // Translate detector shape to origin
// xlf -= xpos;
// xrf -= xpos;
// xrb -= xpos;
// xlb -= xpos;
// ylf -= ypos;
// yrf -= ypos;
// yrb -= ypos;
// ylb -= ypos;
//
// return {Mantid::Geometry::ShapeFactory::createHexahedralShape(
// xlb, xlf, xrf, xrb, ylb, ylf, yrf, yrb),
// Eigen::Vector3d(xpos, ypos, 0)};
//}

//std::vector<
// std::pair<boost::shared_ptr<const Geometry::IObject>, Eigen::Vector3d>>
//parseNexusMeshForBank(const Group &shapeGroup, const size_t numDets,
// int32_t minDetId) {
// 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);
//
// auto vertsPerFace = windingOrder.size() / faceIndices.size();
// std::vector<std::vector<FaceV>> detFaceMap;
// size_t vertStride = 3;
// size_t detFaceIndex = 1;
// detFaceMap.resize(numDets);
// for (size_t i = 0; i < windingOrder.size(); i += vertsPerFace) {
// FaceV fv(vertsPerFace);
// for (size_t v = 0; v < vertsPerFace; ++v) {
// auto vi = windingOrder[i + v] * vertStride;
// fv[v] = Eigen::Vector3d(vertices[vi], vertices[vi + 1], vertices[vi + 2]);
// }
// detFaceMap[detFaces[detFaceIndex] - minDetId].emplace_back(fv);
// detFaceIndex += 2;
// }
//
// std::vector<
// std::pair<boost::shared_ptr<const Geometry::IObject>, Eigen::Vector3d>>
// shapes;
// shapes.reserve(numDets);
// for (auto &faces : detFaceMap) {
// if (faces.size() == 1 && vertsPerFace == 4) { // create hexahedron
// auto &f = faces[0];
// shapes.emplace_back(parseHex(f));
// }
// }
//
// return shapes;
//}

std::vector<boost::shared_ptr<const Geometry::IObject>>
parseNexusMeshForBank(const Group &shapeGroup, const size_t numDets,
int32_t minDetId, const Pixels &pixelOffsets,
const Eigen::Vector3d &bankPosition,
Eigen::Quaterniond &rotation) {
const std::vector<uint16_t> detFaces = convertVector<int32_t, uint16_t>(
get1DDataset<int32_t>("detector_faces", shapeGroup));
get1DDataset<int32_t>("detector_faces", shapeGroup));
const std::vector<uint16_t> faceIndices = convertVector<int32_t, uint16_t>(
get1DDataset<int32_t>("faces", shapeGroup));
get1DDataset<int32_t>("faces", shapeGroup));
const std::vector<uint16_t> windingOrder = convertVector<int32_t, uint16_t>(
get1DDataset<int32_t>("winding_order", shapeGroup));
get1DDataset<int32_t>("winding_order", shapeGroup));
const auto vertices = get1DDataset<float>("vertices", shapeGroup);

auto vertsPerFace = windingOrder.size() / faceIndices.size();
std::vector<std::vector<FaceV>> detFaceMap;
std::vector<std::vector<Eigen::Vector3d>> detFaceVerts(numDets);
std::vector<std::vector<uint16_t>> detFaceIndices(numDets);
std::vector<std::vector<uint16_t>> detWindingOrder(numDets);
size_t vertStride = 3;
size_t detFaceIndex = 1;
detFaceMap.resize(numDets);
for (size_t i = 0; i < windingOrder.size(); i+=vertsPerFace) {
FaceV fv(vertsPerFace);
std::fill(detFaceIndices.begin(), detFaceIndices.end(),
std::vector<uint16_t>(1, 0));
for (size_t i = 0; i < windingOrder.size(); i += vertsPerFace) {
auto &detVerts = detFaceVerts[detFaces[detFaceIndex] - minDetId];
auto &detIndices = detFaceIndices[detFaces[detFaceIndex] - minDetId];
auto &detWinding = detWindingOrder[detFaces[detFaceIndex] - minDetId];
for (size_t v = 0; v < vertsPerFace; ++v) {
auto vi = windingOrder[i + v] * vertStride;
fv[v] = Eigen::Vector3d(vertices[vi], vertices[vi + 1], vertices[vi + 2]);
detVerts.emplace_back(vertices[vi], vertices[vi + 1] ,vertices[vi + 2]);
detWinding.push_back(static_cast<uint16_t>(detWinding.size()));
}
detFaceMap[detFaces[detFaceIndex]-minDetId].emplace_back(fv);
detIndices.push_back(static_cast<uint16_t>(detVerts.size()));
detFaceIndex += 2;
}

std::vector<
std::pair<boost::shared_ptr<const Geometry::IObject>, Eigen::Vector3d>>
shapes;
std::vector<boost::shared_ptr<const Geometry::IObject>> shapes;
shapes.reserve(numDets);
for (auto &faces : detFaceMap) {
if (faces.size() == 1 && vertsPerFace == 4) { // create hexahedron
auto &f = faces[0];
shapes.emplace_back(parseHex(f));
}
for (size_t i = 0; i < numDets; ++i) {
auto &detVerts = detFaceVerts[i];
const auto &detIndices = detFaceIndices[i];
const auto &detWinding = detWindingOrder[i];

std::for_each(detVerts.begin(), detVerts.end(),
[&pixelOffsets, i](Eigen::Vector3d &val) {
val -= pixelOffsets.col(i);
});

auto shape = NexusShapeFactory::createFromOFFMesh(
detIndices, detWinding, detVerts);
shapes.emplace_back(std::move(shape));
}

return shapes;
}

void parseAndAddBank(const H5File &file, const Group &detectorGroup,
InstrumentBuilder &builder, Eigen::Vector3d sourcePosition) {
// Get the pixel offsets
Pixels pixelOffsets = getPixelOffsets(detectorGroup);
// 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
Eigen::Quaterniond bankRotation = Eigen::Quaterniond(transforms.rotation());
builder.addBank(get1DStringDataset(BANK_NAME, detectorGroup), bankPos,
bankRotation);
// Get the pixel detIds
auto detectorIds = getDetectorIds(detectorGroup);
auto minId = std::min(detectorIds.cbegin(), detectorIds.cend());
// Extract shape
auto res = parseNexusMeshForBank(detectorGroup.openGroup(DETECTOR_SHAPE),
detectorIds.size(), *minId, pixelOffsets,
bankPos + sourcePosition, bankRotation);
for (size_t i = 0; i < detectorIds.size(); ++i) {
auto index = static_cast<int>(i);
std::string name = std::to_string(index);
auto &shape = res[i];
builder.addDetectorToLastBank(name, detectorIds[index], pixelOffsets.col(i),
shape);
}
}

/// Choose what shape type to parse
boost::shared_ptr<const Geometry::IObject>
parseNexusShape(const Group &detectorGroup) {
Expand Down Expand Up @@ -523,43 +611,13 @@ parseNexusShape(const Group &detectorGroup) {
if (shapeType == NX_CYLINDER) {
return parseNexusCylinder(shapeGroup);
} else if (shapeType == NX_OFF) {
return parseNexusMesh(shapeGroup);
return parseNexusMesh(detectorGroup, shapeGroup);
} else {
throw std::runtime_error(
"Shape type not recognised by NexusGeometryParser");
}
}

void parseAndAddBank(const H5File &file, const Group &detectorGroup, InstrumentBuilder &builder) {
// Get the pixel offsets
Pixels pixelOffsets = getPixelOffsets(detectorGroup);
// 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
Eigen::Quaterniond bankRotation = Eigen::Quaterniond(transforms.rotation());
builder.addBank(get1DStringDataset(BANK_NAME, detectorGroup), bankPos,
bankRotation);
// Get the pixel detIds
auto detectorIds = getDetectorIds(detectorGroup);
auto minId = std::min(detectorIds.cbegin(), detectorIds.cend());
// Extract shape
auto res = parseNexusMeshForBank(detectorGroup.openGroup(DETECTOR_SHAPE),
detectorIds.size(), *minId);

for (size_t i = 0; i < detectorIds.size(); ++i) {
auto index = static_cast<int>(i);
std::string name = std::to_string(index);
auto &shape = res[i].first;
auto &relativePos = res[i].second;
builder.addDetectorToLastBank(name, detectorIds[index], relativePos, shape);
}
}


// Parse source and add to instrument
void parseAndAddSource(const H5File &file, const Group &root,
InstrumentBuilder &builder) {
Expand All @@ -570,6 +628,7 @@ void parseAndAddSource(const H5File &file, const Group &root,
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) {
Expand All @@ -581,6 +640,13 @@ void parseAndAddSample(const H5File &file, const Group &root,
builder.addSample(sampleName, samplePos);
}

Eigen::Vector3d getSourcePosition(const H5File &file, const Group &root) {
H5std_string sourcePath = "raw_data_1/instrument/source";
Group sourceGroup = root.openGroup(sourcePath);
auto sourceTransformations = getTransformations(file, sourceGroup);
return sourceTransformations * Eigen::Vector3d(0.0, 0.0, 0.0);
}

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

Expand Down Expand Up @@ -610,9 +676,10 @@ extractInstrument(const H5File &file, const Group &root) {
for (auto &detectorGroup : detectorGroups) {
try {
detectorGroup.openGroup(DETECTOR_SHAPE);
parseAndAddBank(file, detectorGroup, builder);
parseAndAddBank(file, detectorGroup, builder,
getSourcePosition(file, root));
continue;
} catch (...) {//No detector_shape group
} catch (...) { // No detector_shape group
}
// Get the pixel offsets
Pixels pixelOffsets = getPixelOffsets(detectorGroup);
Expand All @@ -636,7 +703,7 @@ extractInstrument(const H5File &file, const Group &root) {
for (size_t i = 0; i < detectorIds.size(); ++i) {
auto index = static_cast<int>(i);
std::string name = std::to_string(index);

Eigen::Vector3d relativePos = detectorPixels.col(index);
builder.addDetectorToLastBank(name, detectorIds[index], relativePos,
shape);
Expand Down
Loading

0 comments on commit 7f7bd11

Please sign in to comment.