Introduce a new class and tests for affine coordinate transformations.
Allows to perform a matrix multiplication and shift to coordinates of the form A*x+t with a 3x3 real valued matrix A and a translation vector t.
#include <vector>
#include "gromacs/math/vec.h"
+#include "gromacs/mdspan/extensions.h"
#include "gromacs/utility/arrayref.h"
+#include "matrix.h"
+
namespace gmx
{
+
/********************************************************************
* ScaleCoordinates::Impl
*/
TranslateAndScale& TranslateAndScale::operator=(TranslateAndScale&&) noexcept = default;
+/********************************************************************
+ * AffineTransformation
+ */
+
+AffineTransformation::AffineTransformation(Matrix3x3ConstSpan matrix, const RVec& translation) :
+ translation_{ translation }
+{
+ std::copy(begin(matrix), end(matrix), begin(matrix_));
+}
+
+void AffineTransformation::operator()(ArrayRef<RVec> vectors) const
+{
+ for (RVec& vector : vectors)
+ {
+ matrixVectorMultiply(matrix_.asConstView(), &vector);
+ vector += translation_;
+ }
+}
+
+void AffineTransformation::operator()(RVec* vector) const
+{
+ (*this)({ vector, vector + 1 });
+}
} // namespace gmx
#include "gromacs/math/vectypes.h"
#include "gromacs/utility/classhelpers.h"
+#include "matrix.h"
+
namespace gmx
{
class Impl;
PrivateImplPointer<Impl> impl_;
};
+
+/*! \libinternal
+ * \brief Affine transformation of three-dimensional coordinates.
+ *
+ * Perfoms in-place coordinate transformations.
+ *
+ * Coordinates are first multiplied by a matrix, then translated.
+ */
+class AffineTransformation
+{
+public:
+ /*! \brief Construct a three-dimensional affine transformation.
+ * \param[in] matrix to be applied to the vectors
+ * \param[in] translation to be performed on the vectors
+ */
+ AffineTransformation(Matrix3x3ConstSpan matrix, const RVec& translation);
+
+ /*! \brief Perform an affine transformation on input vectors.
+ * \param[in,out] vectors to be transformed in-place
+ */
+ void operator()(ArrayRef<RVec> vectors) const;
+
+ /*! \brief Perform an affine transformation on a vector.
+ * \param[in,out] vector to be transformed in-place
+ */
+ void operator()(RVec* vector) const;
+
+private:
+ //! The matrix describing the affine transformation A(x) = matrix_ * x + translation_
+ Matrix3x3 matrix_;
+ //! The translation vector describing the affine transformation A(x) = matrix * x + translation
+ RVec translation_;
+};
+
} // namespace gmx
#endif // CoordinateTransformation
--- /dev/null
+/*
+ * This file is part of the GROMACS molecular simulation package.
+ *
+ * Copyright (c) 2020, by the GROMACS development team, led by
+ * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
+ * and including many others, as listed in the AUTHORS file in the
+ * top-level source directory and at http://www.gromacs.org.
+ *
+ * GROMACS is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU Lesser General Public License
+ * as published by the Free Software Foundation; either version 2.1
+ * of the License, or (at your option) any later version.
+ *
+ * GROMACS is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+ * Lesser General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public
+ * License along with GROMACS; if not, see
+ * http://www.gnu.org/licenses, or write to the Free Software Foundation,
+ * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+ *
+ * If you want to redistribute modifications to GROMACS, please
+ * consider that scientific software is very special. Version
+ * control is crucial - bugs must be traceable. We will be happy to
+ * consider code for inclusion in the official distribution, but
+ * derived work must not be called official GROMACS. Details are found
+ * in the README & COPYING files - if they are missing, get the
+ * official version at http://www.gromacs.org.
+ *
+ * To help us fund GROMACS development, we humbly ask that you cite
+ * the research papers on the package. Check out http://www.gromacs.org.
+ */
+/*! \internal \file
+ *
+ * \brief Implements routines in matrix.h .
+ *
+ * \author Christian Blau <blau@kth.se>
+ */
+
+#include "gmxpre.h"
+
+#include "matrix.h"
+
+namespace gmx
+{
+
+
+Matrix3x3 transpose(Matrix3x3ConstSpan matrixView)
+{
+
+ return Matrix3x3({ matrixView(0, 0), matrixView(1, 0), matrixView(2, 0), matrixView(0, 1),
+ matrixView(1, 1), matrixView(2, 1), matrixView(0, 2), matrixView(1, 2),
+ matrixView(2, 2) });
+}
+
+void matrixVectorMultiply(Matrix3x3ConstSpan matrix, RVec* v)
+{
+ const real resultXX =
+ matrix(XX, XX) * (*v)[XX] + matrix(XX, YY) * (*v)[YY] + matrix(XX, ZZ) * (*v)[ZZ];
+ const real resultYY =
+ matrix(YY, XX) * (*v)[XX] + matrix(YY, YY) * (*v)[YY] + matrix(YY, ZZ) * (*v)[ZZ];
+ (*v)[ZZ] = matrix(ZZ, XX) * (*v)[XX] + matrix(ZZ, YY) * (*v)[YY] + matrix(ZZ, ZZ) * (*v)[ZZ];
+ (*v)[XX] = resultXX;
+ (*v)[YY] = resultYY;
+}
+
+
+} // namespace gmx
/*
* This file is part of the GROMACS molecular simulation package.
*
- * Copyright (c) 2019, by the GROMACS development team, led by
+ * Copyright (c) 2019,2020, by the GROMACS development team, led by
* Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
* and including many others, as listed in the AUTHORS file in the
* top-level source directory and at http://www.gromacs.org.
#include <array>
#include "gromacs/math/multidimarray.h"
+#include "gromacs/math/utilities.h"
#include "gromacs/math/vec.h"
#include "gromacs/utility/real.h"
return matrixView(0, 0) + matrixView(1, 1) + matrixView(2, 2);
}
-//! Calculate the transpose of a 3x3 matrix, from its view
-static Matrix3x3 transpose(Matrix3x3ConstSpan matrixView)
+/*! \brief Create an identity matrix of ElementType with N * M elements.
+ *
+ * \tparam ElementType type of matrix elements
+ * \tparam N number of rows
+ * \tparam M number of columns, defaults to number of rows if not set
+ *
+ * \returns a matrix with values one where row equals column index and null
+ * where row does not equal column index
+ */
+template<typename ElementType, int N, int M = N>
+MultiDimArray<std::array<ElementType, N * M>, extents<N, M>> identityMatrix()
{
-
- return Matrix3x3({ matrixView(0, 0), matrixView(1, 0), matrixView(2, 0), matrixView(0, 1),
- matrixView(1, 1), matrixView(2, 1), matrixView(0, 2), matrixView(1, 2),
- matrixView(2, 2) });
+ std::array<ElementType, N * M> matrixEntries{};
+ MultiDimArray<std::array<ElementType, N * M>, extents<N, M>> idMatrix(matrixEntries);
+ for (int i = 0; i < std::min(N, M); i++)
+ {
+ idMatrix(i, i) = 1;
+ }
+ return idMatrix;
}
+//! Calculate the transpose of a 3x3 matrix, from its view
+Matrix3x3 transpose(Matrix3x3ConstSpan matrixView);
+
+//! Multiply matrix with vector.
+void matrixVectorMultiply(Matrix3x3ConstSpan matrix, RVec* v);
+
//! Create new matrix type from legacy type.
static inline Matrix3x3 createMatrix3x3FromLegacyMatrix(const matrix legacyMatrix)
{
#include <gtest/gtest.h>
#include "gromacs/utility/arrayref.h"
+#include "gromacs/utility/real.h"
#include "testutils/testasserts.h"
#include "testutils/testmatchers.h"
{ 3, -6, 2.5 } };
};
+class AffineTransformationTest : public ::testing::Test
+{
+protected:
+ std::vector<RVec> testVectors_ = { { 0, 0, 0 },
+ { 1, 0, 0 },
+ { 0, -1, -1 },
+ { 1e10, 1e1, 1e-2 },
+ { 3, -6, 2.5 } };
+};
+
+
TEST_F(TranslateAndScaleTest, identityTransformation)
{
TranslateAndScale translateAndScale(identityScale_, identityTranslation_);
EXPECT_REAL_EQ(0, test[ZZ]);
}
+TEST_F(AffineTransformationTest, identityTransformYieldsSameVectors)
+{
+ const AffineTransformation identityTransformation(identityMatrix<real, 3>(), { 0, 0, 0 });
+ for (const auto& vector : testVectors_)
+ {
+ RVec vectorTransformed = vector;
+ identityTransformation(&vectorTransformed);
+ EXPECT_REAL_EQ(vector[XX], vectorTransformed[XX]);
+ EXPECT_REAL_EQ(vector[YY], vectorTransformed[YY]);
+ EXPECT_REAL_EQ(vector[ZZ], vectorTransformed[ZZ]);
+ }
+}
+
+TEST_F(AffineTransformationTest, applyTransformationToVectors)
+{
+ const Matrix3x3 transformMatrix({ 0.1, 1, 0.1, 0.4, 1, 0.6, 0.7, 0.8, 0.9 });
+ const RVec transformVector = { 1, -1e5, 1e4 };
+
+ const AffineTransformation affineTransformation(transformMatrix, transformVector);
+
+ const std::vector<RVec> expectedResult = { { 1, -100'000, 10'000 },
+ { 1.1, -99999.6, 10000.7 },
+ { -0.1, -100'002, 9998.3 },
+ { 1e9, 3.9999e9, 7.00001e9 },
+ { -4.45, -100'003, 9'999.5 } };
+
+ auto expected = expectedResult.begin();
+ for (const auto& vector : testVectors_)
+ {
+ RVec vectorTransformed = vector;
+ affineTransformation(&vectorTransformed);
+ // need relaxed tolerance here, due to the number of operations involved
+ EXPECT_REAL_EQ_TOL((*expected)[XX], vectorTransformed[XX],
+ relativeToleranceAsFloatingPoint((*expected)[XX], 1e-5));
+ EXPECT_REAL_EQ_TOL((*expected)[YY], vectorTransformed[YY],
+ relativeToleranceAsFloatingPoint((*expected)[YY], 1e-5));
+ EXPECT_REAL_EQ_TOL((*expected)[ZZ], vectorTransformed[ZZ],
+ relativeToleranceAsFloatingPoint((*expected)[ZZ], 1e-5));
+ ++expected;
+ }
+}
} // namespace test
} // namespace gmx
/*
* This file is part of the GROMACS molecular simulation package.
*
- * Copyright (c) 2019, by the GROMACS development team, led by
+ * Copyright (c) 2019,2020, by the GROMACS development team, led by
* Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
* and including many others, as listed in the AUTHORS file in the
* top-level source directory and at http://www.gromacs.org.
#include <gtest/gtest.h>
#include "gromacs/math/vec.h"
+#include "gromacs/utility/real.h"
#include "testutils/testasserts.h"
}
}
+TEST_F(MatrixTest, IdentityMatrix)
+{
+ const auto realIdMatrix = identityMatrix<real, 2>();
+ EXPECT_REAL_EQ(realIdMatrix(0, 0), 1);
+ EXPECT_REAL_EQ(realIdMatrix(1, 1), 1);
+ EXPECT_REAL_EQ(realIdMatrix(0, 1), 0);
+ EXPECT_REAL_EQ(realIdMatrix(1, 0), 0);
+}
+
+TEST_F(MatrixTest, MatrixVectorMultiplication)
+{
+ const Matrix3x3 matrix({ 0.1, 1, 0.1, 0.4, 1, 0.6, 0.7, 0.8, 0.9 });
+ RVec vector(1, 2, 3);
+ matrixVectorMultiply(matrix, &vector);
+ EXPECT_REAL_EQ(2.4, vector[XX]);
+ EXPECT_REAL_EQ(4.2, vector[YY]);
+ EXPECT_REAL_EQ(5.0, vector[ZZ]);
+}
+
} // namespace
} // namespace test
biod.pnpi.spb.ru / alexxy / gromacs.git / commitdiff