#include <iostream> // // Created by michael on 23/08/17. // //---------------------- // Includes //---------------------- #include "MantidNexusGeometry/NexusGeometryParser.h" #include <string> #include <vector> namespace Mantid { namespace NexusGeometry { using namespace H5; //Eigen typedefs typedef Eigen::Matrix<double, 3, Eigen::Dynamic> Pixels; //Anonymous namespace for H5 constants 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 DETECTOR_IDS = "detector_number"; 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"; //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"; } /// Constructor opens the nexus file NexusGeometryParser::NexusGeometryParser( H5std_string &fileName, iAbstractBuilder_sptr iAbsBuilder_sptr ) { // Disable automatic printing, so Load algorithm can deal with errors // appropriately Exception::dontPrint(); try { H5File file(fileName, H5F_ACC_RDONLY); this->nexusFile = file; this->rootGroup = this->nexusFile.openGroup("/"); } catch (FileIException e) { this->exitStatus = OPENING_FILE_ERROR; } catch (GroupIException e) { this->exitStatus = OPENING_ROOT_GROUP_ERROR; } //Initialize the instrumentAbstractBuilder this->iBuilder_sptr = iAbsBuilder_sptr; } /// OFF NEXUS GEOMETRY PARSER ParsingErrors NexusGeometryParser::ParseNexusGeometry() { // Determine if nexusFile was successfully opened switch (this->exitStatus){ case NO_ERROR: break; default: return this->exitStatus; } // Get path to all detector groups try{ std::vector<Group> detectorGroups = this->openDetectorGroups(); for(std::vector<Group>::iterator iter = detectorGroups.begin(); iter < detectorGroups.end(); ++iter) { //Get the pixel offsets Pixels pixelOffsets = this->getPixelOffsets(*iter); //Get the transformations Eigen::Transform<double, 3, 2> transforms = this->getTransformations(*iter); //Calculate pixel positions Pixels detectorPixels = transforms * pixelOffsets; //Get the pixel detIds std::vector<int> detectorIds = this->getDetectorIds(*iter); // Force call of shape this->parseNexusShape(*iter); for(size_t i = 0; i < detectorIds.size(); ++i) { int index = static_cast<int>(i); std::string name = std::to_string(index); Eigen::Vector3d detPos = detectorPixels.col(index); this->iBuilder_sptr->addDetector(name, detectorIds[index], detPos, this->shape); } } //Sort the detectors this->iBuilder_sptr->sortDetectors(); //Parse source and sample and add to instrument this->parseAndAddSource(); this->parseAndAddSample(); } catch (...) { this->exitStatus = UNKNOWN_ERROR; } return this->exitStatus; } /// Open subgroups of parent group std::vector<Group> NexusGeometryParser::openSubGroups(Group &parentGroup, 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(int i = 0; i < childGroup.getNumAttrs(); ++i){ //Test attribute at current index for NX_class Attribute attribute = childGroup.openAttribute(i); 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> NexusGeometryParser::openDetectorGroups () { std::vector<Group> rawDataGroupPaths = this->openSubGroups (this->rootGroup, NX_ENTRY); //Open all instrument groups within rawDataGroups std::vector<Group> instrumentGroupPaths; for (std::vector<Group>::iterator iter = rawDataGroupPaths.begin (); iter < rawDataGroupPaths.end (); iter++) { std::vector<Group> instrumentGroups = this->openSubGroups (*iter, NX_INSTRUMENT); instrumentGroupPaths.insert (instrumentGroupPaths.end (), instrumentGroups.begin (), instrumentGroups.end ()); } //Open all detector groups within instrumentGroups std::vector<Group> detectorGroupPaths; for(std::vector<Group>::iterator iter = instrumentGroupPaths.begin(); iter < instrumentGroupPaths.end (); iter++) { //Open sub detector groups std::vector<Group> detectorGroups = this->openSubGroups (*iter, NX_DETECTOR); //Append to detectorGroups vector detectorGroupPaths.insert (detectorGroupPaths.end (), detectorGroups.begin (), detectorGroups.end ()); } //Return the detector groups return detectorGroupPaths; } //Function to return the detector ids in the same order as the opffsets std::vector<int> NexusGeometryParser::getDetectorIds ( Group &detectorGroup) { H5std_string detectorName = detectorGroup.getObjName(); std::vector<int> detIds; for (unsigned int i = 0; i < detectorGroup.getNumObjs (); ++i) { H5std_string objName = detectorGroup.getObjnameByIdx (i); H5std_string objPath = detectorName + "/" + objName; if(objName == DETECTOR_IDS) { detIds = this->get1DDataset<int> (objPath); } } return detIds; } //Function to return the (x,y,z) offsets of pixels in the chosen detectorGroup Pixels NexusGeometryParser::getPixelOffsets (Group &detectorGroup) { H5std_string detectorName = detectorGroup.getObjName(); //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); H5std_string objPath = detectorName + "/" + objName; if(objName == X_PIXEL_OFFSET) { xValues = this->get1DDataset<double>(objPath); } if(objName == Y_PIXEL_OFFSET) { yValues = this->get1DDataset<double>(objPath); } if(objName == Z_PIXEL_OFFSET) { zValues = this->get1DDataset<double> (objPath); } } // 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; } //Function to read in a dataset into a vector template<typename valueType> std::vector<valueType> NexusGeometryParser::get1DDataset(H5std_string &dataset) { //Open data set DataSet data = this->nexusFile.openDataSet (dataset); //Get data type DataType dataType = data.getDataType (); //Get data size DataSpace dataSpace = data.getSpace (); //Create vector to hold data std::vector<valueType> values; values.resize (dataSpace.getSelectNpoints ()); //Read data into vector data.read (values.data (), dataType, dataSpace); //Return the data vector return values; } H5std_string NexusGeometryParser::get1DStringDataset(const H5std_string &dataset) { //Open data set DataSet data = this->nexusFile.openDataSet (dataset); //Get size DataSpace dataspace = data.getSpace (); //Initialize string H5std_string value; //Read into string data.read(value, data.getDataType (), dataspace); return value; } //Function to get the transformations from the nexus file, and create the Eigen transform object Eigen::Transform<double, 3, Eigen::Affine> NexusGeometryParser::getTransformations(Group &detectorGroup) { //Get path to depends_on file H5std_string detectorPath = detectorGroup.getObjName (); H5std_string depends_on = detectorPath + "/" + DEPENDS_ON; //Get absolute dependency path H5std_string dependency = this->get1DStringDataset(depends_on); //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 opposite order to direction of discovery while(dependency != NO_DEPENDENCY) { //Open the transformation data set DataSet transformation = this->nexusFile.openDataSet (dependency); //Get magnitude of current transformation double magnitude = this->get1DDataset<double>(dependency)[0]; //Containers for transformation data Eigen::Vector3d transformVector(0.0,0.0,0.0); H5std_string transformType; H5std_string transformUnits; for(int i = 0; i < 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; } 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; } } return transforms; } ///Choose what shape type to parse void NexusGeometryParser::parseNexusShape(Group &detectorGroup) { Group shapeGroup; try{ shapeGroup = detectorGroup.openGroup(PIXEL_SHAPE); } catch(...){ //Placeholder - no group pixel_shape } H5std_string shapeType; for (int i = 0; i < 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) { this->parseNexusCylinder(shapeGroup); } } //Parse cylinder nexus geometry void NexusGeometryParser::parseNexusCylinder(Group &shapeGroup){ H5std_string pointsToVertices = shapeGroup.getObjName() + "/cylinder"; std::vector<int> cPoints = this->get1DDataset<int>(pointsToVertices); H5std_string verticesData = shapeGroup.getObjName() + "/vertices"; // 1D reads row first, then columns std::vector<double> vPoints = this->get1DDataset<double>(verticesData); 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); } this->shape = this->sAbsCreator.createCylinder(vSorted); } //Parse source and add to instrument void NexusGeometryParser::parseAndAddSource() { H5std_string sourcePath = "raw_data_1/instrument/source"; Group sourceGroup = this->rootGroup.openGroup(sourcePath); auto sourceName = this->get1DStringDataset("/name"); auto defaultPos = Eigen::Vector3d(0.0,0.0,0.0); // TODO incomplete } //Parse sample and add to instrument void NexusGeometryParser::parseAndAddSample() { std::string sampleName = "sample"; H5std_string samplePath = "raw_data_1/sample"; Group sampleGroup = this->rootGroup.openGroup(samplePath); auto sampleTransforms = this->getTransformations(sampleGroup); auto samplePos = sampleTransforms*Eigen::Vector3d(0.0,0.0,0.0); this->iBuilder_sptr->addSample(sampleName, samplePos); } } }