correct implementation of alignpoints

CMakeLists.txt

@@ -256,7 +256,6 @@ set(project_headers
 	${PROJECT_SOURCE_DIR}/include/cif++/CifUtils.hpp
 	${PROJECT_SOURCE_DIR}/include/cif++/CifValidator.hpp
 	${PROJECT_SOURCE_DIR}/include/cif++/Compound.hpp
-	${PROJECT_SOURCE_DIR}/include/cif++/Matrix.hpp
 	${PROJECT_SOURCE_DIR}/include/cif++/PDB2Cif.hpp
 	${PROJECT_SOURCE_DIR}/include/cif++/PDB2CifRemark3.hpp
 	${PROJECT_SOURCE_DIR}/include/cif++/Point.hpp

include/cif++/Point.hpp

@@ -380,7 +380,7 @@ PointF<F> Nudge(PointF<F> p, F offset);
 Quaternion Normalize(Quaternion q);
 
 std::tuple<double, Point> QuaternionToAngleAxis(Quaternion q);
-Point Centroid(std::vector<Point> &Points);
+Point Centroid(const std::vector<Point> &Points);
 Point CenterPoints(std::vector<Point> &Points);
 
 /// \brief Returns how the two sets of points \a a and \b b can be aligned

test/unit-test.cpp

@@ -1699,4 +1699,69 @@ BOOST_AUTO_TEST_CASE(reading_file_1)
 
 	cif::File file;
 	BOOST_CHECK_THROW(file.load(is), std::runtime_error);
-}
\ No newline at end of file
+}
+
+// 3d tests
+
+using namespace mmcif;
+
+BOOST_AUTO_TEST_CASE(t1)
+{
+	// std::random_device rnd;
+	// std::mt19937 gen(rnd());
+	// std::uniform_real_distribution<float> dis(0, 1);
+
+	// Quaternion q{ dis(gen), dis(gen), dis(gen), dis(gen) };
+	// q = Normalize(q);
+
+	// Quaternion q{ 0.1, 0.2, 0.3, 0.4 };
+	Quaternion q{ 0.5, 0.5, 0.5, 0.5 };
+	q = Normalize(q);
+
+	// std::cout << "q: " << q << std::endl;
+
+
+	std::vector<Point> p1{
+		{ 16.979, 13.301, 44.555 },
+		{ 18.150, 13.525, 43.680 },
+		{ 18.656, 14.966, 43.784 },
+		{ 17.890, 15.889, 44.078 },
+		{ 17.678, 13.270, 42.255 },
+		{ 16.248, 13.734, 42.347 },
+		{ 15.762, 13.216, 43.724 }
+	};
+
+	auto p2 = p1;
+
+	Point c1 = CenterPoints(p1);
+
+	std::cout << c1 << std::endl;
+
+	for (auto &p : p2)
+		p.rotate(q);
+	
+
+	Point c2 = CenterPoints(p2);
+
+	// std::cout << c2 << std::endl;
+
+	// std::cout << "rmsd: " << RMSd(p1, p2) << std::endl;
+
+	auto q2 = AlignPoints(p1, p2);
+
+	// std::cout << "q2: " << q2 << std::endl;
+
+
+	const auto &&[angle, axis] = QuaternionToAngleAxis(q2);
+
+	// std::cout << "rotated " << angle << " " << std::endl;
+
+	for (auto &p : p1)
+		p.rotate(q2);
+
+	float rmsd = RMSd(p1, p2);
+
+	BOOST_TEST(rmsd < 1e-5);
+
+	// std::cout << "rmsd: " << RMSd(p1, p2) << std::endl;
+}