Beam tracer: draft implementation

src/Math/VoxRaytracer.cpp

@@ -25,14 +25,15 @@
 
 
 #include <iostream>
+#include <unordered_map>
 
 #include "VoxRaytracer.h"
 #include "Utils.h"
 
-#define unlikely(expr) __builtin_expect(!!(expr), 0)
-
 inline float fast_sign(float f) { return 1 - 2 * (f < 0); }
 
+typedef std::unordered_map<uLib::Id_t, uLib::Scalarf> RayTable;
+
 namespace uLib {
 
 ////////////////////////////////////////////////////////////////////////////////
@@ -49,16 +50,17 @@ void VoxRaytracer::RayData::AddElement(Id_t id, float L)
 
 void VoxRaytracer::RayData::AppendRay(const VoxRaytracer::RayData &in)
 {
-    if (unlikely(!in.m_Data.size())) {
+    if (!in.m_Data.size())
+    {
         std::cout << "Warinig: PoCA on exit border!\n";
-        return;
     }
-    else if (unlikely(!m_Data.size())) {
+    else if (!m_Data.size())
+    {
         m_Data = in.m_Data;
         std::cout << "Warinig: PoCA on entrance border!\n";
-        return;
     }
-    else {
+    else
+    {
         // Opzione 1) un voxel in piu' //
         m_Data.reserve(m_Data.size() + in.m_Data.size());
         m_Data.insert(m_Data.end(), in.m_Data.begin(), in.m_Data.end());
@@ -82,16 +84,23 @@ void VoxRaytracer::RayData::AppendRay(const VoxRaytracer::RayData &in)
 void VoxRaytracer::RayData::PrintSelf(std::ostream &o)
 {
     o << "Ray: total lenght " << m_TotalLength << "\n";
-    std::vector<Element>::iterator it;
+    for(std::vector<Element>::iterator it = m_Data.begin(); it < m_Data.end(); ++it)
-    for(it = m_Data.begin(); it < m_Data.end(); ++it)
+    {
         o << "[ " << (*it).vox_id << ", " << (*it).L << "] \n";
     }
+}
 
 
 ////////////////////////////////////////////////////////////////////////////////
 //// RAY TRACER ////////////////////////////////////////////////////////////////
 ////////////////////////////////////////////////////////////////////////////////
 
+VoxRaytracer::VoxRaytracer(StructuredGrid &image) : m_Image(&image)
+{
+    m_scale << (m_Image->GetWorldMatrix() * Vector4f(1,0,0,0)).norm(),
+               (m_Image->GetWorldMatrix() * Vector4f(0,1,0,0)).norm(),
+               (m_Image->GetWorldMatrix() * Vector4f(0,0,1,0)).norm();
+}
 
 bool VoxRaytracer::GetEntryPoint(const HLine3f &line, HPoint3f &pt)
 {
@@ -150,9 +159,7 @@ bool VoxRaytracer::GetExitPoint(const HLine3f &line, HPoint3f &pt)
 }
 
 
-VoxRaytracer::RayData VoxRaytracer::TraceBetweenPoints(const HPoint3f &in,
+VoxRaytracer::RayData VoxRaytracer::TraceBetweenPoints(const HPoint3f &in, const HPoint3f &out) const
-                                                       const HPoint3f &out)
-const
 {
     RayData ray;
     Vector4f pt1 = m_Image->GetLocalPoint(in);
@@ -162,39 +169,33 @@ const
     float l = s.head(3).norm();
     Vector3f L(l/s(0), l/s(1), l/s(2));
 
-    // Vector3f scale; // FIXXX
-    // scale << (m_Image->GetWorldMatrix() * Vector4f(1,0,0,0)).norm(),
-    //          (m_Image->GetWorldMatrix() * Vector4f(0,1,0,0)).norm(),
-    //          (m_Image->GetWorldMatrix() * Vector4f(0,0,1,0)).norm();
-
     Vector3f offset;
-    for(int i=0;i<3;++i) offset(i) = (s(i)>=0) - (pt1(i)-floor(pt1(i))) ;
+    for(int i=0;i<3;++i)
+    {
+        offset(i) = (s(i)>=0) - (pt1(i)-floor(pt1(i)));
+    }
     offset = offset.cwiseProduct(L).cwiseAbs();
     L = L.cwiseAbs();
 
     //---- Check if the ray only crosses one voxel
     Vector3i vid = m_Image->Find(in);
-    if(vid == m_Image->Find(out)){
+    if (vid == m_Image->Find(out))
+    {
         ray.AddElement(m_Image->Map(vid),s.norm());
         return ray;
     }
 
     //---- Otherwise, loop until ray is finished
     int id; float d;
-    while(l>0){
+    while (l>0)
-
+    {
         d = offset.minCoeff(&id);
 
-      if(m_Image->IsInsideGrid(vid)){
+        if (m_Image->IsInsideGrid(vid))
+        {
             ray.AddElement(m_Image->Map(vid), d * m_scale(id));
         }
 
-      // nan check //
-      //        if(unlikely(!isFinite(d * scale(id)))) {
-      //            std:: cout << "NAN in raytracer\n";
-      //            exit(1);
-      //        }
-
         vid(id) += (int)fast_sign(s(id));
 
         l -= d;
@@ -204,11 +205,61 @@ const
     return ray;
 }
 
+
+// Rectangular beam
+VoxRaytracer::RayData VoxRaytracer::BeamBetweenPoints(const HPoint3f &in, const HPoint3f &out,
+                                                      int h_thick, int v_thick) const
+{
+    if (h_thick < 0) h_thick = 0;
+    if (v_thick < 0) v_thick = 0;
+
+    RayData ray = TraceBetweenPoints(in, out);
+    if (h_thick == 0 && v_thick == 0) return ray;
+
+    RayTable rTable = RayTable(ray.Data().size());
+    for(auto it = ray.Data().begin(); it < ray.Data().end(); ++it)
+    {
+        rTable.emplace(it->vox_id, it->L);
+    }
+
+    RayData beam;
+    for(auto it = ray.Data().begin(); it < ray.Data().end(); ++it)
+    {
+        Vector3i rPos = m_Image->UnMap(it->vox_id);
+
+        for (int k = h_thick * -1; k <= h_thick; k++)
+        {
+            for (int j = v_thick * -1; j <= v_thick; j++)
+            {
+                // TODO check data order in StructuredData
+                Vector3i offset { j, k, 0 };
+                Vector3i n_id = rPos + offset;
+                if (!m_Image->IsInsideGrid(n_id)) continue;
+
+                Id_t n_vox_id = m_Image->Map(n_id);
+
+                auto tPair = rTable.find(n_vox_id);
+                if (tPair == rTable.end())
+                {
+                    beam.AddElement(n_vox_id, 0);        //TODO Verify any condition with L==0
+                }
+                else if (tPair->second != 0) 
+                {
+                    beam.AddElement(n_vox_id, tPair->second);
+                    rTable[n_vox_id] = 0;
+                }
+            }
+        }
+    }
+    return beam;
+}
+
+
+
 // 20150528 SV for absorbed muons
 VoxRaytracer::RayData VoxRaytracer::TraceLine(const HLine3f &line) const
 {
     RayData ray;
-
     Vector4f pt = m_Image->GetLocalPoint(line.origin);
     Vector4f s = m_Image->GetLocalPoint(line.direction);
 
@@ -250,4 +301,4 @@ VoxRaytracer::RayData VoxRaytracer::TraceLine(const HLine3f &line) const
     return ray;
 }
 
-}
+} //end of namespace uLib

src/Math/VoxRaytracer.h

@@ -64,12 +64,7 @@ public:
 
 
 public:
- VoxRaytracer(StructuredGrid &image) : m_Image(&image) {
+    VoxRaytracer(StructuredGrid &image);
-      m_scale <<
-	(m_Image->GetWorldMatrix() * Vector4f(1,0,0,0)).norm(),
-	(m_Image->GetWorldMatrix() * Vector4f(0,1,0,0)).norm(),
-	(m_Image->GetWorldMatrix() * Vector4f(0,0,1,0)).norm();
-    }
 
     bool GetEntryPoint(const HLine3f &line, HPoint3f &pt);
 
@@ -77,6 +72,9 @@ public:
 
     RayData TraceBetweenPoints(const HPoint3f &in, const HPoint3f &out) const;
 
+    RayData BeamBetweenPoints(const HPoint3f &in, const HPoint3f &out,
+                              int h_thick = 0, int v_thick = 0) const;
+
     RayData TraceLine(const HLine3f &line) const;
 
     inline StructuredGrid* GetImage() const { return this->m_Image; }