From a2bd38fc2c3403cd6eb6a7f907fe6e6ed92d5bf9 Mon Sep 17 00:00:00 2001 From: Paolo Andreetto Date: Sat, 26 Dec 2020 12:01:36 +0100 Subject: [PATCH] Beam tracer: draft implementation --- src/Math/VoxRaytracer.cpp | 129 ++++++++++++++++++++++++++------------ src/Math/VoxRaytracer.h | 12 ++-- 2 files changed, 95 insertions(+), 46 deletions(-) diff --git a/src/Math/VoxRaytracer.cpp b/src/Math/VoxRaytracer.cpp index ae87c13..7f57c32 100644 --- a/src/Math/VoxRaytracer.cpp +++ b/src/Math/VoxRaytracer.cpp @@ -25,14 +25,15 @@ #include +#include #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 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,9 +84,10 @@ void VoxRaytracer::RayData::AppendRay(const VoxRaytracer::RayData &in) void VoxRaytracer::RayData::PrintSelf(std::ostream &o) { o << "Ray: total lenght " << m_TotalLength << "\n"; - std::vector::iterator it; - for(it = m_Data.begin(); it < m_Data.end(); ++it) + for(std::vector::iterator it = m_Data.begin(); it < m_Data.end(); ++it) + { o << "[ " << (*it).vox_id << ", " << (*it).L << "] \n"; + } } @@ -92,6 +95,12 @@ void VoxRaytracer::RayData::PrintSelf(std::ostream &o) //// 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, - const HPoint3f &out) -const +VoxRaytracer::RayData VoxRaytracer::TraceBetweenPoints(const HPoint3f &in, const HPoint3f &out) const { RayData ray; Vector4f pt1 = m_Image->GetLocalPoint(in); @@ -162,53 +169,97 @@ 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)){ - ray.AddElement(m_Image->Map(vid),s.norm()); - return ray; + 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); - d = offset.minCoeff(&id); + if (m_Image->IsInsideGrid(vid)) + { + ray.AddElement(m_Image->Map(vid), d * m_scale(id)); + } - if(m_Image->IsInsideGrid(vid)){ - ray.AddElement(m_Image->Map(vid), d * m_scale(id) ); - } + vid(id) += (int)fast_sign(s(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; - offset.array() -= d; - offset(id) = fmin(L(id),l); + l -= d; + offset.array() -= d; + offset(id) = fmin(L(id),l); } 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 diff --git a/src/Math/VoxRaytracer.h b/src/Math/VoxRaytracer.h index 39ea821..eba67b4 100644 --- a/src/Math/VoxRaytracer.h +++ b/src/Math/VoxRaytracer.h @@ -63,13 +63,8 @@ public: }; - public: - 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(); - } +public: + VoxRaytracer(StructuredGrid &image); 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; }