Mrc format header to lattice coordinate transform

Extracts a coordinate transformation into lattice coordinates from the mrc density
header, as well as the three dimensional grid extents.

refs #2282

Change-Id: I8c01a370e3010add29b8810e2f3320b9ab428c02

diff --git a/src/gromacs/fileio/mrcdensitymapheader.cpp b/src/gromacs/fileio/mrcdensitymapheader.cpp
index bb8a1483364860a56e477ebe2778bb1e1cb907f5..9244996fd8b8eeba97418dc6c1452dd89f981a3b 100644 (file)
--- a/src/gromacs/fileio/mrcdensitymapheader.cpp
+++ b/src/gromacs/fileio/mrcdensitymapheader.cpp
@@ -64,4 +64,41 @@ size_t numberOfExpectedDataItems(const MrcDensityMapHeader &header)
     return header.numColumnRowSection_[XX] * header.numColumnRowSection_[YY] * header.numColumnRowSection_[ZZ];
 }
 
+TranslateAndScale getCoordinateTransformationToLattice(const MrcDensityMapHeader &header)
+{
+    constexpr real c_AAtoNmConversion = 0.1;
+
+    RVec           scale = {
+        header.extent_[XX] / (header.cellLength_[XX] * c_AAtoNmConversion),
+        header.extent_[YY] / (header.cellLength_[YY] * c_AAtoNmConversion),
+        header.extent_[ZZ] / (header.cellLength_[ZZ] * c_AAtoNmConversion)
+    };
+    const RVec     emdbOrigin {
+        header.userDefinedFloat_[12], header.userDefinedFloat_[13], header.userDefinedFloat_[14]
+    };
+    RVec translation;
+    if (emdbOrigin[XX] == 0. && emdbOrigin[YY] == 0. && emdbOrigin[ZZ] == 0.)
+    {
+        translation = RVec {
+            -header.columnRowSectionStart_[XX] / scale[XX],
+            -header.columnRowSectionStart_[YY] / scale[YY],
+            -header.columnRowSectionStart_[ZZ] / scale[ZZ]
+        };
+    }
+    else
+    {
+        translation = {
+            -emdbOrigin[XX] * c_AAtoNmConversion,
+            -emdbOrigin[YY] * c_AAtoNmConversion,
+            -emdbOrigin[ZZ] * c_AAtoNmConversion
+        };
+    }
+    return {scale, translation};
+}
+
+dynamicExtents3D getDynamicExtents3D(const MrcDensityMapHeader &header)
+{
+    return {header.numColumnRowSection_[XX], header.numColumnRowSection_[YY], header.numColumnRowSection_[ZZ]};
+};
+
 } // namespace gmx

diff --git a/src/gromacs/fileio/mrcdensitymapheader.h b/src/gromacs/fileio/mrcdensitymapheader.h
index 46adf70818b320273d3cf03e87339a545747157a..6f4fb06c11827cfaaf7a9bb61c8c1587ae755c69 100644 (file)
--- a/src/gromacs/fileio/mrcdensitymapheader.h
+++ b/src/gromacs/fileio/mrcdensitymapheader.h
@@ -47,7 +47,9 @@
 #include <array>
 #include <vector>
 
+#include "gromacs/math/coordinatetransformation.h"
 #include "gromacs/math/vectypes.h"
