attempt to use quaternions

CMakeLists.txt

@@ -168,6 +168,8 @@ if(MSVC)
 endif()
 find_package(ZLIB REQUIRED)
 
+find_package(Eigen3)
+
 include(FindFilesystem)
 list(APPEND CIFPP_REQUIRED_LIBRARIES ${STDCPPFS_LIBRARY})
 
@@ -279,6 +281,8 @@ target_include_directories(cifpp
 
 target_link_libraries(cifpp PUBLIC Threads::Threads ZLIB::ZLIB ${CIFPP_REQUIRED_LIBRARIES})
 
+target_link_libraries(cifpp PRIVATE Eigen3::Eigen)
+
 if(CMAKE_CXX_COMPILER_ID STREQUAL "AppleClang")
 	target_link_options(cifpp PRIVATE -undefined dynamic_lookup)
 endif(CMAKE_CXX_COMPILER_ID STREQUAL "AppleClang")

include/cif++/matrix.hpp

@@ -518,7 +518,9 @@ auto eigen(const matrix_expression<M> &mat, bool sort)
 				value_type h = ev[j] - ev[i];
 				value_type t;
 
-				if (std::fabs(a_ij) > 1.0e-12 * std::fabs(h))
+				if (h == 0 and a_ij == 0)
+					t = 1;
+				else if (std::fabs(a_ij) > 1.0e-12 * std::fabs(h))
 				{
 					value_type theta = 0.5 * h / a_ij;
 					t = 1.0 / (std::fabs(theta) + std::sqrt(1 + theta * theta));

include/cif++/symmetry.hpp

@@ -248,10 +248,12 @@ class transformation
 	transformation &operator=(const transformation &) = default;
 	transformation &operator=(transformation &&) = default;
 
-	point operator()(const point &pt) const
-	{
-		return m_rotation * pt + m_translation;
-	}
+	point operator()(const point &pt) const;
+
+	// point operator()(const point &pt) const
+	// {
+	// 	return m_rotation * pt + m_translation;
+	// }
 
 	friend transformation operator*(const transformation &lhs, const transformation &rhs);
 	friend transformation inverse(const transformation &t);
@@ -265,6 +267,7 @@ class transformation
 
   private:
 	matrix3x3<float> m_rotation;
+	quaternion m_q;
 	point m_translation;
 };
 

src/symmetry.cpp

@@ -129,11 +129,54 @@ transformation::transformation(const symop_data &data)
 	m_rotation(2, 1) = d[7];
 	m_rotation(2, 2) = d[8];
 
+	auto &&[ev, em] = eigen(m_rotation, false);
+
+	for (size_t i = 0; i < 3; ++i)
+	{
+		if (ev[i] != 1)
+			continue;
+		
+		auto tr = m_rotation(0, 0) + m_rotation(1, 1) + m_rotation(2, 2);
+		auto angle = std::acos((tr - 1) / 2.0f) / 2;
+
+		auto s = std::sin(angle);
+		auto c = std::cos(angle);
+
+		m_q = normalize(quaternion{
+			static_cast<float>(c),
+			static_cast<float>(s * em(0, i)),
+			static_cast<float>(s * em(1, i)),
+			static_cast<float>(s * em(2, i)) });
+
+		break;
+	}
+
 	m_translation.m_x = d[9] == 0 ? 0 : 1.0 * d[9] / d[10];
 	m_translation.m_y = d[11] == 0 ? 0 : 1.0 * d[11] / d[12];
 	m_translation.m_z = d[13] == 0 ? 0 : 1.0 * d[13] / d[14];
 }
 
+point transformation::operator()(const point &pt) const
+{
+	auto p = pt;
+	p.rotate(m_q);
+	p += m_translation;
+
+	auto p2 = m_rotation * pt + m_translation;
+
+	// return m_rotation * pt + m_translation;
+
+	assert(std::abs(p.m_x - p2.m_x) < 0.01);
+	assert(std::abs(p.m_y - p2.m_y) < 0.01);
+	assert(std::abs(p.m_z - p2.m_z) < 0.01);
+
+	return p;
+
+	// auto p = pt;
+	// p.rotate(m_q);
+	// return p + m_translation;
+}
+
 transformation operator*(const transformation &lhs, const transformation &rhs)
 {
 	auto r = lhs.m_rotation * rhs.m_rotation;

test/unit-3d-test.cpp

@@ -251,6 +251,79 @@ BOOST_AUTO_TEST_CASE(dh_q_1)
 
 // --------------------------------------------------------------------
 
+BOOST_AUTO_TEST_CASE(m2q_1, *utf::tolerance(0.001f))
+{
+	cif::point pts[3] = {
+		{ 1, 2, 1 },
+		{ -2, -1, 1 },
+		{ 0, 0, 0 }
+	};
+
+	auto v1 = pts[0] - pts[2];
+	auto v2 = pts[1] - pts[2];
+
+	float a1 = std::acos(dot_product(v1, v2) / (v1.length() * v2.length()));
+
+	auto s = std::sin(a1 / 2);
+	auto c = std::cos(a1 / 2);
+
+	cif::quaternion q{ c, 0, 0, s };
+
+	auto pr = pts[0];
+	pr.rotate(q);
+	BOOST_TEST(pr.m_x == pts[1].m_x);
+	BOOST_TEST(pr.m_y == pts[1].m_y);
+	BOOST_TEST(pr.m_z == pts[1].m_z);
+
+	cif::matrix3x3<float> rot;
+	rot(0, 0) = 0;
+	rot(0, 1) = -1;
+	rot(0, 2) = 0;
+	rot(1, 0) = 1;
+	rot(1, 1) = -1;
+	rot(1, 2) = 0;
+	rot(2, 0) = 0;
+	rot(2, 1) = 0;
+	rot(2, 2) = 1;
+
+	pr = rot * pts[0];
+	BOOST_TEST(pr.m_x == pts[1].m_x);
+	BOOST_TEST(pr.m_y == pts[1].m_y);
+	BOOST_TEST(pr.m_z == pts[1].m_z);
+
+	auto &&[ev, em] = eigen(rot, false);
+
+	for (size_t i = 0; i < 3; ++i)
+	{
+		if (ev[i] != 1)
+			continue;
+		
+		auto tr = rot(0, 0) + rot(1, 1) + rot(2, 2);
+		auto a2 = std::acos((tr - 1) / 2.0f);
+		BOOST_TEST(2 * std::cos(a2) + 1 == tr);
+
+		auto s2 = std::sin(a2 / 2);
+		auto c2 = std::cos(a2 / 2);
+
+		auto q2 = normalize(cif::quaternion{
+			static_cast<float>(c2),
+			static_cast<float>(s2 * em(0, i)),
+			static_cast<float>(s2 * em(1, i)),
+			static_cast<float>(s2 * em(2, i)) });
+
+		pr = pts[0];
+		pr.rotate(q2);
+		BOOST_TEST(pr.m_x == pts[1].m_x);
+		BOOST_TEST(pr.m_y == pts[1].m_y);
+		BOOST_TEST(pr.m_z == pts[1].m_z);
+
+		break;
+	}
+
+}
+
+// --------------------------------------------------------------------
+
 BOOST_AUTO_TEST_CASE(symm_1)
 {
 	cif::cell c(10, 10, 10);