diff --git a/.github/workflows/build_test.yml b/.github/workflows/build_test.yml new file mode 100644 index 0000000..9ae81d7 --- /dev/null +++ b/.github/workflows/build_test.yml @@ -0,0 +1,25 @@ +on: [push, pull_request] + +jobs: + build_test: + name: Build/Test + strategy: + matrix: + os: [ubuntu, macOS] + defaults: + run: + shell: bash -l {0} + runs-on: ${{ matrix.os }}-latest + steps: + - uses: actions/checkout@v4 + - uses: conda-incubator/setup-miniconda@v3 + with: + channels: conda-forge + - name: Install dependencies + run: | + conda env update --file requirements.txt --name test + conda list + - name: Build + run: mkdir build && cd build && cmake .. && make + - name: Run tests + run: cd build && ctest -L dx2tests --output-on-failure diff --git a/CMakeLists.txt b/CMakeLists.txt index e96b0ab..bcf5570 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -18,6 +18,7 @@ option(BUILD_SHARED_LIBS "Build using shared libraries" ON) # ####################################################################### # External Dependencies find_package(HDF5 REQUIRED) +find_package(gemmi CONFIG REQUIRED) # ####################################################################### # Automatic Dependencies @@ -107,4 +108,4 @@ install( ) export(EXPORT DX2Targets FILE "${CMAKE_CURRENT_BINARY_DIR}/DX2Targets.cmake" -) \ No newline at end of file +) diff --git a/README.md b/README.md index 9e49c19..1e734c2 100644 --- a/README.md +++ b/README.md @@ -2,9 +2,10 @@ ## Quick Start -- Ensure you have an environment with CMake and compilers available +- Ensure you have an environment with CMake and compilers available (see requirements.txt) - Dependencies: - HDF5 (Required) + - gemmi (Required) - [Eigen](https://eigen.tuxfamily.org/) (Optional, downloaded if missing) - [nlohmann/json](https://github.com/nlohmann/json) (Optional, downloaded if missing) - [{fmt}](https://github.com/fmtlib/fmt) (Optional, downloaded if missing) @@ -15,11 +16,11 @@ mkdir build cd build cmake .. ``` -You can then build with `make`, and run the tests with `make test`. +You can then build with `make`, and run the test with `./tests/test_make_models`. ## Writing And Running Tests -See https://google.github.io/googletest/primer.html. The `tests/` subdirectory +See https://google.github.io/googletest/primer.html. The `tests` subdirectory has an example test. You can run test executables directly, run `ctest`, or run `make test` (or diff --git a/dx2/CMakeLists.txt b/dx2/CMakeLists.txt index 64f5ee6..b5c3d73 100644 --- a/dx2/CMakeLists.txt +++ b/dx2/CMakeLists.txt @@ -8,6 +8,7 @@ target_link_libraries( Eigen3::Eigen nlohmann_json::nlohmann_json hdf5::hdf5 + gemmi::gemmi_cpp ) target_include_directories( @@ -15,3 +16,4 @@ target_include_directories( $ $ ) +target_link_libraries(dx2 PUBLIC $) diff --git a/include/dx2/beam.h b/include/dx2/beam.h new file mode 100644 index 0000000..f5e6d39 --- /dev/null +++ b/include/dx2/beam.h @@ -0,0 +1,288 @@ +#ifndef DX2_MODEL_BEAM_H +#define DX2_MODEL_BEAM_H +#include +#include + +using Eigen::Vector3d; +using json = nlohmann::json; + +/* +This file defines a set of Beam classes. + +BeamBase defines most of the class attributes except the wavelength. + +The two subclasses are MonochromaticBeam and PolychromaticBeam, +which define a single wavelength or wavelength range respectively. + +MonochromaticBeam is subclassed further into MonoXrayBeam and +MonoElectronBeam, these simply set the correct probe name +when serializing/deserializing to/from json. +*/ + +class BeamBase { +// A base class for beam objects +public: + BeamBase()=default; + BeamBase(Vector3d direction, double divergence, + double sigma_divergence, Vector3d polarization_normal, + double polarization_fraction, double flux, + double transmission, double sample_to_source_distance); + +protected: + void init_from_json(json beam_data); + void add_to_json(json beam_data) const; + Vector3d sample_to_source_direction_{0.0,0.0,1.0}; //called direction_ in dxtbx + double divergence_{0.0}; // "beam divergence - be more specific with name?" + double sigma_divergence_{0.0}; // standard deviation of the beam divergence + Vector3d polarization_normal_{0.0,1.0,0.0}; + double polarization_fraction_{0.999}; + double flux_{0.0}; + double transmission_{1.0}; + double sample_to_source_distance_{0.0}; // FIXME is this really needed? +}; + +class MonochromaticBeam: public BeamBase { +// A monochromatic beam (i.e. single wavelength value) +public: + MonochromaticBeam()=default; + MonochromaticBeam(double wavelength); + MonochromaticBeam(Vector3d s0); + MonochromaticBeam(double wavelength, Vector3d direction, double divergence, + double sigma_divergence, Vector3d polarization_normal, + double polarization_fraction, double flux, + double transmission, double sample_to_source_distance); + MonochromaticBeam(json beam_data); + json to_json(std::string probe) const; + double get_wavelength() const; + void set_wavelength(double wavelength); + Vector3d get_s0() const; + void set_s0(Vector3d s0); + +protected: + double wavelength_{0.0}; +}; + + +class MonoXrayBeam : public MonochromaticBeam { +// Same as the parent class, except explicitly set a probe type when calling to_json +using MonochromaticBeam::MonochromaticBeam; +public: + json to_json() const; +}; + +class MonoElectronBeam : public MonochromaticBeam { +// Same as the parent class, except explicitly set a probe type when calling to_json +using MonochromaticBeam::MonochromaticBeam; +public: + json to_json() const; +}; + + +class PolychromaticBeam : public BeamBase { +// A polychromatic beam (i.e. a wavelength range) +public: + PolychromaticBeam()=default; + PolychromaticBeam(std::array wavelength_range); + PolychromaticBeam(std::array wavelength_range, Vector3d direction); + PolychromaticBeam(std::array wavelength_range, Vector3d direction, double divergence, + double sigma_divergence, Vector3d polarization_normal, + double polarization_fraction, double flux, + double transmission, double sample_to_source_distance); + PolychromaticBeam(json beam_data); + json to_json(std::string probe) const; + std::array get_wavelength_range() const; + void set_wavelength_range(std::array wavelength_range); + +protected: + std::array wavelength_range_{{0.0, 0.0}}; +}; + + +// BeamBase definitions + +// full constructor +BeamBase::BeamBase( + Vector3d direction, double divergence, + double sigma_divergence, Vector3d polarization_normal, + double polarization_fraction, double flux, + double transmission, double sample_to_source_distance) + : sample_to_source_direction_{direction}, + divergence_{divergence}, + sigma_divergence_{sigma_divergence}, + polarization_normal_{polarization_normal}, + polarization_fraction_{polarization_fraction}, + flux_{flux}, transmission_{transmission}, + sample_to_source_distance_{sample_to_source_distance} {} + +void BeamBase::init_from_json(json beam_data){ + // Load values from a json object. + // Allow these to not be present so we can load a simple json dict without all these items. + if (beam_data.find("direction") != beam_data.end()){ + Vector3d direction{{beam_data["direction"][0], beam_data["direction"][1], beam_data["direction"][2]}}; + sample_to_source_direction_ = direction; + } + if (beam_data.find("divergence") != beam_data.end()){ + divergence_ = beam_data["divergence"]; + } + if (beam_data.find("sigma_divergence") != beam_data.end()){ + sigma_divergence_ = beam_data["sigma_divergence"]; + } + if (beam_data.find("polarization_normal") != beam_data.end()){ + Vector3d pn{ + {beam_data["polarization_normal"][0], + beam_data["polarization_normal"][1], + beam_data["polarization_normal"][2]}}; + polarization_normal_ = pn; + } + if (beam_data.find("polarization_fraction") != beam_data.end()){ + polarization_fraction_ = beam_data["polarization_fraction"]; + } + if (beam_data.find("flux") != beam_data.end()){ + flux_ = beam_data["flux"]; + } + if (beam_data.find("transmission") != beam_data.end()){ + transmission_ = beam_data["transmission"]; + } + if (beam_data.find("sample_to_source_distance") != beam_data.end()){ + sample_to_source_distance_ = beam_data["sample_to_source_distance"]; + } +} + +void BeamBase::add_to_json(json beam_data) const { + // Add the members to the json object to prepare for serialization. + beam_data["direction"] = sample_to_source_direction_; + beam_data["divergence"] = divergence_; + beam_data["sigma_divergence"] = sigma_divergence_; + beam_data["polarization_normal"] = polarization_normal_; + beam_data["polarization_fraction"] = polarization_fraction_; + beam_data["flux"] = flux_; + beam_data["transmission"] = transmission_; + beam_data["sample_to_source_distance"] = sample_to_source_distance_; +} + + +// MonochromaticBeam definitions + +// full constructor +MonochromaticBeam::MonochromaticBeam( + double wavelength, Vector3d direction, double divergence, + double sigma_divergence, Vector3d polarization_normal, + double polarization_fraction, double flux, + double transmission, double sample_to_source_distance) + : BeamBase{direction, divergence, sigma_divergence, + polarization_normal, polarization_fraction, + flux, transmission, sample_to_source_distance}, wavelength_{wavelength} {} + +MonochromaticBeam::MonochromaticBeam(double wavelength) : wavelength_{wavelength} {} + +MonochromaticBeam::MonochromaticBeam(Vector3d s0){ + double len = s0.norm(); + wavelength_ = 1.0 / len; + sample_to_source_direction_ = -1.0 * s0 / len; +} + +// constructor from json data +MonochromaticBeam::MonochromaticBeam(json beam_data) { + // minimal required keys + std::vector required_keys = {"wavelength"}; + for (const auto &key : required_keys) { + if (beam_data.find(key) == beam_data.end()) { + throw std::invalid_argument( + "Key " + key + " is missing from the input beam JSON"); + } + } + wavelength_ = beam_data["wavelength"]; + init_from_json(beam_data); +} + +// serialize to json format +json MonochromaticBeam::to_json(std::string probe="x-ray") const { + // create a json object that conforms to a dials model serialization. + json beam_data = {{"__id__", "monochromatic"}, {"probe", probe}}; + beam_data["wavelength"] = wavelength_; + add_to_json(beam_data); + return beam_data; +} + +double MonochromaticBeam::get_wavelength() const { + return wavelength_; +} +void MonochromaticBeam::set_wavelength(double wavelength){ + wavelength_ = wavelength; +} + +Vector3d MonochromaticBeam::get_s0() const { + return -sample_to_source_direction_ / wavelength_; +} +void MonochromaticBeam::set_s0(Vector3d s0){ + double len = s0.norm(); + wavelength_ = 1.0 / len; + sample_to_source_direction_ = -1.0 * s0 / len; +} + + +// PolychromaticBeam definitions + +// full constructor +PolychromaticBeam::PolychromaticBeam( + std::array wavelength_range, Vector3d direction, double divergence, + double sigma_divergence, Vector3d polarization_normal, + double polarization_fraction, double flux, + double transmission, double sample_to_source_distance) + : BeamBase{direction, divergence, sigma_divergence, + polarization_normal, polarization_fraction, + flux, transmission, sample_to_source_distance}, wavelength_range_{wavelength_range} {} + +PolychromaticBeam::PolychromaticBeam(std::array wavelength_range) : wavelength_range_{wavelength_range} {} + +PolychromaticBeam::PolychromaticBeam(std::array wavelength_range, Vector3d direction) + : wavelength_range_{wavelength_range} { + sample_to_source_direction_ = direction / direction.norm(); +} + +// constructor from json data +PolychromaticBeam::PolychromaticBeam(json beam_data) { + // minimal required keys + std::vector required_keys = {"wavelength_range"}; + for (const auto &key : required_keys) { + if (beam_data.find(key) == beam_data.end()) { + throw std::invalid_argument( + "Key " + key + " is missing from the input beam JSON"); + } + } + wavelength_range_ = beam_data["wavelength_range"]; + init_from_json(beam_data); +} + +// serialize to json format +json PolychromaticBeam::to_json(std::string probe="x-ray") const { + // create a json object that conforms to a dials model serialization. + json beam_data = {{"__id__", "polychromatic"}, {"probe", probe}}; + beam_data["wavelength_range"] = wavelength_range_; + add_to_json(beam_data); + return beam_data; +} + +std::array PolychromaticBeam::get_wavelength_range() const { + return wavelength_range_; +} +void PolychromaticBeam::set_wavelength_range(std::array wavelength_range){ + wavelength_range_ = wavelength_range; +} + + +// MonoXrayBeam definitions + +json MonoXrayBeam::to_json() const { + // call the parent function with the correct probe name (which is the same as the default in this case) + return MonochromaticBeam::to_json("x-ray"); +} + +// MonoElectronBeam definitions + +json MonoElectronBeam::to_json() const { + // call the parent function with the correct probe name + return MonochromaticBeam::to_json("electron"); +} + +#endif //DX2_MODEL_BEAM_H diff --git a/include/dx2/crystal.h b/include/dx2/crystal.h new file mode 100644 index 0000000..af0642c --- /dev/null +++ b/include/dx2/crystal.h @@ -0,0 +1,137 @@ +#ifndef DX2_MODEL_CRYSTAL_H +#define DX2_MODEL_CRYSTAL_H +#include +#include +#include +#include +#include +#include +#include +#include +#include "utils.h" + +using json = nlohmann::json; +using Eigen::Matrix3d; +using Eigen::Vector3d; + +Matrix3d Matrix3d_from_gemmi_cb(gemmi::Op cb){ + std::array, 3> rot = cb.rot; + return Matrix3d{ + {(double)(rot[0][0]/cb.DEN), (double)(rot[1][0]/cb.DEN), (double)(rot[2][0]/cb.DEN)}, + {(double)(rot[0][1]/cb.DEN), (double)(rot[1][1]/cb.DEN), (double)(rot[2][1]/cb.DEN)}, + {(double)(rot[0][2]/cb.DEN), (double)(rot[1][2]/cb.DEN),(double)(rot[2][2]/cb.DEN)}}; +} + + + +class Crystal { + public: + Crystal()=default; + Crystal(Vector3d a, Vector3d b, Vector3d c, gemmi::SpaceGroup space_group); + Crystal(json crystal_data); + gemmi::UnitCell get_unit_cell() const; + gemmi::SpaceGroup get_space_group() const; + Matrix3d get_A_matrix() const; + Matrix3d get_U_matrix() const; + Matrix3d get_B_matrix() const; + void niggli_reduce(); + json to_json() const; + + protected: + void init_from_abc(Vector3d a, Vector3d b, Vector3d c); + gemmi::SpaceGroup space_group_; + gemmi::UnitCell unit_cell_; + Matrix3d B_; + Matrix3d A_; + Matrix3d U_; + +}; + +void Crystal::init_from_abc(Vector3d a, Vector3d b, Vector3d c){ + // calculate B matrix, A matrix, set the input cell values + Matrix3d A{{a[0], a[1], a[2]}, {b[0], b[1], b[2]}, {c[0], c[1], c[2]}}; + A_ = A.inverse(); + double real_space_a = a.norm(); + double real_space_b = b.norm(); + double real_space_c = c.norm(); + double alpha = angle_between_vectors_degrees(b, c); + double beta = angle_between_vectors_degrees(c, a); + double gamma = angle_between_vectors_degrees(a, b); + unit_cell_ = {real_space_a,real_space_b,real_space_c,alpha,beta,gamma}; + gemmi::Mat33 B = unit_cell_.frac.mat; + Matrix3d B_{{B.a[0][0], B.a[0][1], B.a[0][2]}, {B.a[1][0], B.a[1][1], B.a[1][2]}, {B.a[2][0], B.a[2][1], B.a[2][2]}}; + U_ = A_ * B_; +} + +Crystal::Crystal(Vector3d a, Vector3d b, Vector3d c, gemmi::SpaceGroup space_group): space_group_(space_group) { + init_from_abc(a,b,c); +} + +Crystal::Crystal(json crystal_data){ + std::vector required_keys = { + "real_space_a", "real_space_b", "real_space_c", "space_group_hall_symbol" + }; + for (const auto &key : required_keys) { + if (crystal_data.find(key) == crystal_data.end()) { + throw std::invalid_argument( + "Key " + key + " is missing from the input crystal JSON"); + } + } + Vector3d rsa{{crystal_data["real_space_a"][0], crystal_data["real_space_a"][1],crystal_data["real_space_a"][2]}}; + Vector3d rsb{{crystal_data["real_space_b"][0], crystal_data["real_space_b"][1],crystal_data["real_space_b"][2]}}; + Vector3d rsc{{crystal_data["real_space_c"][0], crystal_data["real_space_c"][1],crystal_data["real_space_c"][2]}}; + space_group_ = *gemmi::find_spacegroup_by_name(crystal_data["space_group_hall_symbol"]); + init_from_abc(rsa, rsb, rsc); +} + +void Crystal::niggli_reduce() { + char centering{'P'}; + gemmi::GruberVector gv(unit_cell_, centering, true); + gv.niggli_reduce(); + gemmi::UnitCell niggli_cell = gv.get_cell(); + std::cout << "Input cell:" << std::endl; + std::cout << unit_cell_.a << " " << unit_cell_.b << " " << unit_cell_.c << " " << unit_cell_.alpha << " " << unit_cell_.beta << " " << unit_cell_.gamma << std::endl; + std::cout << "Reduced cell:" << std::endl; + std::cout << niggli_cell.a << " " << niggli_cell.b << " " << niggli_cell.c << " " << niggli_cell.alpha << " " << niggli_cell.beta << " " << niggli_cell.gamma << std::endl; + unit_cell_ = niggli_cell; + gemmi::Op cb = *gv.change_of_basis; + + Matrix3d cb_op = Matrix3d_from_gemmi_cb(cb); + A_ = A_*cb_op.inverse(); +} + +gemmi::UnitCell Crystal::get_unit_cell() const { + return unit_cell_; +} + +gemmi::SpaceGroup Crystal::get_space_group() const { + return space_group_; +} + +Matrix3d Crystal::get_A_matrix() const { + return A_; +} + +Matrix3d Crystal::get_B_matrix() const { + return B_; +} + +Matrix3d Crystal::get_U_matrix() const { + return U_; +} + +json Crystal::to_json() const { + json crystal_data; + crystal_data["__id__"] = "crystal"; + Matrix3d A_inv = A_.inverse(); + Vector3d rsa{{A_inv(0,0), A_inv(0,1), A_inv(0,2)}}; + Vector3d rsb{{A_inv(1,0), A_inv(1,1), A_inv(1,2)}}; + Vector3d rsc{{A_inv(2,0), A_inv(2,1), A_inv(2,2)}}; + crystal_data["real_space_a"] = rsa; + crystal_data["real_space_b"] = rsb; + crystal_data["real_space_c"] = rsc; + crystal_data["space_group_hall_symbol"] = space_group_.hall; + return crystal_data; +} + +#endif // DX2_MODEL_CRYSTAL_H diff --git a/include/dx2/detector.h b/include/dx2/detector.h new file mode 100644 index 0000000..6a657e6 --- /dev/null +++ b/include/dx2/detector.h @@ -0,0 +1,178 @@ +#ifndef DX2_MODEL_DETECTOR_H +#define DX2_MODEL_DETECTOR_H +#include +#include +#include + +using json = nlohmann::json; +using Eigen::Matrix3d; +using Eigen::Vector3d; + +double attenuation_length(double mu, double t0, + Vector3d s1, + Vector3d fast, + Vector3d slow, + Vector3d origin) { + Vector3d normal = fast.cross(slow); + double distance = origin.dot(normal); + if (distance < 0) { + normal = -normal; + } + double cos_t = s1.dot(normal); + //DXTBX_ASSERT(mu > 0 && cos_t > 0); + return (1.0 / mu) - (t0 / cos_t + 1.0 / mu) * exp(-mu * t0 / cos_t); +} + +class Panel { +// A class to represent a single "panel" of a detector (i.e. what data are +// considered to be described by a single set of panel parameters for the +// purposes of data processing, which may consist of several real detector +// modules). +public: + Panel()=default; + Panel(json panel_data); + Matrix3d get_d_matrix() const; + std::array px_to_mm(double x, double y) const; + std::array get_ray_intersection(Vector3d s1) const; + json to_json() const; +protected: + // panel_frame items + Vector3d origin_{{0.0,0.0,100.0}}; //needs to be set + Vector3d fast_axis_{{1.0,0.0,0.0}}; + Vector3d slow_axis_{{0.0,1.0,0.0}}; + Vector3d normal_{{0.0,0.0,1.0}}; + Matrix3d d_{{1,0,0},{0,1,0},{0,0,100.0}}; + Matrix3d D_{{1,0,0},{0,1,0},{0,0,0.01}}; + //double distance_{100.0}; + // panel data + std::array pixel_size_{{0.075,0.075}}; + std::array image_size_{{0,0}}; + std::array trusted_range_{0.0,65536.0}; + std::string type_{"SENSOR_PAD"}; + std::string name_{"module"}; + // also identifier and material present in dxtbx serialization. + double thickness_{0.0}; + double mu_{0.0}; + std::array raw_image_offset_{{0,0}}; //what's this? + // also mask would usually be here - is this what we want still? + //panel + double gain_{1.0}; + double pedestal_{0.0}; + std::string pixel_to_mm_strategy_{"SimplePxMmStrategy"}; // just the name here + bool parallax_correction_ = false; +}; + + +Panel::Panel(json panel_data){ + Vector3d fast{{panel_data["fast_axis"][0], panel_data["fast_axis"][1], panel_data["fast_axis"][2]}}; + Vector3d slow{{panel_data["slow_axis"][0], panel_data["slow_axis"][1], panel_data["slow_axis"][2]}}; + Vector3d origin{{panel_data["origin"][0], panel_data["origin"][1], panel_data["origin"][2]}}; + Matrix3d d_matrix{{fast[0], slow[0], origin[0]},{fast[1], slow[1], origin[1]},{fast[2], slow[2], origin[2]}}; + origin_ = origin; + fast_axis_ = fast; + slow_axis_ = slow; + normal_ = fast_axis_.cross(slow_axis_); + d_ = d_matrix; + D_ = d_.inverse(); + pixel_size_ = {{panel_data["pixel_size"][0], panel_data["pixel_size"][1]}}; + image_size_ = {{panel_data["image_size"][0], panel_data["image_size"][1]}}; + trusted_range_ = {{panel_data["trusted_range"][0], panel_data["trusted_range"][1]}}; + type_ = panel_data["type"]; + name_ = panel_data["name"]; + thickness_ = panel_data["thickness"]; + mu_ = panel_data["mu"]; + raw_image_offset_ = {{panel_data["raw_image_offset"][0], panel_data["raw_image_offset"][1]}}; + gain_ = panel_data["gain"]; + pedestal_ = panel_data["pedestal"]; + pixel_to_mm_strategy_ = panel_data["px_mm_strategy"]["type"]; + if (pixel_to_mm_strategy_ != std::string("SimplePxMmStrategy")){ + parallax_correction_ = true; + } +} + +json Panel::to_json() const { + json panel_data; + panel_data["name"] = name_; + panel_data["type"] = type_; + panel_data["fast_axis"] = fast_axis_; + panel_data["slow_axis"] = slow_axis_; + panel_data["origin"] = origin_; + panel_data["raw_image_offset"] = raw_image_offset_; + panel_data["image_size"] = image_size_; + panel_data["pixel_size"] = pixel_size_; + panel_data["trusted_range"] = trusted_range_; + panel_data["thickness"] = thickness_; + panel_data["mu"] = mu_; + panel_data["mask"] = std::array {}; + panel_data["identifier"]= ""; + panel_data["gain"] = gain_; + panel_data["pedestal"] = pedestal_, + panel_data["px_mm_strategy"] = {{"type", "ParallaxCorrectedPxMmStrategy"}}; + return panel_data; +} + +Matrix3d Panel::get_d_matrix() const { + return d_; +} + +std::array Panel::get_ray_intersection(Vector3d s1) const { + Vector3d v = D_ * s1; + // assert v[2] > 0 + std::array pxy; + pxy[0] = v[0] / v[2]; + pxy[1] = v[1] / v[2]; + // FIXME check is valid + return pxy; //in mmm +} + +std::array Panel::px_to_mm(double x, double y) const { + double x1 = x*pixel_size_[0]; + double x2 = y*pixel_size_[1]; + if (!parallax_correction_){ + return std::array{x1,x2}; + } + Vector3d fast = d_.col(0); + Vector3d slow = d_.col(1); + Vector3d origin = d_.col(2); + Vector3d s1 = origin + x1*fast + x2*slow; + s1.normalize(); + double o = attenuation_length(mu_, thickness_, s1, fast, slow, origin); + double c1 = x1 - (s1.dot(fast))*o; + double c2 = x2 - (s1.dot(slow))*o; + return std::array{c1,c2}; +} + +// Define a simple detector, for now is just a vector of panels without any hierarchy. +class Detector { + public: + Detector()=default; + Detector(json detector_data); + json to_json() const; + std::vector panels() const; + + protected: + std::vector _panels{}; +}; + +Detector::Detector(json detector_data){ + json panel_data = detector_data["panels"]; + for (json::iterator it = panel_data.begin(); it !=panel_data.end(); ++it) { + _panels.push_back(Panel(*it)); + } +} + +json Detector::to_json() const { + json detector_data; + std::vector panels_array; + for (auto p = _panels.begin(); p != _panels.end(); ++p){ + panels_array.push_back(p->to_json()); + } + detector_data["panels"] = panels_array; + return detector_data; +} + +std::vector Detector::panels() const { + return _panels; +} + +#endif // DX2_MODEL_DETECTOR_H \ No newline at end of file diff --git a/include/dx2/experiment.h b/include/dx2/experiment.h new file mode 100644 index 0000000..40d65eb --- /dev/null +++ b/include/dx2/experiment.h @@ -0,0 +1,122 @@ +#ifndef DX2_MODEL_EXPERIMENT_H +#define DX2_MODEL_EXPERIMENT_H +#include +#include +#include +#include +#include +#include +#include + +using Eigen::Vector3d; +using json = nlohmann::json; + +template +class Experiment { + public: + Experiment()=default; + Experiment(json experiment_data); + json to_json() const; + Goniometer goniometer() const; + BeamType beam() const; + Scan scan() const; + Detector detector() const; + Crystal crystal() const; + void set_crystal(Crystal crystal); + + protected: + BeamType _beam{}; + Scan _scan{}; + Goniometer _goniometer{}; + Detector _detector{}; + Crystal _crystal{}; +}; + +template +Experiment::Experiment(json experiment_data){ + json beam_data = experiment_data["beam"][0]; + BeamType beam = BeamType(beam_data); + json scan_data = experiment_data["scan"][0]; + Scan scan(scan_data); + json gonio_data = experiment_data["goniometer"][0]; + Goniometer gonio(gonio_data); + json detector_data = experiment_data["detector"][0]; + Detector detector(detector_data); + this->_beam = beam; + this->_scan = scan; + this->_goniometer = gonio; + this->_detector = detector; + try { // We don't always have a crystal model e.g. before indexing. + json crystal_data = experiment_data["crystal"][0]; + Crystal crystal(crystal_data); + this->_crystal = crystal; + } + catch (...) { + ; + } +} + +template +json Experiment::to_json() const { + // save this experiment as an example experiment list + json elist_out; // a list of potentially multiple experiments + elist_out["__id__"] = "ExperimentList"; + json expt_out; // our single experiment + // no imageset (for now?). + expt_out["__id__"] = "Experiment"; + expt_out["identifier"] = "test"; + expt_out["beam"] = + 0; // the indices of the models that will correspond to our experiment + expt_out["detector"] = 0; + expt_out["goniometer"] = 0; + expt_out["scan"] = 0; + + // add the the actual models + elist_out["scan"] = std::array{_scan.to_json()}; + elist_out["goniometer"] = std::array{_goniometer.to_json()}; + elist_out["beam"] = std::array{_beam.to_json()}; + elist_out["detector"] = std::array{_detector.to_json()}; + + if (_crystal.get_U_matrix().determinant()){ + expt_out["crystal"] = 0; + elist_out["crystal"] = std::array{_crystal.to_json()}; + } + else { + elist_out["crystal"] = std::array{}; + } + + elist_out["experiment"] = std::array{expt_out}; + return elist_out; +} + +template +Scan Experiment::scan() const { + return _scan; +} + +template +Goniometer Experiment::goniometer() const { + return _goniometer; +} + +template +Detector Experiment::detector() const { + return _detector; +} + +template +Crystal Experiment::crystal() const { + return _crystal; +} + +template +void Experiment::set_crystal(Crystal crystal) { + _crystal = crystal; +} + +template +BeamType Experiment::beam() const { + return _beam; +} + +#endif //DX2_MODEL_EXPERIMENT_H \ No newline at end of file diff --git a/include/dx2/goniometer.h b/include/dx2/goniometer.h new file mode 100644 index 0000000..5fd8f32 --- /dev/null +++ b/include/dx2/goniometer.h @@ -0,0 +1,161 @@ +#ifndef DX2_MODEL_GONIO_H +#define DX2_MODEL_GONIO_H +#include +#include +#include + +using json = nlohmann::json; +using Eigen::Matrix3d; +using Eigen::Vector3d; + +Matrix3d axis_and_angle_as_matrix(Vector3d axis, double angle, bool deg=false){ + double q0=0.0; + double q1=0.0; + double q2=0.0; + double q3=0.0; + if (deg) { + angle *= M_PI / 180.0; + } + if (!(std::fmod(angle, 2.0*M_PI))){ + q0=1.0; + } + else { + double h = 0.5 * angle; + q0 = std::cos(h); + double s = std::sin(h); + Vector3d n = axis / axis.norm(); + q1 = n[0]*s; + q2 = n[1]*s; + q3 = n[2]*s; + } + Matrix3d m{ + {2*(q0*q0+q1*q1)-1, 2*(q1*q2-q0*q3), 2*(q1*q3+q0*q2)}, + {2*(q1*q2+q0*q3), 2*(q0*q0+q2*q2)-1, 2*(q2*q3-q0*q1)}, + {2*(q1*q3-q0*q2), 2*(q2*q3+q0*q1), 2*(q0*q0+q3*q3)-1}}; + return m; +} + +class Goniometer { +// A class to represent a multi-axis goniometer. +public: + Goniometer()=default; + Goniometer( + std::vector axes, + std::vector angles, + std::vector names, + std::size_t scan_axis); + Goniometer(json goniometer_data); + Matrix3d get_setting_rotation() const; + Matrix3d get_sample_rotation() const; + Vector3d get_rotation_axis() const; + json to_json() const; + +protected: + void init(); // Sets the matrices from the axes and angles + Matrix3d calculate_setting_rotation(); + Matrix3d calculate_sample_rotation(); + Matrix3d sample_rotation_{{1,0,0},{0,1,0},{0,0,1}}; // F + Vector3d rotation_axis_{{1.0,0.0,0.0}}; // R' + Matrix3d setting_rotation_{{1,0,0},{0,1,0},{0,0,1}}; // S + std::vector axes_{{1.0,0.0,0.0}}; + std::vector angles_{{0.0}}; + std::vector names_{"omega"}; + std::size_t scan_axis_{0}; +}; + +void Goniometer::init(){ + // Sets the matrices from the axes and angles + setting_rotation_ = calculate_setting_rotation(); + sample_rotation_ = calculate_sample_rotation(); + rotation_axis_ = axes_[scan_axis_]; +} + +Matrix3d Goniometer::calculate_setting_rotation(){ + Matrix3d setting_rotation{{1,0,0},{0,1,0},{0,0,1}}; + for (std::size_t i = scan_axis_ + 1; i < axes_.size(); i++) { + Matrix3d R = axis_and_angle_as_matrix(axes_[i], angles_[i], true); + setting_rotation = R * setting_rotation; + } + return setting_rotation; +} + +Matrix3d Goniometer::calculate_sample_rotation(){ + Matrix3d sample_rotation{{1,0,0},{0,1,0},{0,0,1}}; + for (std::size_t i = 0; i < scan_axis_; i++) { + Matrix3d R = axis_and_angle_as_matrix(axes_[i], angles_[i], true); + sample_rotation = R * sample_rotation; + } + return sample_rotation; +} + +Goniometer::Goniometer( + std::vector axes, + std::vector angles, + std::vector names, + std::size_t scan_axis) + : axes_{axes.begin(), axes.end()}, + angles_{angles.begin(), angles.end()}, + names_{names.begin(), names.end()}, + scan_axis_{scan_axis}{ + if (scan_axis_ >= axes_.size()){ + throw std::invalid_argument("Goniometer scan axis number is out of range of axis length"); + } + init(); +} + +Matrix3d Goniometer::get_setting_rotation() const{ + return setting_rotation_; +} + +Matrix3d Goniometer::get_sample_rotation() const { + return sample_rotation_; +} + +Vector3d Goniometer::get_rotation_axis() const { + return rotation_axis_; +} + +Goniometer::Goniometer(json goniometer_data){ + std::vector required_keys = { + "axes", "angles", "names", "scan_axis" + }; + for (const auto &key : required_keys) { + if (goniometer_data.find(key) == goniometer_data.end()) { + throw std::invalid_argument( + "Key " + key + " is missing from the input goniometer JSON"); + } + } + std::vector axes; + for (json::iterator it=goniometer_data["axes"].begin();it != goniometer_data["axes"].end();it++){ + Vector3d axis; + axis[0] = (*it)[0]; + axis[1] = (*it)[1]; + axis[2] = (*it)[2]; + axes.push_back(axis); + } + axes_ = axes; + std::vector angles; + for (json::iterator it=goniometer_data["angles"].begin();it != goniometer_data["angles"].end();it++){ + angles.push_back(*it); + } + angles_ = angles; + std::vector names; + for (json::iterator it=goniometer_data["names"].begin();it != goniometer_data["names"].end();it++){ + names.push_back(*it); + } + names_ = names; + scan_axis_ = goniometer_data["scan_axis"]; + init(); +} + +json Goniometer::to_json() const { + json goniometer_data; + goniometer_data["axes"] = axes_; + goniometer_data["angles"] = angles_; + goniometer_data["names"] = names_; + goniometer_data["scan_axis"] = scan_axis_; + return goniometer_data; +} + + +#endif //DX2_MODEL_GONIO_H diff --git a/include/dx2/h5read_processed.h b/include/dx2/h5read_processed.h new file mode 100644 index 0000000..99228d9 --- /dev/null +++ b/include/dx2/h5read_processed.h @@ -0,0 +1,38 @@ +#include +#include +#include +#include +#include +#include +#include +#include + +template +std::vector read_array_from_h5_file(std::string filename, std::string array_name){ + auto start_time = std::chrono::high_resolution_clock::now(); + hid_t file = H5Fopen(filename.c_str(), H5F_ACC_RDONLY, H5P_DEFAULT); + if (file < 0){ + std::cout << "Error: Unable to open " << filename.c_str() << " as a hdf5 reflection table" < data_out(num_elements); + + H5Dread(dataset, datatype, H5S_ALL, space, H5P_DEFAULT, &data_out[0]); + float total_time = + std::chrono::duration_cast>( + std::chrono::high_resolution_clock::now() - start_time) + .count(); + std::cout << "READ TIME for " << array_name << " : " << total_time << "s" << std::endl; + + H5Dclose(dataset); + H5Fclose(file); + return data_out; +} diff --git a/include/dx2/scan.h b/include/dx2/scan.h new file mode 100644 index 0000000..1008b44 --- /dev/null +++ b/include/dx2/scan.h @@ -0,0 +1,84 @@ +#ifndef DX2_MODEL_SCAN_H +#define DX2_MODEL_SCAN_H +#include +using json = nlohmann::json; + +class Scan { +// A class to represent the physical measurement, consisting of the number of images, +// starting oscillation and a constant oscillation width between sequential images. +// This class MUST NOT be modified during processing or used to track additional metadata. +public: + Scan()=default; + Scan(std::array image_range, std::array oscillation); + Scan(json scan_data); + std::array get_image_range() const; + std::array get_oscillation() const; + json to_json() const; + +protected: + std::array image_range_{{0,0}}; + int num_images_{0}; + double oscillation_width_{0.0}; + double oscillation_start_{0.0}; +}; + +Scan::Scan(std::array image_range, std::array oscillation) + : image_range_{image_range}{ + num_images_ = image_range_[1] - image_range_[0]+1; + oscillation_start_ = oscillation[0]; + oscillation_width_ = oscillation[1]; +} + +Scan::Scan(json scan_data){ + // minimal required keys are image range and ["properties"]:"oscillation" + std::vector required_keys = { + "image_range", "properties" + }; + for (const auto &key : required_keys) { + if (scan_data.find(key) == scan_data.end()) { + throw std::invalid_argument( + "Key " + key + " is missing from the input scan JSON"); + } + } + image_range_ = scan_data["image_range"]; + num_images_ = image_range_[1] - image_range_[0]+1; + if (scan_data["properties"].find("oscillation") == scan_data["properties"].end()){ + throw std::invalid_argument( + "Key oscillation is missing from the input scan['properties'] JSON"); + } + std::vector oscillation; + json osc_array = scan_data["properties"]["oscillation"]; + // Just read the first two oscillation values as that is all that is needed + if (osc_array.size() < 2){ + throw std::invalid_argument( + "scan['properties']['oscillation'] has <2 values in the scan JSON"); + } + for (json::iterator it = osc_array.begin(); it != osc_array.begin()+2; ++it) { + oscillation.push_back(*it); + } + oscillation_start_ = oscillation[0]; + oscillation_width_ = oscillation[1] - oscillation[0]; +} + +json Scan::to_json() const { + json scan_data; + scan_data["image_range"] = image_range_; + scan_data["batch_offset"] = 0; // We MUST NOT use batch offsets in dx2, written out here for compatibility + std::vector oscillation_array(num_images_); + for (int i=0; i Scan::get_image_range() const { + return image_range_; +} + +std::array Scan::get_oscillation() const { + return {oscillation_start_, oscillation_width_}; +} + +#endif //DX2_MODEL_SCAN_H diff --git a/include/dx2/utils.h b/include/dx2/utils.h new file mode 100644 index 0000000..ce6cd6b --- /dev/null +++ b/include/dx2/utils.h @@ -0,0 +1,23 @@ +#ifndef DX2_UTILS_H +#define DX2_UTILS_H +#include +#include + +using Eigen::Vector3d; + +double angle_between_vectors_degrees(Vector3d v1, Vector3d v2) { + double l1 = v1.norm(); + double l2 = v2.norm(); + double dot = v1.dot(v2); + double normdot = dot / (l1 * l2); + if (std::abs(normdot - 1.0) < 1E-6) { + return 0.0; + } + if (std::abs(normdot + 1.0) < 1E-6) { + return 180.0; + } + double angle = std::acos(normdot) * 180.0 / M_PI; + return angle; +} + +#endif // DX2_UTILS_H diff --git a/requirements.txt b/requirements.txt new file mode 100644 index 0000000..257e3fb --- /dev/null +++ b/requirements.txt @@ -0,0 +1,4 @@ +gemmi +cmake +hdf5=1.12.2 +hdf5-external-filter-plugins diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index af440d4..2de1601 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -1,6 +1,11 @@ -add_executable(test_example test_example.cxx) -target_link_libraries(test_example GTest::gtest_main) +add_executable(test_make_models test_make_models.cxx) +target_link_libraries(test_make_models GTest::gtest_main dx2 Eigen3::Eigen nlohmann_json::nlohmann_json gemmi::gemmi_cpp) + +add_executable(test_read_h5_array test_read_h5_array.cxx) +target_link_libraries(test_read_h5_array GTest::gtest_main dx2 hdf5::hdf5) include(GoogleTest) -gtest_discover_tests(test_example) \ No newline at end of file + +gtest_discover_tests(test_make_models PROPERTIES LABELS dx2tests) +gtest_discover_tests(test_read_h5_array WORKING_DIRECTORY "${PROJECT_SOURCE_DIR}/tests" PROPERTIES LABELS dx2tests) diff --git a/tests/data/cut_strong.refl b/tests/data/cut_strong.refl new file mode 100644 index 0000000..4a751da Binary files /dev/null and b/tests/data/cut_strong.refl differ diff --git a/tests/test_example.cxx b/tests/test_example.cxx deleted file mode 100644 index 44b19df..0000000 --- a/tests/test_example.cxx +++ /dev/null @@ -1,3 +0,0 @@ -#include - -TEST(ExampleTests, EmptyTest) { EXPECT_EQ(2 + 2, 4); } \ No newline at end of file diff --git a/tests/test_make_models.cxx b/tests/test_make_models.cxx new file mode 100644 index 0000000..1750e04 --- /dev/null +++ b/tests/test_make_models.cxx @@ -0,0 +1,16 @@ +#include +#include + +using Eigen::Vector3d; +using Eigen::Matrix3d; + +// A very basic test that we can make a crystal, tests tht gemmi & Eigen is installed correctly too. +TEST(ExampleTests, CrystalTest) { + Vector3d a = {1.0, 0.0, 0.0}; + Vector3d b = {0.0, 1.0, 0.0}; + Vector3d c = {0.0, 0.0, 2.0}; + gemmi::SpaceGroup space_group = *gemmi::find_spacegroup_by_name("P1"); + Crystal crystal(a,b,c,space_group); + Matrix3d expected_A{{1.0,0.0,0.0}, {0.0,1.0,0.0}, {0.0,0.0,0.5}}; + EXPECT_EQ(crystal.get_A_matrix(), expected_A); + } diff --git a/tests/test_read_h5_array.cxx b/tests/test_read_h5_array.cxx new file mode 100644 index 0000000..4935466 --- /dev/null +++ b/tests/test_read_h5_array.cxx @@ -0,0 +1,27 @@ +#include +#include +#include +#include + +// A test that we can read data arrays from a h5 processing file. + +TEST(ExampleTests, ReadArrayTest) { + + std::string array_name = "/dials/processing/group_0/xyzobs.px.value"; + std::string flags = "/dials/processing/group_0/flags"; + // The test has been given the tests directory as the working directory + // so we can get the path to the test file. + std::filesystem::path cwd = std::filesystem::current_path(); + std::string fpath = cwd.generic_string(); + fpath.append("/data/cut_strong.refl"); + + std::vector xyzobs_px = read_array_from_h5_file(fpath, array_name); + // check a random value + double expected_value = 528.86470588235295; + EXPECT_EQ(xyzobs_px[10], expected_value); + + std::vector flags_array = read_array_from_h5_file(fpath, flags); + // check a random value + std::size_t expected_flag_value = 32; + EXPECT_EQ(flags_array[5], expected_flag_value); + }