add cylinder

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 &copy)
+        : 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

src/Math/testing/CMakeLists.txt

@@ -3,6 +3,7 @@ set(TESTS
         MathVectorTest
         GeometryTest
         ContainerBoxTest
+        CylinderTest
         VoxImageTest
         VoxRaytracerTest
         VoxRaytracerTestExtended

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 <cmath>
+#include <iostream>
+
+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;
+}