add Normals to meshes

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);
+  }
 }

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 <vector> // 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<double> m_CumulativeAreas;
+    double m_TotalArea;
+
+    void CalculateAreas();
+};
+
+
 
 #endif

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();
+}
+
 
 }

src/Math/QuadMesh.h

@@ -45,6 +45,9 @@ public:
     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;
+
 private:
     std::vector<Vector3f> m_Points;
     std::vector<Vector4i> m_Quads;

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();
+}
+
 
 }

src/Math/TriangleMesh.h

@@ -47,6 +47,9 @@ public:
     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;
+
 private:
     std::vector<Vector3f> m_Points;
     std::vector<Vector3i> m_Triangles;

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_<QuadMesh>(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_<VoxRaytracer::RayData::Element>(m, "VoxRaytracerRayDataElement")
       .def(py::init<>())

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)
 

src/Vtk/Math/testing/CMakeLists.txt

@@ -2,6 +2,7 @@
 set(TESTS
     vtkStructuredGridTest
     vtkTriangleMeshTest
+    vtkQuadMeshTest
     vtkVoxImageTest
 )
 

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 <boost/test/unit_test.hpp>
+
+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);
+}

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 <vtkSmartPointer.h>
+
+#include <vtkOBJReader.h>
+#include <vtkPolyDataReader.h>
+#include <vtkSTLReader.h>
+#include <vtkXMLPolyDataReader.h>
+
+#include <vtkActor.h>
+#include <vtkCellArray.h>
+#include <vtkPoints.h>
+#include <vtkPolyData.h>
+#include <vtkPolyDataMapper.h>
+
+#include "Vtk/Math/vtkQuadMesh.h"
+#include <iostream>
+
+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<vtkIdList> idList = vtkSmartPointer<vtkIdList>::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<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);
+  }
+
+  vtkSmartPointer<vtkCellArray> polys = vtkSmartPointer<vtkCellArray>::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<vtkPolyDataMapper> mapper =
+      vtkSmartPointer<vtkPolyDataMapper>::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<vtkPolyDataReader> reader =
+      vtkSmartPointer<vtkPolyDataReader>::New();
+  reader->SetFileName(filename);
+  reader->Update();
+  m_Poly->DeepCopy(reader->GetOutput());
+  vtk2uLib_update();
+}
+
+void vtkQuadMesh::ReadFromXMLFile(const char *filename) {
+  vtkSmartPointer<vtkXMLPolyDataReader> reader =
+      vtkSmartPointer<vtkXMLPolyDataReader>::New();
+  reader->SetFileName(filename);
+  reader->Update();
+  m_Poly->DeepCopy(reader->GetOutput());
+  vtk2uLib_update();
+}
+
+void vtkQuadMesh::ReadFromObjFile(const char *filename) {
+  vtkSmartPointer<vtkOBJReader> reader = vtkSmartPointer<vtkOBJReader>::New();
+  reader->SetFileName(filename);
+  reader->Update();
+  m_Poly->DeepCopy(reader->GetOutput());
+  vtk2uLib_update();
+}
+
+void vtkQuadMesh::ReadFromStlFile(const char *filename) {
+  vtkSmartPointer<vtkSTLReader> reader = vtkSmartPointer<vtkSTLReader>::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

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