From e5dfb75262aee46a19aca755784b07b3570dee53 Mon Sep 17 00:00:00 2001 From: AndreaRigoni Date: Sat, 14 Mar 2026 13:23:28 +0000 Subject: [PATCH] add Normals to meshes --- src/HEP/Geant/EmitterPrimary.cpp | 99 +++++++++++++ src/HEP/Geant/EmitterPrimary.hh | 29 +++- src/Math/QuadMesh.cpp | 13 ++ src/Math/QuadMesh.h | 3 + src/Math/TriangleMesh.cpp | 15 +- src/Math/TriangleMesh.h | 3 + src/Python/math_bindings.cpp | 8 +- src/Vtk/Math/CMakeLists.txt | 2 + src/Vtk/Math/testing/CMakeLists.txt | 1 + src/Vtk/Math/testing/vtkQuadMeshTest.cpp | 56 ++++++++ src/Vtk/Math/vtkQuadMesh.cpp | 171 +++++++++++++++++++++++ src/Vtk/Math/vtkQuadMesh.h | 69 +++++++++ 12 files changed, 464 insertions(+), 5 deletions(-) create mode 100644 src/Vtk/Math/testing/vtkQuadMeshTest.cpp create mode 100644 src/Vtk/Math/vtkQuadMesh.cpp create mode 100644 src/Vtk/Math/vtkQuadMesh.h diff --git a/src/HEP/Geant/EmitterPrimary.cpp b/src/HEP/Geant/EmitterPrimary.cpp index 41c3989..43d0efe 100644 --- a/src/HEP/Geant/EmitterPrimary.cpp +++ b/src/HEP/Geant/EmitterPrimary.cpp @@ -47,4 +47,103 @@ void EmitterPrimary::GeneratePrimaries(G4Event *anEvent) { // dell'energia. fParticleGun->GeneratePrimaryVertex(anEvent); +} + +// -------------------------------------------------------------------------- // + +QuadMeshEmitterPrimary::QuadMeshEmitterPrimary() + : EmitterPrimary(), m_Mesh(nullptr), m_TotalArea(0.0) { +} + +QuadMeshEmitterPrimary::~QuadMeshEmitterPrimary() { +} + +void QuadMeshEmitterPrimary::SetMesh(uLib::QuadMesh *mesh) { + m_Mesh = mesh; + CalculateAreas(); +} + +void QuadMeshEmitterPrimary::CalculateAreas() { + if (!m_Mesh) return; + 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)]; + + double a1 = 0.5 * (v1 - v0).cross(v2 - v0).norm(); + double a2 = 0.5 * (v2 - v0).cross(v3 - v0).norm(); + m_TotalArea += (a1 + a2); + m_CumulativeAreas.push_back(m_TotalArea); + } +} + +void QuadMeshEmitterPrimary::GeneratePrimaries(G4Event *anEvent) { + if (!m_Mesh || m_TotalArea <= 0.0) return; + + // 1. Choose a quad + double r = G4UniformRand() * m_TotalArea; + auto it = std::lower_bound(m_CumulativeAreas.begin(), m_CumulativeAreas.end(), r); + 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)]; + + // 2. Choose a point on the quad + double a1 = 0.5 * (v1 - v0).cross(v2 - v0).norm(); + double a2 = 0.5 * (v2 - v0).cross(v3 - v0).norm(); + + G4ThreeVector pos; + uLib::Vector3f normal = m_Mesh->GetNormal(quadIdx); + + if (G4UniformRand() < a1 / (a1 + a2)) { + double u = G4UniformRand(); + double v = G4UniformRand(); + if (u + v > 1.0) { u = 1.0 - u; v = 1.0 - v; } + uLib::Vector3f p = v0 + u * (v1 - v0) + v * (v2 - v0); + pos.set(p(0), p(1), p(2)); + } else { + double u = G4UniformRand(); + double v = G4UniformRand(); + if (u + v > 1.0) { u = 1.0 - u; v = 1.0 - v; } + uLib::Vector3f p = v0 + u * (v2 - v0) + v * (v3 - v0); + pos.set(p(0), p(1), p(2)); + } + + // 3. Choose a direction (Cosmic Muon: cos^2(theta)) + G4ThreeVector dir; + bool accepted = false; + int tries = 0; + + while (!accepted && tries < 1000) { + tries++; + double cosTheta = std::pow(G4UniformRand(), 1.0/3.0); + double sinTheta = std::sqrt(1.0 - cosTheta * cosTheta); + double phi = 2.0 * M_PI * G4UniformRand(); + + // Incoming from above (+Z towards -Z) + dir.set(sinTheta * std::cos(phi), sinTheta * std::sin(phi), -cosTheta); + + // Filtering: pointing on the same side of the face normal + if (dir.x() * normal(0) + dir.y() * normal(1) + dir.z() * normal(2) > 0) { + accepted = true; + } + } + + if (accepted) { + fParticleGun->SetParticlePosition(pos); + fParticleGun->SetParticleMomentumDirection(dir); + // Keep energy from base class or set here if needed + fParticleGun->GeneratePrimaryVertex(anEvent); + } } \ No newline at end of file diff --git a/src/HEP/Geant/EmitterPrimary.hh b/src/HEP/Geant/EmitterPrimary.hh index 328139b..df6d8d3 100644 --- a/src/HEP/Geant/EmitterPrimary.hh +++ b/src/HEP/Geant/EmitterPrimary.hh @@ -16,11 +16,38 @@ class EmitterPrimary : public G4VUserPrimaryGeneratorAction // Metodo principale chiamato all'inizio di ogni evento virtual void GeneratePrimaries(G4Event*); - private: + protected: G4ParticleGun* fParticleGun; // Puntatore al cannone di particelle }; +#include "Math/QuadMesh.h" +#include // Added for std::vector + +namespace uLib { + class QuadMesh; +} + +class QuadMeshEmitterPrimary : public EmitterPrimary +{ + public: + QuadMeshEmitterPrimary(); + virtual ~QuadMeshEmitterPrimary(); + + // Metodo principale chiamato all'inizio di ogni evento + virtual void GeneratePrimaries(G4Event*); + + void SetMesh(uLib::QuadMesh* mesh); + + private: + uLib::QuadMesh* m_Mesh; + std::vector m_CumulativeAreas; + double m_TotalArea; + + void CalculateAreas(); +}; + + #endif \ No newline at end of file diff --git a/src/Math/QuadMesh.cpp b/src/Math/QuadMesh.cpp index f17629f..b9b18ac 100644 --- a/src/Math/QuadMesh.cpp +++ b/src/Math/QuadMesh.cpp @@ -60,5 +60,18 @@ void QuadMesh::AddQuad(const Vector4i &id) this->m_Quads.push_back(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)); + + Vector3f edge1 = v1 - v0; + Vector3f edge2 = v3 - v0; + + return edge1.cross(edge2).normalized(); +} + } diff --git a/src/Math/QuadMesh.h b/src/Math/QuadMesh.h index 1a03bf4..6889e78 100644 --- a/src/Math/QuadMesh.h +++ b/src/Math/QuadMesh.h @@ -45,6 +45,9 @@ public: inline std::vector & Points() { return this->m_Points; } inline std::vector & 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; + private: std::vector m_Points; std::vector m_Quads; diff --git a/src/Math/TriangleMesh.cpp b/src/Math/TriangleMesh.cpp index d8828ef..8c3b7df 100644 --- a/src/Math/TriangleMesh.cpp +++ b/src/Math/TriangleMesh.cpp @@ -24,8 +24,6 @@ //////////////////////////////////////////////////////////////////////////////*/ - - #include "TriangleMesh.h" @@ -65,5 +63,18 @@ void TriangleMesh::AddTriangle(const Vector3i &id) this->m_Triangles.push_back(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)); + + Vector3f edge1 = v1 - v0; + Vector3f edge2 = v2 - v0; + + return edge1.cross(edge2).normalized(); +} + } diff --git a/src/Math/TriangleMesh.h b/src/Math/TriangleMesh.h index 70d6fe9..1ebb43b 100644 --- a/src/Math/TriangleMesh.h +++ b/src/Math/TriangleMesh.h @@ -47,6 +47,9 @@ public: inline std::vector & Points() { return this->m_Points; } inline std::vector & 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; + private: std::vector m_Points; std::vector m_Triangles; diff --git a/src/Python/math_bindings.cpp b/src/Python/math_bindings.cpp index d8e0f8d..f3f6b7c 100644 --- a/src/Python/math_bindings.cpp +++ b/src/Python/math_bindings.cpp @@ -435,7 +435,9 @@ void init_math(py::module_ &m) { .def("Points", &TriangleMesh::Points, py::return_value_policy::reference_internal) .def("Triangles", &TriangleMesh::Triangles, - py::return_value_policy::reference_internal); + py::return_value_policy::reference_internal) + .def("GetTriangle", &TriangleMesh::GetTriangle) + .def("GetNormal", &TriangleMesh::GetNormal); py::class_(m, "QuadMesh") .def(py::init<>()) @@ -445,7 +447,9 @@ void init_math(py::module_ &m) { .def("Points", &QuadMesh::Points, py::return_value_policy::reference_internal) .def("Quads", &QuadMesh::Quads, - py::return_value_policy::reference_internal); + py::return_value_policy::reference_internal) + .def("GetQuad", &QuadMesh::GetQuad) + .def("GetNormal", &QuadMesh::GetNormal); py::class_(m, "VoxRaytracerRayDataElement") .def(py::init<>()) diff --git a/src/Vtk/Math/CMakeLists.txt b/src/Vtk/Math/CMakeLists.txt index 810a8a8..ca56349 100644 --- a/src/Vtk/Math/CMakeLists.txt +++ b/src/Vtk/Math/CMakeLists.txt @@ -6,6 +6,7 @@ set(MATH_SOURCES ${CMAKE_CURRENT_SOURCE_DIR}/vtkStructuredGrid.cpp ${CMAKE_CURRENT_SOURCE_DIR}/vtkTriangleMesh.cpp + ${CMAKE_CURRENT_SOURCE_DIR}/vtkQuadMesh.cpp ${CMAKE_CURRENT_SOURCE_DIR}/vtkVoxImage.cpp PARENT_SCOPE) @@ -13,6 +14,7 @@ set(MATH_HEADERS ${CMAKE_CURRENT_SOURCE_DIR}/vtkHLineRepresentation.h ${CMAKE_CURRENT_SOURCE_DIR}/vtkStructuredGrid.h ${CMAKE_CURRENT_SOURCE_DIR}/vtkTriangleMesh.h + ${CMAKE_CURRENT_SOURCE_DIR}/vtkQuadMesh.h ${CMAKE_CURRENT_SOURCE_DIR}/vtkVoxImage.h PARENT_SCOPE) diff --git a/src/Vtk/Math/testing/CMakeLists.txt b/src/Vtk/Math/testing/CMakeLists.txt index 96a9131..aae16bb 100644 --- a/src/Vtk/Math/testing/CMakeLists.txt +++ b/src/Vtk/Math/testing/CMakeLists.txt @@ -2,6 +2,7 @@ set(TESTS vtkStructuredGridTest vtkTriangleMeshTest + vtkQuadMeshTest vtkVoxImageTest ) diff --git a/src/Vtk/Math/testing/vtkQuadMeshTest.cpp b/src/Vtk/Math/testing/vtkQuadMeshTest.cpp new file mode 100644 index 0000000..1ecb57b --- /dev/null +++ b/src/Vtk/Math/testing/vtkQuadMeshTest.cpp @@ -0,0 +1,56 @@ +/*////////////////////////////////////////////////////////////////////////////// +// 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 "Vtk/Math/vtkQuadMesh.h" +#include "Math/QuadMesh.h" +#include "Vtk/uLibVtkViewer.h" + +#define BOOST_TEST_MODULE VtkQuadMeshTest +#include + +using namespace uLib; + +BOOST_AUTO_TEST_CASE(vtkQuadMeshConstruction) { + QuadMesh mesh; + + mesh.AddPoint(Vector3f(0, 0, 0)); + mesh.AddPoint(Vector3f(1, 0, 0)); + mesh.AddPoint(Vector3f(1, 1, 0)); + mesh.AddPoint(Vector3f(0, 1, 0)); + + mesh.AddQuad(Vector4i(0, 1, 2, 3)); + + Vtk::vtkQuadMesh v_mesh(mesh); + v_mesh.Update(); + + if (std::getenv("CTEST_PROJECT_NAME") == nullptr) { + Vtk::Viewer viewer; + viewer.AddPuppet(v_mesh); + viewer.Start(); + } + + BOOST_CHECK_EQUAL(mesh.Points().size(), 4u); + BOOST_CHECK_EQUAL(mesh.Quads().size(), 1u); +} diff --git a/src/Vtk/Math/vtkQuadMesh.cpp b/src/Vtk/Math/vtkQuadMesh.cpp new file mode 100644 index 0000000..dcb164b --- /dev/null +++ b/src/Vtk/Math/vtkQuadMesh.cpp @@ -0,0 +1,171 @@ +/*////////////////////////////////////////////////////////////////////////////// +// 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. + +//////////////////////////////////////////////////////////////////////////////*/ + +#ifdef HAVE_CONFIG_H +#include "config.h" +#endif + +#include + +#include +#include +#include +#include + +#include +#include +#include +#include +#include + +#include "Vtk/Math/vtkQuadMesh.h" +#include + +namespace uLib { +namespace Vtk { + +void vtkQuadMesh::vtk2uLib_update() { + vtkIdType number_of_points = m_Poly->GetNumberOfPoints(); + vtkIdType number_of_quads = m_Poly->GetNumberOfPolys(); + + std::cout << "//////\n" + << "number of points = " << number_of_points << "\n" + << "number of quads = " << number_of_quads << "\n" + << "//////\n"; + + m_content.Points().resize(number_of_points); + 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.Quads().resize(number_of_quads); + m_Poly->GetPolys()->InitTraversal(); + vtkSmartPointer idList = vtkSmartPointer::New(); + for (int i = 0; i < number_of_quads; ++i) { + m_Poly->GetPolys()->GetNextCell(idList); + if (idList->GetNumberOfIds() == 4) { + m_content.Quads()[i](0) = idList->GetId(0); + m_content.Quads()[i](1) = idList->GetId(1); + m_content.Quads()[i](2) = idList->GetId(2); + m_content.Quads()[i](3) = idList->GetId(3); + } + } + m_Poly->Modified(); + m_Actor->GetMapper()->Update(); +} + +void vtkQuadMesh::uLib2vtk_update() { + vtkIdType number_of_points = m_content.Points().size(); + vtkIdType number_of_quads = m_content.Quads().size(); + + vtkSmartPointer points = vtkSmartPointer::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); + } + + vtkSmartPointer polys = vtkSmartPointer::New(); + for (vtkIdType i = 0; i < number_of_quads; i++) { + vtkIdType a, b, c, d; + a = m_content.Quads().at(i)(0); + b = m_content.Quads().at(i)(1); + c = m_content.Quads().at(i)(2); + d = m_content.Quads().at(i)(3); + polys->InsertNextCell(4); + polys->InsertCellPoint(a); + polys->InsertCellPoint(b); + polys->InsertCellPoint(c); + polys->InsertCellPoint(d); + } + + m_Poly->SetPoints(points); + m_Poly->SetPolys(polys); + m_Poly->Modified(); + m_Actor->GetMapper()->Update(); +} + +// -------------------------------------------------------------------------- // + +vtkQuadMesh::vtkQuadMesh(vtkQuadMesh::Content &content) + : m_content(content), m_Poly(vtkPolyData::New()), m_Actor(vtkActor::New()) { + vtkSmartPointer mapper = + vtkSmartPointer::New(); + mapper->SetInputData(m_Poly); + m_Actor->SetMapper(mapper); + this->SetProp(m_Actor); +} + +vtkQuadMesh::~vtkQuadMesh() { + m_Poly->Delete(); + m_Actor->Delete(); +} + +void vtkQuadMesh::ReadFromFile(const char *filename) { + vtkSmartPointer reader = + vtkSmartPointer::New(); + reader->SetFileName(filename); + reader->Update(); + m_Poly->DeepCopy(reader->GetOutput()); + vtk2uLib_update(); +} + +void vtkQuadMesh::ReadFromXMLFile(const char *filename) { + vtkSmartPointer reader = + vtkSmartPointer::New(); + reader->SetFileName(filename); + reader->Update(); + m_Poly->DeepCopy(reader->GetOutput()); + vtk2uLib_update(); +} + +void vtkQuadMesh::ReadFromObjFile(const char *filename) { + vtkSmartPointer reader = vtkSmartPointer::New(); + reader->SetFileName(filename); + reader->Update(); + m_Poly->DeepCopy(reader->GetOutput()); + vtk2uLib_update(); +} + +void vtkQuadMesh::ReadFromStlFile(const char *filename) { + vtkSmartPointer reader = vtkSmartPointer::New(); + reader->SetFileName(filename); + reader->Update(); + m_Poly->DeepCopy(reader->GetOutput()); + vtk2uLib_update(); +} + +vtkPolyData *vtkQuadMesh::GetPolyData() const { return m_Poly; } + +void vtkQuadMesh::Update() { uLib2vtk_update(); } + +} // namespace Vtk +} // namespace uLib diff --git a/src/Vtk/Math/vtkQuadMesh.h b/src/Vtk/Math/vtkQuadMesh.h new file mode 100644 index 0000000..b238ed7 --- /dev/null +++ b/src/Vtk/Math/vtkQuadMesh.h @@ -0,0 +1,69 @@ +/*////////////////////////////////////////////////////////////////////////////// +// 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 VTKQUADMESH_H +#define VTKQUADMESH_H + +#include "Math/QuadMesh.h" +#include "Vtk/uLibVtkInterface.h" + +class vtkPolyData; +class vtkActor; + +namespace uLib { +namespace Vtk { + +class vtkQuadMesh : public Puppet, public Polydata { + typedef QuadMesh Content; + +public: + vtkQuadMesh(Content &content); + ~vtkQuadMesh(); + + void ReadFromFile(const char *filename); + + void ReadFromXMLFile(const char *filename); + + void ReadFromObjFile(const char *filename); + + void ReadFromStlFile(const char *filename); + + virtual class vtkPolyData *GetPolyData() const; + + void Update(); + +private: + void vtk2uLib_update(); + void uLib2vtk_update(); + + QuadMesh &m_content; + vtkPolyData *m_Poly; + vtkActor *m_Actor; +}; + +} // namespace Vtk +} // namespace uLib + +#endif // VTKQUADMESH_H