add triangle mesh affine transform

src/Core/Object.cpp

@@ -145,6 +145,8 @@ GenericMFPtr *Object::findSlotImpl(const char *name) const {
   return NULL;
 }
 
+void Object::Updated() { ULIB_SIGNAL_EMIT(Object::Updated); }
+
 // std::ostream &
 // operator << (std::ostream &os, uLib::Object &ob)
 // {

src/Core/Object.h

@@ -97,6 +97,9 @@ public:
   ////////////////////////////////////////////////////////////////////////////
   // SIGNALS //
 
+  signals:
+  virtual void Updated();
+
   // Qt4 style connector //
   static bool connect(const Object *ob1, const char *signal_name,
                       const Object *receiver, const char *slot_name) {
@@ -121,6 +124,25 @@ public:
     return true;
   }
 
+  // Lambda/Function object connector //
+  template <typename Func1, typename SlotT>
+  static bool connect(typename FunctionPointer<Func1>::Object *sender,
+                      Func1 sigf, SlotT slof) {
+    SignalBase *sigb = sender->findOrAddSignal(sigf);
+    typedef typename FunctionPointer<Func1>::SignalSignature SigSignature;
+    typedef typename Signal<SigSignature>::type SigT;
+    reinterpret_cast<SigT *>(sigb)->connect(slof);
+    return true;
+  }
+
+  template <typename Func1, typename Func2>
+  static bool
+  disconnect(typename FunctionPointer<Func1>::Object *sender, Func1 sigf,
+             typename FunctionPointer<Func2>::Object *receiver, Func2 slof) {
+    // TODO: implement actual disconnect in Signal.h //
+    return true;
+  }
+
   template <typename FuncT>
   static inline bool connect(SignalBase *sigb, FuncT slof, Object *receiver) {
     ConnectSignal<typename FunctionPointer<FuncT>::SignalSignature>(sigb, slof,

src/HEP/Detectors/CMakeLists.txt

@@ -11,19 +11,31 @@ set(HEADERS
     MuonScatter.h
 )
 
+set(SOURCES
+    DetectorChamber.cpp
+)
+
+set(LIBRARIES
+       ${PACKAGE_LIBPREFIX}Core
+       ${PACKAGE_LIBPREFIX}Math
+)
+
 set(libname ${PACKAGE_LIBPREFIX}Detectors)
 set(ULIB_SHARED_LIBRARIES ${ULIB_SHARED_LIBRARIES} ${libname} PARENT_SCOPE)
 set(ULIB_SELECTED_MODULES ${ULIB_SELECTED_MODULES} Detectors PARENT_SCOPE)
 
-## Headers-only INTERFACE library
-add_library(${libname} INTERFACE)
-target_include_directories(${libname} INTERFACE
-    $<BUILD_INTERFACE:${SRC_DIR}>
-    $<INSTALL_INTERFACE:${INSTALL_INC_DIR}>
-)
+## SHARED library
+add_library(${libname} SHARED ${SOURCES})
+set_target_properties(${libname} PROPERTIES
+                      VERSION ${PROJECT_VERSION}
+                      SOVERSION ${PROJECT_SOVERSION}
+                      CXX_STANDARD 17)
+target_link_libraries(${libname} PRIVATE ${LIBRARIES})
 
 install(TARGETS ${libname}
-        EXPORT "uLibTargets")
+        EXPORT "uLibTargets"
+        RUNTIME DESTINATION ${INSTALL_BIN_DIR} COMPONENT bin
+        LIBRARY DESTINATION ${INSTALL_LIB_DIR} COMPONENT lib)
 
 install(FILES ${HEADERS}
         DESTINATION ${INSTALL_INC_DIR}/HEP/Detectors)

src/HEP/Detectors/DetectorChamber.cpp

@@ -0,0 +1,53 @@
+#include "HEP/Detectors/DetectorChamber.h"
+#include <cmath>
+
+namespace uLib {
+
+MuonEvent DetectorChamber::ProjectMuonEvent(const MuonEvent &muon) const {
+  MuonEvent projectedMuon = muon;
+
+  // Transform the local projection plane to world coordinates
+  HLine3f worldPlane = this->GetWorldProjectionPlane();
+  HPoint3f P = worldPlane.origin;
+  HVector3f N = worldPlane.direction;
+
+  HPoint3f X_in = muon.LineIn().origin;
+  HPoint3f X_out = muon.LineOut().origin;
+
+  // Calculate squared distances to the plane normal point for comparison
+  // Actually, we should probably follow the user's description literally:
+  // "closest ... with the projection plane ( so the colsest direction point with the point of the normal defining the plane )"
+  // This could mean point-to-plane or point-to-point. 
+  // Given "closest with the projection plane", point-to-plane is more natural.
+  // However, "closest direction point with the point of the normal" strongly suggests point-to-point distance.
+  // Let's use distance to the plane for the first part and keep the logic consistent.
+  
+  float dist_in = std::abs((X_in - P).dot(N));
+  float dist_out = std::abs((X_out - P).dot(N));
+
+  const HLine3f &chosenLine = (dist_in <= dist_out) ? muon.LineIn() : muon.LineOut();
+  HPoint3f X_chosen = chosenLine.origin;
+
+  // Project X_chosen into the plane defined by P and normal N
+  // X_proj = X_chosen - ((X_chosen - P) . N / (N . N)) * N
+  float dot = (X_chosen - P).dot(N);
+  float n_sq = N.dot(N);
+
+  HPoint3f X_proj = X_chosen;
+  if (n_sq > 0) {
+    X_proj = X_chosen - (dot / n_sq) * N;
+  }
+
+  // Define the projected line with projected origin and original direction
+  HLine3f projectedLine;
+  projectedLine.origin = X_proj;
+  projectedLine.direction = chosenLine.direction;
+
+  // Set both input and output lines of the projected muon to the same projected line
+  projectedMuon.LineIn() = projectedLine;
+  projectedMuon.LineOut() = projectedLine;
+
+  return projectedMuon;
+}
+
+} // namespace uLib

src/HEP/Detectors/DetectorChamber.h

@@ -31,14 +31,48 @@
 #include "Core/Types.h"
 #include "Math/ContainerBox.h"
 
+#include "HEP/Detectors/Hit.h"
+#include "HEP/Detectors/HitMC.h"
+#include "HEP/Detectors/MuonEvent.h"
+
 namespace uLib {
 
 
 class DetectorChamber : public ContainerBox {
 
+  typedef ContainerBox BaseClass;
+  
+
 public:
 
+  DetectorChamber() : BaseClass() {
+    m_ProjectionPlane.origin = HPoint3f(0, 0, 0);
+    m_ProjectionPlane.direction = HVector3f(0, 0, 1);
+  }
+
+  DetectorChamber(const Vector3f &size) : BaseClass(size) {
+    m_ProjectionPlane.origin = HPoint3f(0, 0, 0);
+    m_ProjectionPlane.direction = HVector3f(0, 0, 1);
+  }
+
+  // set the plane where muons hit is projected 
+  // coordinates are local to the container box
+  void SetProjectionPlane(const HLine3f &normal) { m_ProjectionPlane = normal; }
+
+  const HLine3f &GetProjectionPlane() const { return m_ProjectionPlane; }
+
+  HLine3f GetWorldProjectionPlane() const {
+    HLine3f worldPlane;
+    Matrix4f M = this->GetWorldMatrix();
+    worldPlane.origin = M * m_ProjectionPlane.origin;
+    worldPlane.direction = M * m_ProjectionPlane.direction;
+    return worldPlane;
+  }
+
+  MuonEvent ProjectMuonEvent(const MuonEvent &muon) const;
+
 private:
+  HLine3f m_ProjectionPlane;
 };
 
 

src/HEP/Detectors/MuonEvent.h

@@ -48,6 +48,10 @@ protected:
 
 class MuonEvent  : public MuonEventData {
 public:
+    using MuonEventData::LineIn;
+    using MuonEventData::LineOut;
+    using MuonEventData::GetMomentum;
+
     inline HLine3f & LineIn() { return this->m_LineIn; }
     inline HLine3f & LineOut() { return this->m_LineOut; }
     inline Scalarf & Momentum() { return this->m_Momentum; }

src/HEP/Detectors/testing/CMakeLists.txt

@@ -1,12 +1,13 @@
 # TESTS
 set( TESTS
-#            GDMLSolidTest
             HierarchicalEncodingTest
+            DetectorChamberTest
 )
 
 set(LIBRARIES
        ${PACKAGE_LIBPREFIX}Core
        ${PACKAGE_LIBPREFIX}Math
+       ${PACKAGE_LIBPREFIX}Detectors
        Boost::serialization
        Boost::program_options
        Eigen3::Eigen

src/HEP/Detectors/testing/DetectorChamberTest.cpp

@@ -0,0 +1,139 @@
+/*//////////////////////////////////////////////////////////////////////////////
+// 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 "HEP/Detectors/DetectorChamber.h"
+#include "testing-prototype.h"
+#include <iostream>
+
+using namespace uLib;
+
+int main() {
+    BEGIN_TESTING(DetectorChamber Projection);
+
+    DetectorChamber chamber;
+
+    // Define a horizontal plane at z = 100
+    HLine3f plane;
+    plane.origin = HPoint3f(0, 0, 100);
+    plane.direction = HVector3f(0, 0, 1); // Normal to the plane
+    chamber.SetProjectionPlane(plane);
+
+    // Create a muon with two segments
+    MuonEvent muon;
+    
+    // Segment 1 (muon.LineIn()): origin is at (10, 20, 50), direction along Z
+    muon.LineIn().origin = HPoint3f(10, 20, 50);
+    muon.LineIn().direction = HVector3f(0, 0, 1);
+    
+    // Segment 2 (muon.LineOut()): origin is at (10, 20, 200), direction along Z
+    muon.LineOut().origin = HPoint3f(10, 20, 200);
+    muon.LineOut().direction = HVector3f(0, 0, 1);
+    
+    // distance_in = |50 - 100| = 50
+    // distance_out = |200 - 100| = 100
+    // LineIn is closer to the plane (z=100)
+    
+    MuonEvent projected = chamber.ProjectMuonEvent(muon);
+
+    // Expected:
+    // chosenLine = LineIn (it is closer)
+    // X_chosen = (10, 20, 50)
+    // X_proj = (10, 20, 50) - (( (10, 20, 50) - (0, 0, 100) ) . (0, 0, 1)) * (0, 0, 1)
+    // X_proj = (10, 20, 50) - (-50) * (0, 0, 1) = (10, 20, 100)
+    
+    HPoint3f expectedPos(10, 20, 100);
+    
+    std::cout << "Test Case 1: LineIn is closer" << std::endl;
+    std::cout << "Projected Position: " << projected.LineIn().origin.transpose() << std::endl;
+    std::cout << "Expected Position:  " << expectedPos.transpose() << std::endl;
+
+    // Check if the projected position is correct
+    // norm() includes the 4th component, which for HVector3f (diff of points) should be 0.
+    bool posOk1 = (projected.LineIn().origin - expectedPos).norm() < 1e-5;
+    TEST1(posOk1);
+
+    // Check if LineIn and LineOut are the same
+    bool linesMatch1 = (projected.LineIn().origin == projected.LineOut().origin) &&
+                       (projected.LineIn().direction == projected.LineOut().direction);
+    TEST1(linesMatch1);
+    
+    // Test Case 2: LineOut is closer
+    muon.LineIn().origin = HPoint3f(30, 40, 0);      // dist = 100
+    muon.LineOut().origin = HPoint3f(30, 40, 110);   // dist = 10
+    
+    projected = chamber.ProjectMuonEvent(muon);
+    expectedPos = HPoint3f(30, 40, 100); // projection of (30,40,110) onto z=100
+    
+    std::cout << "\nTest Case 2: LineOut is closer" << std::endl;
+    std::cout << "Projected Position: " << projected.LineIn().origin.transpose() << std::endl;
+    std::cout << "Expected Position:  " << expectedPos.transpose() << std::endl;
+
+    bool posOk2 = (projected.LineIn().origin - expectedPos).norm() < 1e-5;
+    TEST1(posOk2);
+    
+    // Test Case 3: Oblique plane
+    // Plane through (0,0,0) with normal (1,1,1) (normalized)
+    plane.origin = HPoint3f(0, 0, 0);
+    plane.direction = HVector3f(1, 1, 1).normalized();
+    chamber.SetProjectionPlane(plane);
+    
+    muon.LineIn().origin = HPoint3f(1, 1, 1); // dist = (1,1,1) . (1,1,1).norm()
+    muon.LineIn().direction = HVector3f(0, 0, 1);
+    muon.LineOut().origin = HPoint3f(10, 10, 10); // dist = 10 * sqrt(3)
+    
+    projected = chamber.ProjectMuonEvent(muon);
+    // X_chosen = (1,1,1)
+    // X_proj = (1,1,1) - ((1,1,1) . N) * N   where N = (1,1,1)/sqrt(3)
+    // X_proj = (1,1,1) - (sqrt(3)) * (1,1,1)/sqrt(3) = (1,1,1) - (1,1,1) = (0,0,0)
+    expectedPos = HPoint3f(0, 0, 0);
+    
+    std::cout << "\nTest Case 3: Oblique plane" << std::endl;
+    std::cout << "Projected Position: " << projected.LineIn().origin.transpose() << std::endl;
+    std::cout << "Expected Position:  " << expectedPos.transpose() << std::endl;
+    
+    bool posOk3 = (projected.LineIn().origin - expectedPos).norm() < 1e-5;
+    TEST1(posOk3);
+
+    // Test Case 4: Transformed DetectorChamber
+    DetectorChamber chamber2;
+    chamber2.SetPosition(Vector3f(0, 0, 100)); // Move chamber to z=100
+    // chamber2.GetProjectionPlane has default origin (0,0,0) and direction (0,0,1)
+    // In world coordinates, this plane is at z = 100 + 0 = 100.
+    
+    muon.LineIn().origin = HPoint3f(50, 60, 50); // dist to world plane (z=100) is 50
+    muon.LineOut().origin = HPoint3f(50, 60, 200); // dist to world plane (z=100) is 100
+    
+    projected = chamber2.ProjectMuonEvent(muon);
+    expectedPos = HPoint3f(50, 60, 100);
+    
+    std::cout << "\nTest Case 4: Transformed DetectorChamber (active world matrix)" << std::endl;
+    std::cout << "Projected Position: " << projected.LineIn().origin.transpose() << std::endl;
+    std::cout << "Expected Position:  " << expectedPos.transpose() << std::endl;
+    
+    bool posOk4 = (projected.LineIn().origin - expectedPos).norm() < 1e-5;
+    TEST1(posOk4);
+
+    END_TESTING;
+}

src/HEP/Detectors/testing/Makefile.am

@@ -1,16 +0,0 @@
-
-include $(top_srcdir)/Common.am
-
-#AM_DEFAULT_SOURCE_EXT = .cpp
-
-# if HAVE_CHECK
-TESTS = GDMLSolidTest
-
-# else
-# TEST =
-# endif
-
-LDADD = $(top_srcdir)/libmutom-${PACKAGE_VERSION}.la
-
-check_PROGRAMS = $(TESTS)
-

src/HEP/Geant/EmitterPrimary.cpp

@@ -71,14 +71,13 @@ void QuadMeshEmitterPrimary::CalculateAreas() {
   m_CumulativeAreas.clear();
   m_TotalArea = 0.0;
 
-  const auto &points = m_Mesh->Points();
   const auto &quads = m_Mesh->Quads();
 
   for (const auto &q : quads) {
-    uLib::Vector3f v0 = points[q(0)];
-    uLib::Vector3f v1 = points[q(1)];
-    uLib::Vector3f v2 = points[q(2)];
-    uLib::Vector3f v3 = points[q(3)];
+    uLib::Vector3f v0 = m_Mesh->GetPoint(q(0));
+    uLib::Vector3f v1 = m_Mesh->GetPoint(q(1));
+    uLib::Vector3f v2 = m_Mesh->GetPoint(q(2));
+    uLib::Vector3f v3 = m_Mesh->GetPoint(q(3));
 
     double a1 = 0.5 * (v1 - v0).cross(v2 - v0).norm();
     double a2 = 0.5 * (v2 - v0).cross(v3 - v0).norm();
@@ -96,11 +95,10 @@ void QuadMeshEmitterPrimary::GeneratePrimaries(G4Event *anEvent) {
   int quadIdx = std::distance(m_CumulativeAreas.begin(), it);
 
   const auto &q = m_Mesh->Quads()[quadIdx];
-  const auto &points = m_Mesh->Points();
-  uLib::Vector3f v0 = points[q(0)];
-  uLib::Vector3f v1 = points[q(1)];
-  uLib::Vector3f v2 = points[q(2)];
-  uLib::Vector3f v3 = points[q(3)];
+  uLib::Vector3f v0 = m_Mesh->GetPoint(q(0));
+  uLib::Vector3f v1 = m_Mesh->GetPoint(q(1));
+  uLib::Vector3f v2 = m_Mesh->GetPoint(q(2));
+  uLib::Vector3f v3 = m_Mesh->GetPoint(q(3));
 
   // 2. Choose a point on the quad
   double a1 = 0.5 * (v1 - v0).cross(v2 - v0).norm();

src/Math/ContainerBox.h

@@ -180,7 +180,7 @@ public:
 signals:
 
   // signal to emit when the box is updated //
-  void Updated() { ULIB_SIGNAL_EMIT(ContainerBox::Updated); }
+  virtual void Updated() override { ULIB_SIGNAL_EMIT(ContainerBox::Updated); }
 
 private:
   AffineTransform m_LocalT;

src/Math/QuadMesh.cpp

@@ -35,10 +35,10 @@ void QuadMesh::PrintSelf(std::ostream &o)
     o << "    #Quads  : " << m_Quads.size() << "\n";
     for(int i=0; i < m_Quads.size(); ++i ) {
         o << "     - quad[" << i << "]" <<
-             " " << m_Quads[i](0) << "->(" << m_Points[m_Quads[i](0)].transpose() << ") " <<
-             " " << m_Quads[i](1) << "->(" << m_Points[m_Quads[i](1)].transpose() << ") " <<
-             " " << m_Quads[i](2) << "->(" << m_Points[m_Quads[i](2)].transpose() << ") " <<
-             " " << m_Quads[i](3) << "->(" << m_Points[m_Quads[i](3)].transpose() << ") " <<
+             " " << m_Quads[i](0) << "->(" << GetPoint(m_Quads[i](0)).transpose() << ") " <<
+             " " << m_Quads[i](1) << "->(" << GetPoint(m_Quads[i](1)).transpose() << ") " <<
+             " " << m_Quads[i](2) << "->(" << GetPoint(m_Quads[i](2)).transpose() << ") " <<
+             " " << m_Quads[i](3) << "->(" << GetPoint(m_Quads[i](3)).transpose() << ") " <<
              " \n";
     }
     o << " // ------------------------- // \n";
@@ -46,7 +46,16 @@ void QuadMesh::PrintSelf(std::ostream &o)
 
 void QuadMesh::AddPoint(const Vector3f &pt)
 {
-    this->m_Points.push_back(pt);
+    Vector4f p(pt.x(), pt.y(), pt.z(), 1.0f);
+    Vector4f localP = this->GetWorldMatrix().inverse() * p;
+    this->m_Points.push_back(localP.head<3>());
+}
+
+Vector3f QuadMesh::GetPoint(const Id_t id) const
+{
+    Vector4f p(m_Points.at(id).x(), m_Points.at(id).y(), m_Points.at(id).z(), 1.0f);
+    Vector4f worldP = this->GetWorldMatrix() * p;
+    return worldP.head<3>();
 }
 
 void QuadMesh::AddQuad(const Id_t *id)
@@ -63,9 +72,9 @@ void QuadMesh::AddQuad(const Vector4i &id)
 Vector3f QuadMesh::GetNormal(const Id_t id) const
 {
     const Vector4i &quad = m_Quads.at(id);
-    const Vector3f &v0 = m_Points.at(quad(0));
-    const Vector3f &v1 = m_Points.at(quad(1));
-    const Vector3f &v3 = m_Points.at(quad(3));
+    const Vector3f v0 = this->GetPoint(quad(0));
+    const Vector3f v1 = this->GetPoint(quad(1));
+    const Vector3f v3 = this->GetPoint(quad(3));
     
     Vector3f edge1 = v1 - v0;
     Vector3f edge2 = v3 - v0;

src/Math/QuadMesh.h

@@ -29,25 +29,33 @@
 #include <vector>
 
 #include "Math/Dense.h"
+#include "Core/Object.h"
+#include "Math/Transform.h"
 
 namespace uLib {
 
-class QuadMesh
+class QuadMesh : public AffineTransform, public Object
 {
 public:
     void PrintSelf(std::ostream &o);
 
+    /** @brief Adds a point in global coordinates. Stored in local coordinates. */
     void AddPoint(const Vector3f &pt);
 
     void AddQuad(const Id_t *id);
     void AddQuad(const Vector4i &id);
 
+    /** @brief Returns point in global coordinates. */
+    Vector3f GetPoint(const Id_t id) const;
+
     inline std::vector<Vector3f> & Points() { return this->m_Points; }
     inline std::vector<Vector4i> & Quads() { return this->m_Quads; }
 
     const Vector4i & GetQuad(const Id_t id) const { return m_Quads.at(id); }
     Vector3f GetNormal(const Id_t id) const;
 
+    virtual void Updated() override { ULIB_SIGNAL_EMIT(QuadMesh::Updated); }
+
 private:
     std::vector<Vector3f> m_Points;
     std::vector<Vector4i> m_Quads;

src/Math/Transform.h

@@ -101,7 +101,7 @@ public:
 
     inline Matrix3f GetRotation() const { return this->m_T.rotation(); }
 
-    inline void Translate(const Vector3f v) { this->m_T.translate(v); }
+    inline void Translate(const Vector3f v) { this->m_T.pretranslate(v); }
 
     inline void Scale(const Vector3f v) { this->m_T.scale(v); }
 

src/Math/TriangleMesh.cpp

@@ -34,14 +34,14 @@ void TriangleMesh::PrintSelf(std::ostream &o)
     o << " // ------- TRIANGLE MESH ------- // \n" ;
     o << "    #Points : " << m_Points.size() << "\n";
     o << "    #Triang : " << m_Triangles.size() << "\n";
-    for(int i=0; i < m_Triangles.size(); ++i ) {
+    for(int i=0; i < (int)m_Triangles.size(); ++i ) {
         o << "     - triangle[" << i << "]" <<
              " " << m_Triangles[i](0) <<
-             "->(" << m_Points[m_Triangles[i](0)].transpose() << ") " <<
+             "->(" << GetPoint(m_Triangles[i](0)).transpose() << ") " <<
              " " << m_Triangles[i](1) <<
-             "->(" << m_Points[m_Triangles[i](1)].transpose() << ") " <<
+             "->(" << GetPoint(m_Triangles[i](1)).transpose() << ") " <<
              " " << m_Triangles[i](2) <<
-             "->(" << m_Points[m_Triangles[i](2)].transpose() << ") " <<
+             "->(" << GetPoint(m_Triangles[i](2)).transpose() << ") " <<
              " \n";
     }
     o << " // ----------------------------- // \n";
@@ -49,7 +49,16 @@ void TriangleMesh::PrintSelf(std::ostream &o)
 
 void TriangleMesh::AddPoint(const Vector3f &pt)
 {
-    this->m_Points.push_back(pt);
+    Vector4f p(pt.x(), pt.y(), pt.z(), 1.0f);
+    Vector4f localP = this->GetWorldMatrix().inverse() * p;
+    this->m_Points.push_back(localP.head<3>());
+}
+
+Vector3f TriangleMesh::GetPoint(const Id_t id) const
+{
+    Vector4f p(m_Points.at(id).x(), m_Points.at(id).y(), m_Points.at(id).z(), 1.0f);
+    Vector4f worldP = this->GetWorldMatrix() * p;
+    return worldP.head<3>();
 }
 
 void TriangleMesh::AddTriangle(const Id_t *id)
@@ -66,9 +75,9 @@ void TriangleMesh::AddTriangle(const Vector3i &id)
 Vector3f TriangleMesh::GetNormal(const Id_t id) const
 {
     const Vector3i &trg = m_Triangles.at(id);
-    const Vector3f &v0 = m_Points.at(trg(0));
-    const Vector3f &v1 = m_Points.at(trg(1));
-    const Vector3f &v2 = m_Points.at(trg(2));
+    const Vector3f v0 = this->GetPoint(trg(0));
+    const Vector3f v1 = this->GetPoint(trg(1));
+    const Vector3f v2 = this->GetPoint(trg(2));
     
     Vector3f edge1 = v1 - v0;
     Vector3f edge2 = v2 - v0;

src/Math/TriangleMesh.h

@@ -32,24 +32,33 @@
 
 #include "Math/Dense.h"
 
+#include "Core/Object.h"
+#include "Math/Transform.h"
+
 namespace uLib {
 
-class TriangleMesh
+class TriangleMesh : public AffineTransform, public Object
 {
 public:
     void PrintSelf(std::ostream &o);
 
+    /** @brief Adds a point in global coordinates. Stored in local coordinates. */
     void AddPoint(const Vector3f &pt);
 
     void AddTriangle(const Id_t *id);
     void AddTriangle(const Vector3i &id);
 
+    /** @brief Returns point in global coordinates. */
+    Vector3f GetPoint(const Id_t id) const;
+
     inline std::vector<Vector3f> & Points() { return this->m_Points; }
     inline std::vector<Vector3i> & Triangles() { return this->m_Triangles; }
 
     const Vector3i & GetTriangle(const Id_t id) const { return m_Triangles.at(id); }
     Vector3f GetNormal(const Id_t id) const;
 
+    virtual void Updated() override { ULIB_SIGNAL_EMIT(TriangleMesh::Updated); }
+
 private:
     std::vector<Vector3f> m_Points;
     std::vector<Vector3i> m_Triangles;

src/Math/testing/QuadMeshTest.cpp

@@ -47,5 +47,17 @@ int main() {
   ASSERT_EQ(mesh.Points().size(), 4);
   ASSERT_EQ(mesh.Quads().size(), 1);
 
+  // Test transformation
+  mesh.Translate(Vector3f(10, 20, 30));
+  Vector3f p0 = mesh.GetPoint(0);
+  TEST1( (p0 - Vector3f(10, 20, 30)).norm() < 1e-6 );
+  
+  // Test AddPoint during transformation
+  mesh.AddPoint(Vector3f(11, 21, 31)); // Should be stored as (1, 1, 1) locally
+  Id_t lastId = mesh.Points().size() - 1;
+  TEST1( (mesh.Points().at(lastId) - Vector3f(1, 1, 1)).norm() < 1e-5 );
+
+  mesh.PrintSelf(std::cout);
+
   END_TESTING;
 }

src/Math/testing/TriangleMeshTest.cpp

@@ -43,5 +43,20 @@ int main() {
 
   mesh.PrintSelf(std::cout);
 
+  TEST1(mesh.Points().size() == 3);
+  TEST1(mesh.Triangles().size() == 1);
+
+  // Test transformation
+  mesh.Translate(Vector3f(10, 20, 30));
+  Vector3f p0 = mesh.GetPoint(0);
+  TEST1( (p0 - Vector3f(10, 20, 30)).norm() < 1e-6 );
+  
+  // Test AddPoint during transformation
+  mesh.AddPoint(Vector3f(11, 21, 31)); // Should be stored as (1, 1, 1) locally
+  Id_t lastId = mesh.Points().size() - 1;
+  TEST1( (mesh.Points().at(lastId) - Vector3f(1, 1, 1)).norm() < 1e-5 );
+
+  mesh.PrintSelf(std::cout);
+
   END_TESTING;
 }

src/Vtk/HEP/Detectors/testing/vtkDetectorChamberTest.cpp

@@ -36,16 +36,12 @@ using namespace uLib;
 
 BOOST_AUTO_TEST_CASE(vtkDetectorChamberTest) {
   DetectorChamber d1, d2;
-  // d1.SetSize(Vector3f(1, 1, 1));
-  // d1.SetPosition(Vector3f(0, 0, 0));
-  d1.Scale(Vector3f(5_m, 10_m, 2_m));
+  d1.Scale(Vector3f(1_m, 2_m, 20_cm));
   d1.Translate(Vector3f(0, 0, 0));
 
-  // d2.SetSize(Vector3f(1, 1, 1));
-  // d2.SetPosition(Vector3f(0, 0, 0));
-  d2.Scale(Vector3f(5_m, 10_m, 2_m));
-  d2.Translate(Vector3f(0, 0, 10_m));
-
+  d2.Rotate(180_deg, Vector3f(0, 1, 0));
+  d2.Scale(Vector3f(1_m, 2_m, 20_cm));
+  d2.Translate(Vector3f(1_m, 0, 10_m));
 
   Vtk::vtkDetectorChamber v_d1(&d1);
   Vtk::vtkDetectorChamber v_d2(&d2);

src/Vtk/HEP/Detectors/vtkDetectorChamber.cxx

@@ -41,20 +41,89 @@
 #include "Vtk/HEP/Detectors/vtkDetectorChamber.h"
 #include <vtkBoxWidget.h>
 #include <vtkTransformPolyDataFilter.h>
+#include <vtkPlaneSource.h>
+#include <vtkProperty.h>
+#include <vtkNew.h>
 
 namespace uLib {
 namespace Vtk {
 
 vtkDetectorChamber::vtkDetectorChamber(DetectorChamber *content)
     : vtkContainerBox(content) {
+  m_PlaneSource = vtkPlaneSource::New();
+  
+  vtkNew<vtkPolyDataMapper> mapper;
+  mapper->SetInputConnection(m_PlaneSource->GetOutputPort());
+  
+  m_PlaneActor = vtkActor::New();
+  m_PlaneActor->SetMapper(mapper);
+  m_PlaneActor->GetProperty()->SetColor(0.2, 0.8, 0.2); // Light green
+  m_PlaneActor->GetProperty()->SetOpacity(0.3);
+  m_PlaneActor->GetProperty()->SetAmbient(1.0);
+  m_PlaneActor->GetProperty()->SetDiffuse(0.0);
+  
+  this->SetProp(m_PlaneActor);
+  
+  this->contentUpdate();
 }
 
 vtkDetectorChamber::~vtkDetectorChamber() {
+  m_PlaneSource->Delete();
+  m_PlaneActor->Delete();
 }
 
 DetectorChamber *vtkDetectorChamber::GetContent() {
   return static_cast<DetectorChamber *>(m_Content);
 }
 
+void vtkDetectorChamber::contentUpdate() {
+  this->BaseClass::contentUpdate();
+
+  if (!m_Content) return;
+  DetectorChamber *c = this->GetContent();
+  Vector3f size = c->GetSize();
+  HLine3f plane = c->GetProjectionPlane();
+
+  // Normalized local space of the chamber
+  float Lx = plane.origin.x() / size.x();
+  float Ly = plane.origin.y() / size.y();
+  float Lz = plane.origin.z() / size.z();
+  
+  float Dx = plane.direction.x() * size.x();
+  float Dy = plane.direction.y() * size.y();
+  float Dz = plane.direction.z() * size.z();
+
+  Vector3f normal(Dx, Dy, Dz);
+  normal.normalize();
+
+  // Find spans
+  Vector3f v1;
+  if (std::abs(normal.z()) < 0.9) v1 = Vector3f(0,0,1);
+  else v1 = Vector3f(1,0,0);
+  
+  Vector3f p1 = normal.cross(v1).normalized();
+  Vector3f p2 = normal.cross(p1).normalized();
+  
+  // Center the visual representation on the box extents (unit box [0,1]^3)
+  // instead of the plane origin (which might be a corner).
+  Vector3f boxCenter(0.5f, 0.5f, 0.5f);
+  Vector3f planePt(Lx, Ly, Lz);
+  Vector3f visualCenter = boxCenter - (boxCenter - planePt).dot(normal) * normal;
+  
+  float scale = 1.5; // Slightly larger than the diagonal of a face (1.41)
+  
+  m_PlaneSource->SetOrigin(visualCenter.x() - 0.5f*scale*p1.x() - 0.5f*scale*p2.x(),
+                           visualCenter.y() - 0.5f*scale*p1.y() - 0.5f*scale*p2.y(),
+                           visualCenter.z() - 0.5f*scale*p1.z() - 0.5f*scale*p2.z());
+                           
+  m_PlaneSource->SetPoint1(visualCenter.x() + 0.5f*scale*p1.x() - 0.5f*scale*p2.x(),
+                           visualCenter.y() + 0.5f*scale*p1.y() - 0.5f*scale*p2.y(),
+                           visualCenter.z() + 0.5f*scale*p1.z() - 0.5f*scale*p2.z());
+                           
+  m_PlaneSource->SetPoint2(visualCenter.x() - 0.5f*scale*p1.x() + 0.5f*scale*p2.x(),
+                           visualCenter.y() - 0.5f*scale*p1.y() + 0.5f*scale*p2.y(),
+                           visualCenter.z() - 0.5f*scale*p1.z() + 0.5f*scale*p2.z());
+}
+
 } // namespace Vtk
 } // namespace uLib

src/Vtk/HEP/Detectors/vtkDetectorChamber.h

@@ -39,6 +39,8 @@
 #include <vtkBoxWidget.h>
 #include <vtkTransformPolyDataFilter.h>
 
+class vtkPlaneSource;
+
 namespace uLib {
 namespace Vtk {
 
@@ -53,6 +55,12 @@ public:
   virtual ~vtkDetectorChamber();
 
   Content *GetContent();
+
+  virtual void contentUpdate() override;
+
+protected:
+  vtkActor *m_PlaneActor;
+  vtkPlaneSource *m_PlaneSource;
 };
 
 } // namespace Vtk

src/Vtk/Math/testing/vtkQuadMeshTest.cpp

@@ -43,6 +43,17 @@ BOOST_AUTO_TEST_CASE(vtkQuadMeshConstruction) {
   mesh.AddQuad(Vector4i(0, 1, 2, 3));
 
   Vtk::vtkQuadMesh v_mesh(mesh);
+
+  Object::connect(&mesh, &QuadMesh::Updated, [&mesh]() {
+    Vector3f points[4];
+    points[0] = mesh.GetPoint(0);
+    points[1] = mesh.GetPoint(1);
+    points[2] = mesh.GetPoint(2);
+    points[3] = mesh.GetPoint(3);
+    std::cout << "mesh updated: " << points[0] << " " << points[1] 
+              << " " << points[2] << " " << points[3] << std::endl;
+  });
+
   v_mesh.Update();
 
   if (std::getenv("CTEST_PROJECT_NAME") == nullptr) {

src/Vtk/Math/testing/vtkTriangleMeshTest.cpp

@@ -41,6 +41,16 @@ BOOST_AUTO_TEST_CASE(vtkTriangleMeshConstruction) {
   mesh.AddTriangle(Vector3i(0, 1, 2));
 
   Vtk::vtkTriangleMesh v_mesh(mesh);
+
+  Object::connect(&mesh, &TriangleMesh::Updated, [&mesh]() {
+    Vector3f points[3];
+    points[0] = mesh.GetPoint(0);
+    points[1] = mesh.GetPoint(1);
+    points[2] = mesh.GetPoint(2);
+    std::cout << "mesh updated: " << points[0].transpose() << " " << points[1].transpose() 
+              << " " << points[2].transpose() << std::endl;
+  });
+
   v_mesh.Update();
 
   if (std::getenv("CTEST_PROJECT_NAME") == nullptr) {

src/Vtk/Math/vtkQuadMesh.cpp

@@ -40,6 +40,9 @@
 #include <vtkPolyData.h>
 #include <vtkPolyDataMapper.h>
 
+#include <vtkMatrix4x4.h>
+#include <vtkNew.h>
+#include "Math/vtkDense.h"
 #include "Vtk/Math/vtkQuadMesh.h"
 #include <iostream>
 
@@ -55,12 +58,10 @@ void vtkQuadMesh::vtk2uLib_update() {
             << "number of quads  = " << number_of_quads << "\n"
             << "//////\n";
 
-  m_content.Points().resize(number_of_points);
+  m_content.Points().clear();
   for (int i = 0; i < number_of_points; ++i) {
     double *point = m_Poly->GetPoint(i);
-    m_content.Points()[i](0) = point[0];
-    m_content.Points()[i](1) = point[1];
-    m_content.Points()[i](2) = point[2];
+    m_content.Points().push_back(Vector3f(point[0], point[1], point[2]));
   }
 
   m_content.Quads().resize(number_of_quads);
@@ -86,11 +87,8 @@ void vtkQuadMesh::uLib2vtk_update() {
   vtkSmartPointer<vtkPoints> points = vtkSmartPointer<vtkPoints>::New();
   points->SetNumberOfPoints(number_of_points);
   for (vtkIdType i = 0; i < number_of_points; i++) {
-    double x, y, z;
-    x = m_content.Points().at(i)(0);
-    y = m_content.Points().at(i)(1);
-    z = m_content.Points().at(i)(2);
-    points->SetPoint(i, x, y, z);
+    Vector3f p = m_content.Points().at(i);
+    points->SetPoint(i, p(0), p(1), p(2));
   }
 
   vtkSmartPointer<vtkCellArray> polys = vtkSmartPointer<vtkCellArray>::New();
@@ -110,7 +108,33 @@ void vtkQuadMesh::uLib2vtk_update() {
   m_Poly->SetPoints(points);
   m_Poly->SetPolys(polys);
   m_Poly->Modified();
+}
+
+void vtkQuadMesh::contentUpdate() {
+  vtkMatrix4x4 *vmat = m_Actor->GetUserMatrix();
+  if (!vmat) {
+    vtkNew<vtkMatrix4x4> mat;
+    m_Actor->SetUserMatrix(mat);
+    vmat = mat;
+  }
+
+  Matrix4f transform = m_content.GetWorldMatrix();
+  Matrix4fToVtk(transform, vmat);
+
+  uLib2vtk_update();
+
+  m_Poly->Modified();
   m_Actor->GetMapper()->Update();
+  Puppet::Update();
+}
+
+void vtkQuadMesh::Update() {
+  vtkMatrix4x4 *vmat = m_Actor->GetUserMatrix();
+  if (!vmat) return;
+
+  Matrix4f transform = VtkToMatrix4f(vmat);
+  m_content.SetMatrix(transform);
+  m_content.Updated();
 }
 
 // -------------------------------------------------------------------------- //
@@ -121,10 +145,18 @@ vtkQuadMesh::vtkQuadMesh(vtkQuadMesh::Content &content)
       vtkSmartPointer<vtkPolyDataMapper>::New();
   mapper->SetInputData(m_Poly);
   m_Actor->SetMapper(mapper);
+
+  vtkNew<vtkMatrix4x4> vmat;
+  Matrix4fToVtk(m_content.GetWorldMatrix(), vmat);
+  m_Actor->SetUserMatrix(vmat);
+
   this->SetProp(m_Actor);
+  Object::connect(&m_content, &Content::Updated, this, &vtkQuadMesh::contentUpdate);
+  this->contentUpdate();
 }
 
 vtkQuadMesh::~vtkQuadMesh() {
+  Object::disconnect(&m_content, &Content::Updated, this, &vtkQuadMesh::contentUpdate);
   m_Poly->Delete();
   m_Actor->Delete();
 }
@@ -165,7 +197,5 @@ void vtkQuadMesh::ReadFromStlFile(const char *filename) {
 
 vtkPolyData *vtkQuadMesh::GetPolyData() const { return m_Poly; }
 
-void vtkQuadMesh::Update() { uLib2vtk_update(); }
-
 } // namespace Vtk
 } // namespace uLib

src/Vtk/Math/vtkQuadMesh.h

@@ -53,7 +53,9 @@ public:
 
   virtual class vtkPolyData *GetPolyData() const;
 
-  void Update();
+  virtual void contentUpdate();
+
+  virtual void Update();
 
 private:
   void vtk2uLib_update();

src/Vtk/Math/vtkStructuredGrid.cpp

@@ -84,9 +84,6 @@ void vtkStructuredGrid::Update() {
   }
 
   m_Content->Updated(); // Notify others (like raytracer)
-  
-  // Debug output
-  std::cout << "vtkStructuredGrid::Update matrix:\n" << transform << std::endl;
 }
 
 void vtkStructuredGrid::InstallPipe() {

src/Vtk/Math/vtkTriangleMesh.cpp

@@ -40,6 +40,9 @@
 #include <vtkPolyData.h>
 #include <vtkPolyDataMapper.h>
 
+#include <vtkMatrix4x4.h>
+#include <vtkNew.h>
+#include "Math/vtkDense.h"
 #include "Vtk/Math/vtkTriangleMesh.h"
 #include <iostream>
 
@@ -55,12 +58,10 @@ void vtkTriangleMesh::vtk2uLib_update() {
             << "number of polys  = " << number_of_triangles << "\n"
             << "//////\n";
 
-  m_content.Points().resize(number_of_points);
+  m_content.Points().clear();
   for (int i = 0; i < number_of_points; ++i) {
     double *point = m_Poly->GetPoint(i);
-    m_content.Points()[i](0) = point[0];
-    m_content.Points()[i](1) = point[1];
-    m_content.Points()[i](2) = point[2];
+    m_content.Points().push_back(Vector3f(point[0], point[1], point[2]));
   }
 
   m_content.Triangles().resize(number_of_triangles);
@@ -83,11 +84,8 @@ void vtkTriangleMesh::uLib2vtk_update() {
   vtkSmartPointer<vtkPoints> points = vtkSmartPointer<vtkPoints>::New();
   points->SetNumberOfPoints(number_of_points);
   for (vtkIdType i = 0; i < number_of_points; i++) {
-    double x, y, z;
-    x = m_content.Points().at(i)(0);
-    y = m_content.Points().at(i)(1);
-    z = m_content.Points().at(i)(2);
-    points->SetPoint(i, x, y, z);
+    Vector3f p = m_content.Points().at(i);
+    points->SetPoint(i, p(0), p(1), p(2));
   }
 
   vtkSmartPointer<vtkCellArray> polys = vtkSmartPointer<vtkCellArray>::New();
@@ -105,7 +103,33 @@ void vtkTriangleMesh::uLib2vtk_update() {
   m_Poly->SetPoints(points);
   m_Poly->SetPolys(polys);
   m_Poly->Modified();
+}
+
+void vtkTriangleMesh::contentUpdate() {
+  vtkMatrix4x4 *vmat = m_Actor->GetUserMatrix();
+  if (!vmat) {
+    vtkNew<vtkMatrix4x4> mat;
+    m_Actor->SetUserMatrix(mat);
+    vmat = mat;
+  }
+
+  Matrix4f transform = m_content.GetWorldMatrix();
+  Matrix4fToVtk(transform, vmat);
+
+  uLib2vtk_update();
+
+  m_Poly->Modified();
   m_Actor->GetMapper()->Update();
+  Puppet::Update();
+}
+
+void vtkTriangleMesh::Update() {
+  vtkMatrix4x4 *vmat = m_Actor->GetUserMatrix();
+  if (!vmat) return;
+
+  Matrix4f transform = VtkToMatrix4f(vmat);
+  m_content.SetMatrix(transform);
+  m_content.Updated();
 }
 
 // -------------------------------------------------------------------------- //
@@ -116,10 +140,18 @@ vtkTriangleMesh::vtkTriangleMesh(vtkTriangleMesh::Content &content)
       vtkSmartPointer<vtkPolyDataMapper>::New();
   mapper->SetInputData(m_Poly);
   m_Actor->SetMapper(mapper);
+
+  vtkNew<vtkMatrix4x4> vmat;
+  Matrix4fToVtk(m_content.GetWorldMatrix(), vmat);
+  m_Actor->SetUserMatrix(vmat);
+
   this->SetProp(m_Actor);
+  Object::connect(&m_content, &Content::Updated, this, &vtkTriangleMesh::contentUpdate);
+  this->contentUpdate();
 }
 
 vtkTriangleMesh::~vtkTriangleMesh() {
+  Object::disconnect(&m_content, &Content::Updated, this, &vtkTriangleMesh::contentUpdate);
   m_Poly->Delete();
   m_Actor->Delete();
 }
@@ -160,7 +192,5 @@ void vtkTriangleMesh::ReadFromStlFile(const char *filename) {
 
 vtkPolyData *vtkTriangleMesh::GetPolyData() const { return m_Poly; }
 
-void vtkTriangleMesh::Update() { uLib2vtk_update(); }
-
 } // namespace Vtk
 } // namespace uLib

src/Vtk/Math/vtkTriangleMesh.h

@@ -53,7 +53,9 @@ public:
 
   virtual class vtkPolyData *GetPolyData() const;
 
-  void Update();
+  virtual void contentUpdate();
+
+  virtual void Update();
 
 private:
   void vtk2uLib_update();

src/Vtk/vtkContainerBox.cpp

@@ -115,9 +115,6 @@ void vtkContainerBox::Update() {
   }
 
   m_Content->Updated(); // Notify change
-  
-  // Debug output
-  std::cout << "vtkContainerBox::Update matrix:\n" << transform << std::endl;
 }
 
 

src/Vtk/vtkHandlerWidget.cpp

@@ -261,7 +261,7 @@ void vtkHandlerWidget::OnMouseMove() {
   double dy = Y - this->StartEventPosition[1];
 
   if (dx == 0 && dy == 0) return;
-  std::cout << "Interaction " << this->Interaction << " dx=" << dx << " dy=" << dy << std::endl;
+  // std::cout << "Interaction " << this->Interaction << " dx=" << dx << " dy=" << dy << std::endl;
 
   // Get current gizmo properties from its actors
   vtkMatrix4x4 *gizmo_mat = m_AxesX->GetUserMatrix();
@@ -438,9 +438,6 @@ void vtkHandlerWidget::OnMouseMove() {
   this->Prop3D->Modified();
   this->UpdateGizmoPosition();
 
-  // HIGH VISIBILITY LOG
-  std::printf("--- WIDGET Interaction: %d, Mag: %f, Pos: %f %f %f\n", Interaction, mag, gpos[0], gpos[1], gpos[2]);
-
   this->InvokeEvent(::vtkCommand::InteractionEvent, nullptr);
   this->Interactor->Render();
 }

src/Vtk/vtkViewport.cpp

@@ -215,6 +215,9 @@ void Viewport::SetupPipeline(vtkRenderWindowInteractor* iren)
           self->m_HandlerWidget->SetReferenceFrame(vtkHandlerWidget::CENTER_LOCAL);
           std::cout << "Widget Frame: CENTER_LOCAL" << std::endl;
         }
+        else if (key == "s") {
+          self->ZoomSelected();
+        }
         
         iren->Render();
     });
@@ -235,6 +238,45 @@ void Viewport::ZoomAuto()
     }
 }
 
+void Viewport::ZoomSelected()
+{
+    if (!m_Renderer) return;
+
+    Puppet* selected = nullptr;
+    for (auto* p : m_Puppets) {
+        if (p->IsSelected()) {
+            selected = p;
+            break;
+        }
+    }
+    if (!selected) return;
+
+    vtkProp* prop = selected->GetProp();
+    if (!prop) return;
+
+    double* b = prop->GetBounds();
+    if (!b) return;
+    
+    double bounds[6];
+    std::copy(b, b + 6, bounds);
+    
+    if (bounds[0] > bounds[1]) return; // Invalid bounds
+
+    // Expand bounds by 1.5 from center
+    double center[3] = {(bounds[0] + bounds[1]) / 2.0, (bounds[2] + bounds[3]) / 2.0, (bounds[4] + bounds[5]) / 2.0};
+    double h_ext[3] =  {(bounds[1] - bounds[0]) / 2.0, (bounds[3] - bounds[2]) / 2.0, (bounds[5] - bounds[4]) / 2.0};
+    
+    double newBounds[6];
+    for (int i=0; i<3; ++i) {
+        newBounds[2*i]   = center[i] - 1.5 * h_ext[i];
+        newBounds[2*i+1] = center[i] + 1.5 * h_ext[i];
+    }
+    
+    m_Renderer->ResetCamera(newBounds);
+    m_Renderer->ResetCameraClippingRange();
+    this->Render();
+}
+
 void Viewport::AddPuppet(Puppet& prop)
 {
     m_Puppets.push_back(&prop);

src/Vtk/vtkViewport.h

@@ -43,6 +43,7 @@ public:
     virtual void Render() = 0;
     void Reset();
     void ZoomAuto();
+    void ZoomSelected();
 
     // Puppet / prop management
     void AddPuppet(Puppet &prop);