diff --git a/src/Math/Cylinder.h b/src/Math/Cylinder.h new file mode 100644 index 0000000..334c31a --- /dev/null +++ b/src/Math/Cylinder.h @@ -0,0 +1,141 @@ +/*////////////////////////////////////////////////////////////////////////////// +// CMT Cosmic Muon Tomography project ////////////////////////////////////////// +//////////////////////////////////////////////////////////////////////////////// + + Copyright (c) 2014, Universita' degli Studi di Padova, INFN sez. di Padova + All rights reserved + + Authors: Andrea Rigoni Garola < andrea.rigoni@pd.infn.it > + + ------------------------------------------------------------------ + This library 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 3.0 of the License, or (at your option) any later version. + + This library 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 this library. + +//////////////////////////////////////////////////////////////////////////////*/ + +#ifndef U_CYLINDER_H +#define U_CYLINDER_H + +#include "Geometry.h" +#include "Core/Object.h" +#include "Math/Dense.h" +#include "Math/Transform.h" + +namespace uLib { + +/** + * @brief Represents a cylindrical volume centered in the base circle. + * + * Cylinder inherits from AffineTransform, which defines its parent + * coordinate system. It contains an internal local transformation (m_LocalT) + * that defines the cylinder's actual volume (radius and height) + * relative to the emitter's origin (base circle center). + */ +class Cylinder : public AffineTransform, public Object { + typedef AffineTransform BaseClass; + +public: + /** + * @brief Default constructor. + * Initializes with radius 1 and height 1. + */ + Cylinder() : m_LocalT(this), m_Radius(1.0), m_Height(1.0) { + UpdateLocalMatrix(); + } + + /** + * @brief Constructor with radius and height. + */ + Cylinder(float radius, float height) : m_LocalT(this), m_Radius(radius), m_Height(height) { + UpdateLocalMatrix(); + } + + /** + * @brief Copy constructor. + */ + Cylinder(const Cylinder ©) + : m_LocalT(this), AffineTransform(copy) { + this->SetRadius(copy.GetRadius()); + this->SetHeight(copy.GetHeight()); + } + + /** Sets the radius of the cylinder */ + inline void SetRadius(float r) { + m_Radius = r; + UpdateLocalMatrix(); + } + + /** Gets the radius of the cylinder */ + inline float GetRadius() const { return m_Radius; } + + /** Sets the height of the cylinder */ + inline void SetHeight(float h) { + m_Height = h; + UpdateLocalMatrix(); + } + + /** Gets the height of the cylinder */ + inline float GetHeight() const { return m_Height; } + + /** + * @brief Returns the world transformation matrix of the cylinder's volume. + */ + Matrix4f GetWorldMatrix() const { return m_LocalT.GetWorldMatrix(); } + + /** + * @brief Returns the local transformation matrix of the cylinder's volume. + */ + Matrix4f GetLocalMatrix() const { return m_LocalT.GetMatrix(); } + + /** + * @brief Transforms local cylindrical coordinates to world space. + * @param r Local radius (absolute). + * @param theta Local angle in radians. + * @param z Local height (absolute, relative to base circle). + * @return Transformed point in world space. + */ + inline Vector4f GetWorldPoint(float r, float theta, float z) const { + return BaseClass::GetWorldMatrix() * Vector4f(r * std::cos(theta), r * std::sin(theta), z, 1.0f); + } + + /** + * @brief Transforms a world point to cylindrical local space. + * @return Vector3f(r, theta, z) + */ + inline Vector3f GetCylindricalLocal(const Vector4f &world_v) const { + Vector4f local_v = BaseClass::GetWorldMatrix().inverse() * world_v; + float r = std::sqrt(local_v.x() * local_v.x() + local_v.y() * local_v.y()); + float theta = std::atan2(local_v.y(), local_v.x()); + return Vector3f(r, theta, local_v.z()); + } + +signals: + /** Signal emitted when the cylinder geometry or transform is updated */ + virtual void Updated() override { ULIB_SIGNAL_EMIT(Cylinder::Updated); } + +private: + /** Recalculates the internal local matrix based on radius and height */ + void UpdateLocalMatrix() { + m_LocalT = AffineTransform(this); // BaseClass is parent + m_LocalT.Scale(Vector3f(m_Radius, m_Radius, m_Height)); + this->Updated(); + } + + float m_Radius; + float m_Height; + AffineTransform m_LocalT; +}; + +} // namespace uLib + +#endif // U_CYLINDER_H diff --git a/src/Math/testing/CMakeLists.txt b/src/Math/testing/CMakeLists.txt index 9fc8aa2..50f2197 100644 --- a/src/Math/testing/CMakeLists.txt +++ b/src/Math/testing/CMakeLists.txt @@ -3,6 +3,7 @@ set(TESTS MathVectorTest GeometryTest ContainerBoxTest + CylinderTest VoxImageTest VoxRaytracerTest VoxRaytracerTestExtended diff --git a/src/Math/testing/CylinderTest.cpp b/src/Math/testing/CylinderTest.cpp new file mode 100644 index 0000000..9f8575f --- /dev/null +++ b/src/Math/testing/CylinderTest.cpp @@ -0,0 +1,120 @@ +/*////////////////////////////////////////////////////////////////////////////// +// CMT Cosmic Muon Tomography project ////////////////////////////////////////// +//////////////////////////////////////////////////////////////////////////////// + + Copyright (c) 2014, Universita' degli Studi di Padova, INFN sez. di Padova + All rights reserved + + Authors: Andrea Rigoni Garola < andrea.rigoni@pd.infn.it > + + ------------------------------------------------------------------ + This library 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 3.0 of the License, or (at your option) any later version. + + This library 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 this library. + +//////////////////////////////////////////////////////////////////////////////*/ + +#include "testing-prototype.h" + +#include "Math/Dense.h" +#include "Math/Cylinder.h" + +#include +#include + +using namespace uLib; + +/** + * @brief Utility function to check if a 4D vector (spatial part) is zero within a threshold. + * Returns 0 if zero, 1 otherwise. + */ +int Vector4f0(Vector4f c) +{ + c(3) = 0; + if ( fabs(c(0)) < 0.001 && fabs(c(1)) < 0.001 && fabs(c(2)) < 0.001 ) + return 0; + else + return 1; +} + +int main() +{ + BEGIN_TESTING(Math Cylinder); + + // Test 1: Basic identity transformation and cylinder parameters + { + Cylinder cyl(2.0, 10.0); + std::cout << "Cyl R=" << cyl.GetRadius() << " H=" << cyl.GetHeight() << std::endl; + std::cout << "Cyl World Matrix:\n" << cyl.GetWorldMatrix() << std::endl; + + TEST0( std::abs(cyl.GetRadius() - 2.0) > 0.001 ); + TEST0( std::abs(cyl.GetHeight() - 10.0) > 0.001 ); + + // Point on the base circle center (Origin) + Vector4f p0 = cyl.GetWorldPoint(0, 0, 0); + std::cout << "p0: " << p0.transpose() << std::endl; + TEST0( Vector4f0(p0 - Vector4f(0, 0, 0, 1)) ); + + // Point on the top circle center (0, 0, Height) + Vector4f p1 = cyl.GetWorldPoint(0, 0, 10.0); + std::cout << "p1: " << p1.transpose() << std::endl; + TEST0( Vector4f0(p1 - Vector4f(0, 0, 10.0, 1)) ); + + // Point on the edge of the base circle at theta=0 (Radius, 0, 0) + Vector4f p2 = cyl.GetWorldPoint(2.0, 0, 0); + std::cout << "p2: " << p2.transpose() << std::endl; + TEST0( Vector4f0(p2 - Vector4f(2.0, 0, 0, 1)) ); + + // Point at 90 degrees on the side wall at middle height + Vector4f p3 = cyl.GetWorldPoint(2.0, M_PI/2.0, 5.0); + std::cout << "p3: " << p3.transpose() << std::endl; + TEST0( Vector4f0(p3 - Vector4f(0, 2.0, 5.0, 1)) ); + } + + // Test 2: Translation + { + Cylinder cyl(1.0, 2.0); + cyl.SetPosition(Vector3f(10, 20, 30)); + + // Local base origin (0, 0, 0) -> World (10, 20, 30) + Vector4f p0 = cyl.GetWorldPoint(0, 0, 0); + TEST0( Vector4f0(p0 - Vector4f(10, 20, 30, 1)) ); + + // Local top edge (1, 0, 2) -> World (11, 20, 32) + Vector4f p1 = cyl.GetWorldPoint(1, 0, 2); + TEST0( Vector4f0(p1 - Vector4f(11, 20, 32, 1)) ); + } + + // Test 3: Rotation and complex mapping + { + Cylinder cyl(5.0, 20.0); + cyl.SetPosition(Vector3f(1.0, 2.0, 3.0)); + // Rotate 90 degrees around X: Local Y becomes World Z, Local Z becomes World -Y + cyl.Rotate(M_PI/2.0, Vector3f(1, 0, 0)); + + // Let's take a local point: r=5, theta=pi/2, z=10 -> (0, 5, 10) in local cartesian + // Transformed: + // Position: (1, 2, 3) + // Point relative to position: (0, -10, 5) [since Z local -> -Y world, Y local -> Z world] + // Final World: (1, 2-10, 3+5) = (1, -8, 8) + Vector4f world_p = cyl.GetWorldPoint(5.0, M_PI/2.0, 10.0); + TEST0( Vector4f0(world_p - Vector4f(1.0, -8.0, 8.0, 1)) ); + + // Test inverse mapping + Vector3f cyl_local = cyl.GetCylindricalLocal(world_p); + TEST0( std::abs(cyl_local.x() - 5.0) > 0.001 ); + TEST0( std::abs(cyl_local.y() - M_PI/2.0) > 0.001 ); + TEST0( std::abs(cyl_local.z() - 10.0) > 0.001 ); + } + + END_TESTING; +}