+#include "gromacs/mdspan/extensions.h"
 
 namespace gmx
 {
@@ -166,5 +168,26 @@ struct MrcDensityMapHeader{
  */
 size_t numberOfExpectedDataItems(const MrcDensityMapHeader &header);
 
+/*! \brief Extract the transformation into lattice coordinates.
+ * \note Transformation into lattice coordinates is not treated uniformly
+ *       in different implementations for the mrc format,e.g., vmd, pymol and
+ *       chimera have different conventions. Following the vmd implementation here.
+ *
+ * In determining the density origin coordinates, explicit ORIGIN records
+ * (also called origin2k) in the user defined floats 13 - 15, corresponding to
+ * words 50,51 and 52 in the mrc header, precedence over ColumnRowSectionStart.
+ * Only if above values are zero, using the column, row and section start to
+ * determine the translation vector.
+ *
+ * \param[in] header from which the coordinate transformation is to be extracted
+ * \returns a functor that transforms real space coordinates into the lattice
+ */
+TranslateAndScale getCoordinateTransformationToLattice(const MrcDensityMapHeader &header);
+
+/*! \brief Extract the extents of the density data
+ * \param[in] header from which the extents are to be extracted
+ * \returns density data extents in three dimensions.
+ */
+dynamicExtents3D getDynamicExtents3D(const MrcDensityMapHeader &header);
 }      // namespace gmx
 #endif /* end of include guard: GMX_FILEIO_MRCDENSITYMAPHEADER_H */

diff --git a/src/gromacs/fileio/tests/mrcdensitymapheader.cpp b/src/gromacs/fileio/tests/mrcdensitymapheader.cpp
index 8359a84c652e3340a77c191b2d1e64894d5bfc26..8e90e99c2599117b6e4dd12e25aa6e9042ef761f 100644 (file)
--- a/src/gromacs/fileio/tests/mrcdensitymapheader.cpp
+++ b/src/gromacs/fileio/tests/mrcdensitymapheader.cpp
@@ -44,10 +44,14 @@
 
 #include "gromacs/fileio/mrcdensitymapheader.h"
 
+#include <array>
+#include <vector>
+
 #include <gmock/gmock.h>
 #include <gtest/gtest.h>
 
 #include "testutils/testasserts.h"
+#include "testutils/testmatchers.h"
 
 namespace gmx
 {
@@ -71,5 +75,66 @@ TEST(MrcDensityMapHeaderTest, DataSizeThrowsWhenInvalid)
     header.numColumnRowSection_ = {-1, 2, 3};
     EXPECT_THROW(numberOfExpectedDataItems(header), InternalError);
 }
+TEST(MrcDensityMapHeaderTest, GetsCorrectCoordinateTransformNoOriginGiven)
+{
+    MrcDensityMapHeader header;
+
+    header.extent_                = {100, 200, 300};
+    header.columnRowSectionStart_ = {50, 200, 0};
+    header.cellLength_            = {10, 20, 15};
+
+    std::array<RVec, 2>                 testVectors = {RVec {0., 0., 0.}, RVec {1., 1., 1.}};
+    getCoordinateTransformationToLattice(header)(testVectors);
+
+    std::vector<RVec>                   expectedVectors = {{-50, -200, 0}, {50, -100, 200}};
+    EXPECT_THAT(expectedVectors, testing::Pointwise(test::RVecEq(test::defaultFloatTolerance()), testVectors));
+}
+
+TEST(MrcDensityMapHeaderTest, GetsCorrectCoordinateTransformWithOriginDefined)
+{
+    MrcDensityMapHeader header;
+    header.userDefinedFloat_[12]  = 1.;
+    header.userDefinedFloat_[13]  = 2.;
+    header.userDefinedFloat_[14]  = 3.;
+    header.extent_                = {100, 200, 300};
+    // setting the columnRowSectionStart values that are to be ignored if userDefinedFloat_ is not zero
+    header.columnRowSectionStart_ = {50, 200, 0};
+    header.cellLength_            = {10, 20, 15};
+
+    std::array<RVec, 2>                 testVectors = {RVec {0., 0., 0.}, RVec {1., 1., 1.}};
+    getCoordinateTransformationToLattice(header)(testVectors);
+
+    std::vector<RVec>                   expectedVectors = {{-10, -20, -60}, {90, 80, 140}};
+    EXPECT_THAT(expectedVectors, testing::Pointwise(test::RVecEq(test::defaultFloatTolerance()), testVectors));
+}
+
+TEST(MrcDensityMapHeaderTest, GetsCorrectCoordinateTransformWithStartValues)
+{
+    MrcDensityMapHeader header;
+    header.userDefinedFloat_[12]  = 0;
+    header.userDefinedFloat_[13]  = 0;
+    header.userDefinedFloat_[14]  = 0;
+    header.extent_                = {100, 200, 300};
+    header.columnRowSectionStart_ = {50, 200, 0};
+    header.cellLength_            = {10, 20, 15};
+
+    std::array<RVec, 2>                 testVectors = {RVec {0., 0., 0.}, RVec {1., 1., 1.}};
+    getCoordinateTransformationToLattice(header)(testVectors);
+
+    std::vector<RVec>                   expectedVectors = {{-50, -200, 0}, {50, -100, 200}};
+    EXPECT_THAT(expectedVectors, testing::Pointwise(test::RVecEq(test::defaultFloatTolerance()), testVectors));
+}
+
+TEST(MrcDensityMapHeaderTest, GetsCorrectExtents)
+{
+    MrcDensityMapHeader header;
+    header.numColumnRowSection_ = {100, 200, 300};
+
+    const auto extents = getDynamicExtents3D(header);
+    std::array<std::ptrdiff_t, DIM> expectedExtents = {100, 200, 300};
+    EXPECT_EQ(expectedExtents[XX], extents.extent(XX));
+    EXPECT_EQ(expectedExtents[YY], extents.extent(YY));
+    EXPECT_EQ(expectedExtents[ZZ], extents.extent(ZZ));
+}
 
 } // namespace gmx