OpenCMT/uLib
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

feat: add CUDA raytracing benchmark and refactor VoxRaytracer::RayData to use DataAllocator for host/device memory management.

Browse Source
This commit is contained in:
AndreaRigoni
2026-03-04 17:47:18 +00:00
parent eb76521060
commit adedbcc37c
12 changed files with 1013 additions and 364 deletions
Download Patch File Download Diff File Expand all files Collapse all files

2
src/Core/CMakeLists.txt
View File

@@ -3,7 +3,7 @@ set(HEADERS Archives.h Array.h Collection.h Debug.h Export.h Function.h Macros.h
set(SOURCES Archives.cpp Debug.cpp Object.cpp Options.cpp Serializable.cpp Signal.cpp Uuid.cpp)
set(LIBRARIES Boost::program_options)
set(LIBRARIES Boost::program_options Boost::serialization)
set(libname ${PACKAGE_LIBPREFIX}Core)
set(ULIB_SHARED_LIBRARIES ${ULIB_SHARED_LIBRARIES} ${libname} PARENT_SCOPE)

2
src/Math/VoxImage.h
View File

@@ -132,7 +132,7 @@ public:
inline void SetDims(const Vector3i &size) {
this->m_Data.resize(size.prod());
BaseClass::BaseClass::SetDims(size); // FIX horrible coding style !
StructuredGrid::SetDims(size);
}
inline VoxImage<T> clipImage(const Vector3i begin, const Vector3i end) const;

71
src/Math/VoxRaytracer.cpp
View File

@@ -39,48 +39,48 @@ namespace uLib {
////////////////////////////////////////////////////////////////////////////////
void VoxRaytracer::RayData::AddElement(Id_t id, float L) {
if (m_Count >= m_Data.size()) {
size_t new_size = m_Data.size() == 0 ? 128 : m_Data.size() * 2;
m_Data.resize(new_size);
}
Element el = {id, L};
m_Data.push_back(el);
m_Data[m_Count] = el;
m_Count++;
m_TotalLength += L;
}
void VoxRaytracer::RayData::AppendRay(const VoxRaytracer::RayData &in) {
if (unlikely(!in.m_Data.size())) {
if (unlikely(in.m_Count == 0)) {
std::cout << "Warinig: PoCA on exit border!\n";
return;
} else if (unlikely(!m_Data.size())) {
m_Data = in.m_Data;
} else if (unlikely(m_Count == 0)) {
m_Data.resize(in.m_Count);
for (size_t i = 0; i < in.m_Count; ++i) {
m_Data[i] = in.m_Data[i];
}
m_Count = in.m_Count;
m_TotalLength = in.m_TotalLength;
std::cout << "Warinig: PoCA on entrance border!\n";
return;
} else {
// Opzione 1) un voxel in piu' //
if (in.m_Data.size() > 0) {
m_Data.insert(m_Data.end(), in.m_Data.begin(), in.m_Data.end());
if (in.m_Count > 0) {
if (m_Count + in.m_Count > m_Data.size()) {
m_Data.resize(m_Count + in.m_Count);
}
for (size_t i = 0; i < in.m_Count; ++i) {
m_Data[m_Count + i] = in.m_Data[i];
}
m_Count += in.m_Count;
}
// Opzione 2) merge dei voxel nel poca.
// RayData::Element &e1 = m_Data.back();
// const RayData::Element &e2 = in.m_Data.front();
// if(e1.vox_id == e2.vox_id)
// {
// m_Data.reserve(m_Data.size() + in.m_Data.size() - 1);
// e1.L += e2.L; //fix//
// m_Data.insert(m_Data.end(), in.m_Data.begin()+1,
// in.m_Data.end());
// }
// else {
// m_Data.reserve(m_Data.size() + in.m_Data.size());
// m_Data.insert(m_Data.end(), in.m_Data.begin(),
// in.m_Data.end());
// }
m_TotalLength += in.m_TotalLength;
}
}
void VoxRaytracer::RayData::PrintSelf(std::ostream &o) {
o << "Ray: total lenght " << m_TotalLength << "\n";
std::vector<Element>::iterator it;
for (it = m_Data.begin(); it < m_Data.end(); ++it)
o << "[ " << (*it).vox_id << ", " << (*it).L << "] \n";
for (size_t i = 0; i < m_Count; ++i)
o << "[ " << m_Data[i].vox_id << ", " << m_Data[i].L << "] \n";
}
////////////////////////////////////////////////////////////////////////////////
@@ -144,14 +144,21 @@ VoxRaytracer::RayData
VoxRaytracer::TraceBetweenPoints(const HPoint3f &in,
const HPoint3f &out) const {
RayData ray;
// get the local points and the direction vector
// local to image means in the normalized voxel space where the size
// of the voxel is 1 in all dimensions
Vector4f pt1 = m_Image->GetLocalPoint(in);
Vector4f pt2 = m_Image->GetLocalPoint(out);
Vector4f s = pt2 - pt1;
// l is the total length of the ray in normalized voxel space
float l = s.head(3).norm();
// L is the length of the ray between two grid lines in grid
Vector3f L(l / s(0), l / s(1), l / s(2));
// Vector3f scale; // FIXXX
// Vector3f scale; // TODO: FIX Scaling
// 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();
@@ -174,21 +181,23 @@ VoxRaytracer::TraceBetweenPoints(const HPoint3f &in,
float d;
while (l > 0) {
// find which is the minimum of the offsets to the next grid line
// it will be also the actual normalized voxel ray length
d = offset.minCoeff(&id);
// see if the voxel is inside the grid (we are still inside image)
if (m_Image->IsInsideGrid(vid)) {
// add the voxel to the ray with mapping id and length scaled
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);
// }
// move to the next voxel
vid(id) += (int)fast_sign(s(id));
// update the remaining length
l -= d;
// update the offsets
offset.array() -= d;
offset(id) = fmin(L(id), l);
}

25
src/Math/VoxRaytracer.h
View File

@@ -26,6 +26,7 @@
#ifndef VOXRAYTRACER_H
#define VOXRAYTRACER_H
#include "Math/DataAllocator.h"
#include <math.h>
#include <vector>
@@ -43,7 +44,7 @@ class VoxRaytracer {
public:
class RayData {
public:
RayData() : m_TotalLength(0) {}
RayData() : m_TotalLength(0), m_Count(0) {}
struct Element {
Id_t vox_id;
@@ -55,15 +56,24 @@ public:
void AppendRay(const RayData &in);
inline const std::vector<Element> &Data() const { return this->m_Data; }
inline DataAllocator<Element> &Data() { return this->m_Data; }
inline const DataAllocator<Element> &Data() const { return this->m_Data; }
inline size_t Count() const { return this->m_Count; }
inline const Scalarf &TotalLength() const { return this->m_TotalLength; }
inline void SetCount(size_t c) { this->m_Count = c; }
inline void SetTotalLength(Scalarf tl) { this->m_TotalLength = tl; }
void PrintSelf(std::ostream &o);
private:
std::vector<Element> m_Data;
DataAllocator<Element> m_Data;
Scalarf m_TotalLength;
size_t m_Count;
};
public:
@@ -87,6 +97,15 @@ public:
template <typename VoxelT>
void AccumulateLinesCUDA(const HLine3f *lines, size_t num_lines,
VoxImage<VoxelT> &image);
void TraceLineCUDA(const HLine3f *lines, size_t num_lines, RayData *out_rays,
int max_elements_per_ray = 128,
float *kernel_time_ms = nullptr);
void TraceBetweenPointsCUDA(const HPoint3f *in_pts, const HPoint3f *out_pts,
size_t num_lines, RayData *out_rays,
int max_elements_per_ray = 128,
float *kernel_time_ms = nullptr);
#endif
private:

410
src/Math/VoxRaytracerCUDA.hpp
View File

@@ -131,6 +131,416 @@ void VoxRaytracer::AccumulateLinesCUDA(const HLine3f *lines, size_t num_lines,
cudaFree(d_inv_world);
}
#ifdef __CUDACC__
__global__ void TraceBetweenPointsKernel(
const float *in_pts_data, const float *out_pts_data, int num_lines,
VoxRaytracer::RayData::Element **d_out_elements, size_t *d_out_counts,
float *d_out_lengths, int max_elements, int dim0, int dim1, int dim2,
const float *inv_world_matrix_data, float scale0, float scale1,
float scale2, int inc0, int inc1, int inc2) {
int idx = blockIdx.x * blockDim.x + threadIdx.x;
if (idx >= num_lines)
return;
VoxRaytracer::RayData::Element *ray_out = d_out_elements[idx];
size_t count = 0;
float tot_len = 0.0f;
const float *in_ptr = &in_pts_data[idx * 4];
const float *out_ptr = &out_pts_data[idx * 4];
float pt1[4] = {0, 0, 0, 0}, pt2[4] = {0, 0, 0, 0};
for (int i = 0; i < 4; ++i) {
for (int j = 0; j < 4; ++j) {
float m_val = inv_world_matrix_data[i + j * 4];
pt1[i] += m_val * in_ptr[j];
pt2[i] += m_val * out_ptr[j];
}
}
float s[4];
for (int i = 0; i < 4; ++i)
s[i] = pt2[i] - pt1[i];
float l = sqrtf(s[0] * s[0] + s[1] * s[1] + s[2] * s[2]);
if (l == 0) {
d_out_counts[idx] = count;
d_out_lengths[idx] = tot_len;
return;
}
float L[3];
L[0] = fabsf(l / s[0]);
L[1] = fabsf(l / s[1]);
L[2] = fabsf(l / s[2]);
float offset[3];
for (int i = 0; i < 3; ++i) {
float fpt_i = floorf(pt1[i]);
offset[i] = (s[i] >= 0) ? (1.0f - (pt1[i] - fpt_i)) : (pt1[i] - fpt_i);
offset[i] = fabsf(offset[i] * L[i]);
}
int vid[3] = {(int)floorf(pt1[0]), (int)floorf(pt1[1]), (int)floorf(pt1[2])};
int vid_out[3] = {(int)floorf(pt2[0]), (int)floorf(pt2[1]),
(int)floorf(pt2[2])};
float scale_arr[3] = {scale0, scale1, scale2};
if (vid[0] == vid_out[0] && vid[1] == vid_out[1] && vid[2] == vid_out[2]) {
if (vid[0] >= 0 && vid[0] < dim0 && vid[1] >= 0 && vid[1] < dim1 &&
vid[2] >= 0 && vid[2] < dim2) {
if (count < max_elements) {
int map_id = vid[0] * inc0 + vid[1] * inc1 + vid[2] * inc2;
ray_out[count].vox_id = map_id;
ray_out[count].L = l;
tot_len += l;
count++;
}
}
d_out_counts[idx] = count;
d_out_lengths[idx] = tot_len;
return;
}
int id;
float d;
while (l > 0) {
d = offset[0];
id = 0;
if (offset[1] < d) {
d = offset[1];
id = 1;
}
if (offset[2] < d) {
d = offset[2];
id = 2;
}
if (vid[0] >= 0 && vid[0] < dim0 && vid[1] >= 0 && vid[1] < dim1 &&
vid[2] >= 0 && vid[2] < dim2) {
if (count < max_elements) {
int map_id = vid[0] * inc0 + vid[1] * inc1 + vid[2] * inc2;
ray_out[count].vox_id = map_id;
ray_out[count].L = d * scale_arr[id];
tot_len += d * scale_arr[id];
count++;
}
}
float sign_s = (s[id] >= 0) ? 1.0f : -1.0f;
vid[id] += (int)sign_s;
l -= d;
offset[0] -= d;
offset[1] -= d;
offset[2] -= d;
offset[id] = fminf(L[id], l);
}
d_out_counts[idx] = count;
d_out_lengths[idx] = tot_len;
}
__global__ void TraceLineKernel(const float *lines_data, int num_lines,
VoxRaytracer::RayData::Element **d_out_elements,
size_t *d_out_counts, float *d_out_lengths,
int max_elements, int dim0, int dim1, int dim2,
const float *inv_world_matrix_data,
float scale0, float scale1, float scale2,
int inc0, int inc1, int inc2) {
int idx = blockIdx.x * blockDim.x + threadIdx.x;
if (idx >= num_lines)
return;
VoxRaytracer::RayData::Element *ray_out = d_out_elements[idx];
size_t count = 0;
float tot_len = 0.0f;
const float *line_ptr = &lines_data[idx * 8];
float o_vec[4] = {line_ptr[0], line_ptr[1], line_ptr[2], line_ptr[3]};
float d_vec[4] = {line_ptr[4], line_ptr[5], line_ptr[6], line_ptr[7]};
float pt[4] = {0, 0, 0, 0}, s[4] = {0, 0, 0, 0};
for (int i = 0; i < 4; ++i) {
for (int j = 0; j < 4; ++j) {
float m_val = inv_world_matrix_data[i + j * 4];
pt[i] += m_val * o_vec[j];
s[i] += m_val * d_vec[j];
}
}
float l = sqrtf(s[0] * s[0] + s[1] * s[1] + s[2] * s[2]);
if (l == 0) {
d_out_counts[idx] = count;
d_out_lengths[idx] = tot_len;
return;
}
float L[3];
L[0] = fabsf(l / s[0]);
L[1] = fabsf(l / s[1]);
L[2] = fabsf(l / s[2]);
float offset[3];
for (int i = 0; i < 3; ++i) {
float fpt_i = floorf(pt[i]);
offset[i] = (s[i] >= 0) ? (1.0f - (pt[i] - fpt_i)) : (pt[i] - fpt_i);
offset[i] = fabsf(offset[i] * L[i]);
}
int id;
float d;
int vid[3] = {(int)floorf(pt[0]), (int)floorf(pt[1]), (int)floorf(pt[2])};
float scale_arr[3] = {scale0, scale1, scale2};
while (vid[0] >= 0 && vid[0] < dim0 && vid[1] >= 0 && vid[1] < dim1 &&
vid[2] >= 0 && vid[2] < dim2) {
d = offset[0];
id = 0;
if (offset[1] < d) {
d = offset[1];
id = 1;
}
if (offset[2] < d) {
d = offset[2];
id = 2;
}
if (count < max_elements) {
int map_id = vid[0] * inc0 + vid[1] * inc1 + vid[2] * inc2;
ray_out[count].vox_id = map_id;
ray_out[count].L = d * scale_arr[id];
tot_len += d * scale_arr[id];
count++;
}
float sign_s = (s[id] >= 0) ? 1.0f : -1.0f;
vid[id] += (int)sign_s;
offset[0] -= d;
offset[1] -= d;
offset[2] -= d;
offset[id] = L[id];
}
d_out_counts[idx] = count;
d_out_lengths[idx] = tot_len;
}
#endif // __CUDACC__
inline void VoxRaytracer::TraceLineCUDA(const HLine3f *lines, size_t num_lines,
RayData *out_rays,
int max_elements_per_ray,
float *kernel_time_ms) {
if (num_lines == 0)
return;
float *d_lines = nullptr;
bool alloc_lines = false;
cudaPointerAttributes ptr_attr;
cudaPointerGetAttributes(&ptr_attr, lines);
if (ptr_attr.type == cudaMemoryTypeDevice) {
d_lines = (float *)lines;
} else {
alloc_lines = true;
size_t lines_size = num_lines * sizeof(HLine3f);
cudaMalloc(&d_lines, lines_size);
cudaMemcpy(d_lines, lines, lines_size, cudaMemcpyHostToDevice);
}
std::vector<RayData::Element *> h_out_elements(num_lines);
for (size_t i = 0; i < num_lines; ++i) {
out_rays[i].Data().resize(max_elements_per_ray);
out_rays[i].Data().MoveToVRAM();
h_out_elements[i] = out_rays[i].Data().GetVRAMData();
}
RayData::Element **d_out_elements;
cudaMalloc(&d_out_elements, num_lines * sizeof(RayData::Element *));
cudaMemcpy(d_out_elements, h_out_elements.data(),
num_lines * sizeof(RayData::Element *), cudaMemcpyHostToDevice);
size_t *d_out_counts;
float *d_out_lengths;
cudaMalloc(&d_out_counts, num_lines * sizeof(size_t));
cudaMalloc(&d_out_lengths, num_lines * sizeof(float));
int threadsPerBlock = 256;
int blocksPerGrid = (num_lines + threadsPerBlock - 1) / threadsPerBlock;
Vector3i dims = m_Image->GetDims();
Vector3i incs = m_Image->GetIncrements();
Matrix4f inv_world_matrix = m_Image->GetWorldMatrix().inverse();
float *d_inv_world;
cudaMalloc(&d_inv_world, 16 * sizeof(float));
cudaMemcpy(d_inv_world, inv_world_matrix.data(), 16 * sizeof(float),
cudaMemcpyHostToDevice);
#ifdef __CUDACC__
cudaEvent_t start, stop;
if (kernel_time_ms) {
cudaEventCreate(&start);
cudaEventCreate(&stop);
cudaEventRecord(start);
}
TraceLineKernel<<<blocksPerGrid, threadsPerBlock>>>(
d_lines, num_lines, d_out_elements, d_out_counts, d_out_lengths,
max_elements_per_ray, dims(0), dims(1), dims(2), d_inv_world, m_scale(0),
m_scale(1), m_scale(2), incs(0), incs(1), incs(2));
if (kernel_time_ms) {
cudaEventRecord(stop);
cudaEventSynchronize(stop);
cudaEventElapsedTime(kernel_time_ms, start, stop);
cudaEventDestroy(start);
cudaEventDestroy(stop);
} else {
cudaDeviceSynchronize();
}
cudaError_t err = cudaGetLastError();
if (err != cudaSuccess) {
std::cerr << "CUDA Error in TraceLineCUDA: " << cudaGetErrorString(err)
<< std::endl;
}
#else
std::cerr << "TraceLineKernel requires NVCC!" << std::endl;
#endif
std::vector<size_t> h_out_counts(num_lines);
std::vector<float> h_out_lengths(num_lines);
cudaMemcpy(h_out_counts.data(), d_out_counts, num_lines * sizeof(size_t),
cudaMemcpyDeviceToHost);
cudaMemcpy(h_out_lengths.data(), d_out_lengths, num_lines * sizeof(float),
cudaMemcpyDeviceToHost);
for (size_t i = 0; i < num_lines; ++i) {
out_rays[i].SetCount(h_out_counts[i]);
out_rays[i].SetTotalLength(h_out_lengths[i]);
}
if (alloc_lines) {
cudaFree(d_lines);
}
cudaFree(d_out_elements);
cudaFree(d_out_counts);
cudaFree(d_out_lengths);
cudaFree(d_inv_world);
}
inline void VoxRaytracer::TraceBetweenPointsCUDA(
const HPoint3f *in_pts, const HPoint3f *out_pts, size_t num_lines,
RayData *out_rays, int max_elements_per_ray, float *kernel_time_ms) {
if (num_lines == 0)
return;
float *d_in_pts = nullptr;
float *d_out_pts = nullptr;
bool alloc_pts = false;
cudaPointerAttributes ptr_attr;
cudaPointerGetAttributes(&ptr_attr, in_pts);
if (ptr_attr.type == cudaMemoryTypeDevice) {
d_in_pts = (float *)in_pts;
d_out_pts = (float *)out_pts;
} else {
alloc_pts = true;
size_t pts_size = num_lines * sizeof(HPoint3f);
cudaMalloc(&d_in_pts, pts_size);
cudaMalloc(&d_out_pts, pts_size);
cudaMemcpy(d_in_pts, in_pts, pts_size, cudaMemcpyHostToDevice);
cudaMemcpy(d_out_pts, out_pts, pts_size, cudaMemcpyHostToDevice);
}
std::vector<RayData::Element *> h_out_elements(num_lines);
for (size_t i = 0; i < num_lines; ++i) {
out_rays[i].Data().resize(max_elements_per_ray);
out_rays[i].Data().MoveToVRAM();
h_out_elements[i] = out_rays[i].Data().GetVRAMData();
}
RayData::Element **d_out_elements;
cudaMalloc(&d_out_elements, num_lines * sizeof(RayData::Element *));
cudaMemcpy(d_out_elements, h_out_elements.data(),
num_lines * sizeof(RayData::Element *), cudaMemcpyHostToDevice);
size_t *d_out_counts;
float *d_out_lengths;
cudaMalloc(&d_out_counts, num_lines * sizeof(size_t));
cudaMalloc(&d_out_lengths, num_lines * sizeof(float));
int threadsPerBlock = 256;
int blocksPerGrid = (num_lines + threadsPerBlock - 1) / threadsPerBlock;
Vector3i dims = m_Image->GetDims();
Vector3i incs = m_Image->GetIncrements();
Matrix4f inv_world_matrix = m_Image->GetWorldMatrix().inverse();
float *d_inv_world;
cudaMalloc(&d_inv_world, 16 * sizeof(float));
cudaMemcpy(d_inv_world, inv_world_matrix.data(), 16 * sizeof(float),
cudaMemcpyHostToDevice);
#ifdef __CUDACC__
cudaEvent_t start, stop;
if (kernel_time_ms) {
cudaEventCreate(&start);
cudaEventCreate(&stop);
cudaEventRecord(start);
}
TraceBetweenPointsKernel<<<blocksPerGrid, threadsPerBlock>>>(
d_in_pts, d_out_pts, num_lines, d_out_elements, d_out_counts,
d_out_lengths, max_elements_per_ray, dims(0), dims(1), dims(2),
d_inv_world, m_scale(0), m_scale(1), m_scale(2), incs(0), incs(1),
incs(2));
if (kernel_time_ms) {
cudaEventRecord(stop);
cudaEventSynchronize(stop);
cudaEventElapsedTime(kernel_time_ms, start, stop);
cudaEventDestroy(start);
cudaEventDestroy(stop);
} else {
cudaDeviceSynchronize();
}
cudaError_t err = cudaGetLastError();
if (err != cudaSuccess) {
std::cerr << "CUDA Error in TraceBetweenPointsCUDA: "
<< cudaGetErrorString(err) << std::endl;
}
#else
std::cerr << "TraceBetweenPointsKernel requires NVCC!" << std::endl;
#endif
std::vector<size_t> h_out_counts(num_lines);
std::vector<float> h_out_lengths(num_lines);
cudaMemcpy(h_out_counts.data(), d_out_counts, num_lines * sizeof(size_t),
cudaMemcpyDeviceToHost);
cudaMemcpy(h_out_lengths.data(), d_out_lengths, num_lines * sizeof(float),
cudaMemcpyDeviceToHost);
for (size_t i = 0; i < num_lines; ++i) {
out_rays[i].SetCount(h_out_counts[i]);
out_rays[i].SetTotalLength(h_out_lengths[i]);
}
if (alloc_pts) {
cudaFree(d_in_pts);
cudaFree(d_out_pts);
}
cudaFree(d_out_elements);
cudaFree(d_out_counts);
cudaFree(d_out_lengths);
cudaFree(d_inv_world);
}
} // namespace uLib
#endif // USE_CUDA

5
src/Math/testing/CMakeLists.txt
View File

@@ -5,6 +5,7 @@ set(TESTS
ContainerBoxTest
VoxImageTest
VoxRaytracerTest
VoxRaytracerTestExtended
StructuredDataTest
VoxImageFilterTest
PolicyTest
@@ -24,6 +25,6 @@ set(LIBRARIES
uLib_add_tests(Math)
if(USE_CUDA)
set_source_files_properties(VoxImageTest.cpp VoxImageCopyTest.cpp VoxImageFilterTest.cpp VoxRaytracerTest.cpp PROPERTIES LANGUAGE CUDA)
set_source_files_properties(VoxRaytracerTest.cpp PROPERTIES CXX_STANDARD 17 CUDA_STANDARD 17)
set_source_files_properties(VoxImageTest.cpp VoxImageCopyTest.cpp VoxImageFilterTest.cpp VoxRaytracerTest.cpp VoxRaytracerTestExtended.cpp PROPERTIES LANGUAGE CUDA)
set_source_files_properties(VoxRaytracerTest.cpp VoxRaytracerTestExtended.cpp PROPERTIES CXX_STANDARD 17 CUDA_STANDARD 17)
endif()

49
src/Math/testing/VoxRaytracerTest.cpp
View File

@@ -94,7 +94,8 @@ int main() {
Raytracer::RayData rdata =
ray.TraceBetweenPoints(HPoint3f(-3, -3, -3), HPoint3f(3, 3, 3));
for (const Raytracer::RayData::Element &el : rdata.Data()) {
for (size_t i = 0; i < rdata.Count(); ++i) {
const Raytracer::RayData::Element &el = rdata.Data()[i];
std::cout << " " << el.vox_id << " , " << el.L << "\n";
}
}
@@ -105,7 +106,7 @@ int main() {
Raytracer rt(img);
Raytracer::RayData ray = rt.TraceBetweenPoints(pt1, pt2);
TEST1(ray.Data().size() == 2);
TEST1(ray.Count() == 2);
TEST1(ray.Data().at(0).vox_id == 6);
TEST1(ray.Data().at(1).vox_id == 7);
ray.PrintSelf(std::cout);
@@ -117,7 +118,7 @@ int main() {
Raytracer rt(img);
Raytracer::RayData ray = rt.TraceBetweenPoints(pt1, pt2);
TEST1(ray.Data().size() == 2);
TEST1(ray.Count() == 2);
TEST1(ray.Data().at(0).vox_id == 6);
TEST1(ray.Data().at(1).vox_id == 4);
ray.PrintSelf(std::cout);
@@ -129,7 +130,7 @@ int main() {
Raytracer rt(img);
Raytracer::RayData ray = rt.TraceBetweenPoints(pt1, pt2);
TEST1(ray.Data().size() == 4);
TEST1(ray.Count() == 4);
TEST1(ray.Data().at(0).vox_id == 6);
TEST1(ray.Data().at(1).vox_id == 4);
TEST1(ray.Data().at(2).vox_id == 5);
@@ -141,6 +142,46 @@ int main() {
{
std::cout << "\n--- Testing CUDA Raytracer Accumulator ---\n";
Raytracer rt(img);
{
HPoint3f pt1(1, -0.5, 1);
HPoint3f pt2(1, 4.5, 1);
HPoint3f pts1[1] = {pt1};
HPoint3f pts2[1] = {pt2};
Raytracer::RayData ray_cuda[1];
rt.TraceBetweenPointsCUDA(pts1, pts2, 1, ray_cuda);
TEST1(ray_cuda[0].Count() == 2);
TEST1(ray_cuda[0].Data().at(0).vox_id == 6);
TEST1(ray_cuda[0].Data().at(1).vox_id == 7);
}
{
HPoint3f pt1(5, 1, 1);
HPoint3f pt2(-3, 1, 1);
HPoint3f pts1[1] = {pt1};
HPoint3f pts2[1] = {pt2};
Raytracer::RayData ray_cuda[1];
rt.TraceBetweenPointsCUDA(pts1, pts2, 1, ray_cuda);
TEST1(ray_cuda[0].Count() == 2);
TEST1(ray_cuda[0].Data().at(0).vox_id == 6);
TEST1(ray_cuda[0].Data().at(1).vox_id == 4);
}
{
HPoint3f pt1(1, 1, 1);
HPoint3f pt2(-1, 3, -1);
HPoint3f pts1[1] = {pt1};
HPoint3f pts2[1] = {pt2};
Raytracer::RayData ray_cuda[1];
rt.TraceBetweenPointsCUDA(pts1, pts2, 1, ray_cuda);
TEST1(ray_cuda[0].Count() == 4);
TEST1(ray_cuda[0].Data().at(0).vox_id == 6);
TEST1(ray_cuda[0].Data().at(1).vox_id == 4);
TEST1(ray_cuda[0].Data().at(2).vox_id == 5);
TEST1(ray_cuda[0].Data().at(3).vox_id == 1);
}
VoxImage<TestVoxel> img_cuda(Vector3i(4, 4, 4));
img_cuda.SetSpacing(Vector3f(2, 2, 2));
img_cuda.SetPosition(Vector3f(-4, -4, -4));

211
src/Math/testing/VoxRaytracerTestExtended.cpp Normal file
View File

@@ -0,0 +1,211 @@
/*//////////////////////////////////////////////////////////////////////////////
// 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 "Math/StructuredGrid.h"
#include "Math/VoxRaytracer.h"
#include "testing-prototype.h"
#include <chrono>
#include <iostream>
#include <random>
using namespace uLib;
typedef VoxRaytracer Raytracer;
int main() {
BEGIN_TESTING(Math VoxRaytracer Extended Benchmark);
std::cout << "\n=============================================\n";
std::cout << " VoxRaytracer CPU vs CUDA Benchmark Test\n";
std::cout << "=============================================\n\n";
// Create a 100x100x100 grid (1 million voxels)
StructuredGrid img(Vector3i(100, 100, 100));
img.SetSpacing(Vector3f(1.0f, 1.0f, 1.0f));
img.SetPosition(Vector3f(-50.0f, -50.0f, -50.0f));
Raytracer rt(img);
const size_t NUM_RAYS = 1000000;
std::cout << "Generating " << NUM_RAYS
<< " random ray pairs across a 100x100x100 grid...\n";
std::vector<HPoint3f> in_pts(NUM_RAYS);
std::vector<HPoint3f> out_pts(NUM_RAYS);
// Use a fixed seed for reproducible tests
std::random_device rd;
std::mt19937 gen(rd());
// The grid spans from -50 to 50 on each axis
std::uniform_real_distribution<float> dist(-49.9f, 49.9f);
// Pick a random face for in/out to ensure rays cross the volume
std::uniform_int_distribution<int> face_dist(0, 5);
for (size_t i = 0; i < NUM_RAYS; ++i) {
HPoint3f p1, p2;
// Generate point 1 on a random face
int f1 = face_dist(gen);
p1(0) = (f1 == 0) ? -50.0f : (f1 == 1) ? 50.0f : dist(gen);
p1(1) = (f1 == 2) ? -50.0f : (f1 == 3) ? 50.0f : dist(gen);
p1(2) = (f1 == 4) ? -50.0f : (f1 == 5) ? 50.0f : dist(gen);
p1(3) = 1.0f;
// Generate point 2 on a different face
int f2;
do {
f2 = face_dist(gen);
} while (
f1 == f2 ||
f1 / 2 ==
f2 / 2); // Avoid same face or opposite face trivially if desired
p2(0) = (f2 == 0) ? -50.0f : (f2 == 1) ? 50.0f : dist(gen);
p2(1) = (f2 == 2) ? -50.0f : (f2 == 3) ? 50.0f : dist(gen);
p2(2) = (f2 == 4) ? -50.0f : (f2 == 5) ? 50.0f : dist(gen);
p2(3) = 1.0f;
in_pts[i] = p1;
out_pts[i] = p2;
}
std::vector<Raytracer::RayData> cpu_results(NUM_RAYS);
std::cout << "\nRunning CPU Raytracing...\n";
auto start_cpu = std::chrono::high_resolution_clock::now();
for (size_t i = 0; i < NUM_RAYS; ++i) {
cpu_results[i] = rt.TraceBetweenPoints(in_pts[i], out_pts[i]);
}
auto end_cpu = std::chrono::high_resolution_clock::now();
std::chrono::duration<double, std::milli> cpu_ms = end_cpu - start_cpu;
std::cout << "CPU Execution Time: " << cpu_ms.count() << " ms\n";
#ifdef USE_CUDA
std::vector<Raytracer::RayData> cuda_results(NUM_RAYS);
int max_elements_per_ray =
400; // 100x100x100 grid max trace length usually ~300 items
std::cout << "\nPre-Allocating Data to VRAM...\n";
// Pre-allocate input and output points to VRAM
HPoint3f *d_in_pts;
HPoint3f *d_out_pts;
size_t pts_size = NUM_RAYS * sizeof(HPoint3f);
cudaMalloc(&d_in_pts, pts_size);
cudaMalloc(&d_out_pts, pts_size);
cudaMemcpy(d_in_pts, in_pts.data(), pts_size, cudaMemcpyHostToDevice);
cudaMemcpy(d_out_pts, out_pts.data(), pts_size, cudaMemcpyHostToDevice);
// Pre-allocate elements output arrays in VRAM via DataAllocator
for (size_t i = 0; i < NUM_RAYS; ++i) {
cuda_results[i].Data().resize(max_elements_per_ray);
cuda_results[i].Data().MoveToVRAM();
}
std::cout << "Running CUDA Raytracing...\n";
auto start_cuda = std::chrono::high_resolution_clock::now();
float kernel_time_ms = 0.0f;
rt.TraceBetweenPointsCUDA(d_in_pts, d_out_pts, NUM_RAYS, cuda_results.data(),
max_elements_per_ray, &kernel_time_ms);
auto end_cuda = std::chrono::high_resolution_clock::now();
std::chrono::duration<double, std::milli> cuda_ms = end_cuda - start_cuda;
// Free explicit input pointers
cudaFree(d_in_pts);
cudaFree(d_out_pts);
// Also query memory usage info
size_t free_byte;
size_t total_byte;
cudaMemGetInfo(&free_byte, &total_byte);
double free_db = (double)free_byte / (1024.0 * 1024.0);
double total_db = (double)total_byte / (1024.0 * 1024.0);
double used_db = total_db - free_db;
std::cout << "CUDA Kernel Exec Time: " << kernel_time_ms << " ms\n";
std::cout << "CUDA Total Time (API): " << cuda_ms.count() << " ms\n";
std::cout << "CUDA Total Time Spdup: " << (cpu_ms.count() / cuda_ms.count())
<< "x\n";
if (kernel_time_ms > 0.0f) {
std::cout << "CUDA Kernel Speedup : " << (cpu_ms.count() / kernel_time_ms)
<< "x\n";
}
std::cout << "CUDA VRAM Usage Est. : " << used_db << " MB out of " << total_db
<< " MB total\n";
std::cout << "\nVerifying CUDA results against CPU...\n";
size_t mismatches = 0;
for (size_t i = 0; i < NUM_RAYS; ++i) {
const auto &cpu_ray = cpu_results[i];
const auto &cuda_ray = cuda_results[i];
if (cpu_ray.Count() != cuda_ray.Count() ||
std::abs(cpu_ray.TotalLength() - cuda_ray.TotalLength()) > 1e-3) {
if (mismatches < 5) {
std::cout << "Mismatch at ray " << i
<< ": CPU count=" << cpu_ray.Count()
<< ", len=" << cpu_ray.TotalLength()
<< " vs CUDA count=" << cuda_ray.Count()
<< ", len=" << cuda_ray.TotalLength() << "\n";
}
mismatches++;
continue;
}
// Check elements
for (size_t j = 0; j < cpu_ray.Count(); ++j) {
if (cpu_ray.Data()[j].vox_id != cuda_ray.Data()[j].vox_id ||
std::abs(cpu_ray.Data()[j].L - cuda_ray.Data()[j].L) > 1e-3) {
if (mismatches < 5) {
std::cout << "Mismatch at ray " << i << ", element " << j
<< ": CPU id=" << cpu_ray.Data()[j].vox_id
<< ", L=" << cpu_ray.Data()[j].L
<< " vs CUDA id=" << cuda_ray.Data()[j].vox_id
<< ", L=" << cuda_ray.Data()[j].L << "\n";
}
mismatches++;
break;
}
}
}
if (mismatches == 0) {
std::cout << "SUCCESS! All " << NUM_RAYS
<< " rays perfectly match between CPU and CUDA.\n";
} else {
std::cout << "FAILED! " << mismatches << " rays contain mismatched data.\n";
}
TEST1(mismatches == 0);
#else
std::cout << "\nUSE_CUDA is not defined. Skipping CUDA benchmarking.\n";
#endif
std::cout << "=============================================\n";
END_TESTING
}

5
src/Vtk/CMakeLists.txt
View File

@@ -19,6 +19,11 @@ set(LIBRARIES Eigen3::Eigen
${VTK_LIBRARIES}
${PACKAGE_LIBPREFIX}Math)
if(USE_CUDA)
find_package(CUDAToolkit REQUIRED)
list(APPEND LIBRARIES CUDA::cudart)
endif()
set(libname ${PACKAGE_LIBPREFIX}Vtk)
set(ULIB_SHARED_LIBRARIES ${ULIB_SHARED_LIBRARIES} ${libname} PARENT_SCOPE)
set(ULIB_SELECTED_MODULES ${ULIB_SELECTED_MODULES} Vtk PARENT_SCOPE)

590
src/Vtk/vtkVoxRaytracerRepresentation.cpp
View File

@@ -23,21 +23,17 @@
//////////////////////////////////////////////////////////////////////////////*/
#ifdef HAVE_CONFIG_H
#include "config.h"
#endif
#include "vtkVoxRaytracerRepresentation.h"
#include "Math/VoxRaytracer.h"
//#include "vtkMuonEvent.h"
// #include "vtkMuonEvent.h"
#include "vtkMuonScatter.h"
namespace uLib {
namespace Vtk {
@@ -45,345 +41,297 @@ namespace Vtk {
////// VOX RAYTRACER REPRESENTATION ///////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////
vtkVoxRaytracerRepresentation::vtkVoxRaytracerRepresentation(Content &content)
: m_Content(&content), m_Assembly(vtkAssembly::New()),
m_Sphere1(vtkSphereSource::New()), m_Sphere2(vtkSphereSource::New()),
m_Line1(vtkLineSource::New()), m_Line2(vtkLineSource::New()),
m_Line3(vtkLineSource::New()), m_RayLine(vtkAppendPolyData::New()),
m_RayLineActor(vtkActor::New()),
m_RayRepresentation(vtkAppendPolyData::New()),
m_RayRepresentationActor(vtkActor::New()),
m_Transform(vtkTransform::New()) {
default_radius = content.GetImage()->GetSpacing()(0) / 4;
m_Sphere1->SetRadius(default_radius);
m_Sphere2->SetRadius(default_radius);
m_SelectedElement = m_RayLine;
InstallPipe();
}
vtkVoxRaytracerRepresentation::vtkVoxRaytracerRepresentation(Content &content) :
m_Content(&content),
m_Assembly(vtkAssembly::New()),
m_Sphere1(vtkSphereSource::New()),
m_Sphere2(vtkSphereSource::New()),
m_Line1(vtkLineSource::New()),
m_Line2(vtkLineSource::New()),
m_Line3(vtkLineSource::New()),
m_RayLine(vtkAppendPolyData::New()),
m_RayLineActor(vtkActor::New()),
m_RayRepresentation(vtkAppendPolyData::New()),
m_RayRepresentationActor(vtkActor::New()),
m_Transform(vtkTransform::New())
{
default_radius = content.GetImage()->GetSpacing()(0)/4;
m_Sphere1->SetRadius(default_radius);
m_Sphere2->SetRadius(default_radius);
vtkVoxRaytracerRepresentation::~vtkVoxRaytracerRepresentation() {
m_Assembly->Delete();
m_RayLine->Delete();
m_RayLineActor->Delete();
m_RayRepresentationActor->Delete();
m_Transform->Delete();
}
VoxRaytracer *vtkVoxRaytracerRepresentation::GetRaytracerAlgorithm() {
return m_Content;
}
vtkProp *vtkVoxRaytracerRepresentation::GetProp() { return m_Assembly; }
vtkPolyData *vtkVoxRaytracerRepresentation::GetPolyData() const {
std::cout << "get Raytracer polydata\n";
m_SelectedElement->Update();
return m_SelectedElement->GetOutput();
}
void vtkVoxRaytracerRepresentation::SetRepresentationElements(
vtkVoxRaytracerRepresentation::RepresentationElements el) {
switch (el) {
case Vtk::vtkVoxRaytracerRepresentation::RayElements:
m_SelectedElement = m_RayLine;
InstallPipe();
break;
case Vtk::vtkVoxRaytracerRepresentation::VoxelsElements:
m_SelectedElement = m_RayRepresentation;
break;
default:
m_SelectedElement = m_RayLine;
break;
}
}
vtkVoxRaytracerRepresentation::~vtkVoxRaytracerRepresentation()
{
m_Assembly->Delete();
m_RayLine->Delete();
m_RayLineActor->Delete();
m_RayRepresentationActor->Delete();
m_Transform->Delete();
void vtkVoxRaytracerRepresentation::SetMuon(MuonScatter &muon) {
HPoint3f pt1, pt2, src;
src = muon.LineIn().origin;
m_Content->GetEntryPoint(muon.LineIn(), pt1);
m_Sphere1->SetCenter(pt1(0), pt1(1), pt1(2));
m_Line1->SetPoint1(src(0), src(1), src(2));
m_Line1->SetPoint2(pt1(0), pt1(1), pt1(2));
HLine3f line_out = muon.LineOut();
src = line_out.origin;
float *direction = line_out.direction.data();
for (int i = 0; i < 3; ++i)
direction[i] *= -1;
m_Content->GetEntryPoint(line_out, pt2);
m_Sphere2->SetCenter(pt2(0), pt2(1), pt2(2));
m_Line2->SetPoint1(src(0), src(1), src(2));
m_Line2->SetPoint2(pt2(0), pt2(1), pt2(2));
m_Line3->SetPoint1(pt1(0), pt1(1), pt1(2));
m_Line3->SetPoint2(pt2(0), pt2(1), pt2(2));
// Create a vtkPoints object and store the points in it
vtkSmartPointer<vtkPoints> points = vtkSmartPointer<vtkPoints>::New();
points->InsertNextPoint(pt1(0), pt1(1), pt1(2));
points->InsertNextPoint(pt2(0), pt2(1), pt2(2));
// Create a cell array to store the lines in and add the lines to it
vtkSmartPointer<vtkCellArray> lines = vtkSmartPointer<vtkCellArray>::New();
vtkSmartPointer<vtkLine> line = vtkSmartPointer<vtkLine>::New();
line->GetPointIds()->SetId(0, 0);
line->GetPointIds()->SetId(1, 1);
lines->InsertNextCell(line);
// Create a polydata to store everything in
vtkSmartPointer<vtkPolyData> linesPolyData =
vtkSmartPointer<vtkPolyData>::New();
// Add the points to the dataset
linesPolyData->SetPoints(points);
// Add the lines to the dataset
linesPolyData->SetLines(lines);
m_RayLine->RemoveAllInputs();
#if VTK_MAJOR_VERSION <= 5
m_RayLine->AddInputConnection(linesPolyData->GetProducerPort());
#endif
m_RayLine->AddInputConnection(m_Line1->GetOutputPort());
m_RayLine->AddInputConnection(m_Sphere1->GetOutputPort());
m_RayLine->AddInputConnection(m_Line2->GetOutputPort());
m_RayLine->AddInputConnection(m_Sphere2->GetOutputPort());
vtkSmartPointer<vtkMatrix4x4> vmat = vtkSmartPointer<vtkMatrix4x4>::New();
Matrix4f mat = m_Content->GetImage()->GetWorldMatrix();
for (int i = 0; i < 4; ++i)
for (int j = 0; j < 4; ++j)
vmat->SetElement(i, j, mat(i, j));
m_Transform->SetMatrix(vmat);
this->SetRay(pt1, pt2);
}
VoxRaytracer *vtkVoxRaytracerRepresentation::GetRaytracerAlgorithm()
{
return m_Content;
void vtkVoxRaytracerRepresentation::SetMuon(MuonScatter &muon, HPoint3f poca) {
HPoint3f pt1, pt2, src;
src = muon.LineIn().origin;
m_Content->GetEntryPoint(muon.LineIn(), pt1);
m_Sphere1->SetCenter(pt1(0), pt1(1), pt1(2));
m_Line1->SetPoint1(src(0), src(1), src(2));
m_Line1->SetPoint2(pt1(0), pt1(1), pt1(2));
HLine3f line_out = muon.LineOut();
src = line_out.origin;
float *direction = line_out.direction.data();
for (int i = 0; i < 3; ++i)
direction[i] *= -1;
m_Content->GetEntryPoint(line_out, pt2);
m_Sphere2->SetCenter(pt2(0), pt2(1), pt2(2));
m_Line2->SetPoint1(src(0), src(1), src(2));
m_Line2->SetPoint2(pt2(0), pt2(1), pt2(2));
m_Line3->SetPoint1(pt1(0), pt1(1), pt1(2));
m_Line3->SetPoint2(pt2(0), pt2(1), pt2(2));
// Create a vtkPoints object and store the points in it
vtkSmartPointer<vtkPoints> points = vtkSmartPointer<vtkPoints>::New();
points->InsertNextPoint(pt1(0), pt1(1), pt1(2));
points->InsertNextPoint(poca(0), poca(1), poca(2));
points->InsertNextPoint(pt2(0), pt2(1), pt2(2));
// Create a cell array to store the lines in and add the lines to it
vtkSmartPointer<vtkCellArray> lines = vtkSmartPointer<vtkCellArray>::New();
for (unsigned int i = 0; i < 2; i++) {
vtkSmartPointer<vtkLine> line = vtkSmartPointer<vtkLine>::New();
line->GetPointIds()->SetId(0, i);
line->GetPointIds()->SetId(1, i + 1);
lines->InsertNextCell(line);
}
// Create a polydata to store everything in
vtkSmartPointer<vtkPolyData> linesPolyData =
vtkSmartPointer<vtkPolyData>::New();
// Add the points to the dataset
linesPolyData->SetPoints(points);
// Add the lines to the dataset
linesPolyData->SetLines(lines);
m_RayLine->RemoveAllInputs();
#if VTK_MAJOR_VERSION <= 5
m_RayLine->AddInputConnection(linesPolyData->GetProducerPort());
#endif
m_RayLine->AddInputConnection(m_Line1->GetOutputPort());
m_RayLine->AddInputConnection(m_Sphere1->GetOutputPort());
m_RayLine->AddInputConnection(m_Line2->GetOutputPort());
m_RayLine->AddInputConnection(m_Sphere2->GetOutputPort());
vtkSmartPointer<vtkSphereSource> poca_sphere =
vtkSmartPointer<vtkSphereSource>::New();
poca_sphere->SetRadius(default_radius);
poca_sphere->SetCenter(poca(0), poca(1), poca(2));
m_RayLine->AddInputConnection(poca_sphere->GetOutputPort());
vtkSmartPointer<vtkMatrix4x4> vmat = vtkSmartPointer<vtkMatrix4x4>::New();
Matrix4f mat = m_Content->GetImage()->GetWorldMatrix();
for (int i = 0; i < 4; ++i)
for (int j = 0; j < 4; ++j)
vmat->SetElement(i, j, mat(i, j));
m_Transform->SetMatrix(vmat);
if (m_Content->GetImage()->IsInsideBounds(poca))
this->SetRay(pt1, poca, pt2);
else
this->SetRay(pt1, pt2);
}
vtkProp *vtkVoxRaytracerRepresentation::GetProp()
{
return m_Assembly;
void vtkVoxRaytracerRepresentation::SetMuon(vtkMuonScatter &muon) {
HPoint3f poca = muon.GetPocaPoint();
MuonScatter &mu = muon.GetContent();
this->SetMuon(mu, poca);
}
vtkPolyData *vtkVoxRaytracerRepresentation::GetPolyData() const
{
std::cout << "get Raytracer polydata\n";
m_SelectedElement->Update();
return m_SelectedElement->GetOutput();
VoxRaytracer::RayData vtkVoxRaytracerRepresentation::GetRay() { return m_Ray; }
void vtkVoxRaytracerRepresentation::SetRay(HPoint3f in, HPoint3f out) {
m_Ray = m_Content->TraceBetweenPoints(in, out);
this->SetRay(&m_Ray);
}
void vtkVoxRaytracerRepresentation::SetRepresentationElements(vtkVoxRaytracerRepresentation::RepresentationElements el)
{
switch(el) {
case Vtk::vtkVoxRaytracerRepresentation::RayElements:
m_SelectedElement = m_RayLine;
break;
case Vtk::vtkVoxRaytracerRepresentation::VoxelsElements:
m_SelectedElement = m_RayRepresentation;
break;
default:
m_SelectedElement = m_RayLine;
break;
}
void vtkVoxRaytracerRepresentation::SetRay(HPoint3f in, HPoint3f mid,
HPoint3f out) {
m_Ray = m_Content->TraceBetweenPoints(in, mid);
m_Ray.AppendRay(m_Content->TraceBetweenPoints(mid, out));
this->SetRay(&m_Ray);
}
void vtkVoxRaytracerRepresentation::SetRay(VoxRaytracer::RayData *ray) {
vtkAppendPolyData *appender = m_RayRepresentation;
appender->RemoveAllInputs();
void vtkVoxRaytracerRepresentation::SetMuon(MuonScatter &muon)
{
HPoint3f pt1,pt2,src;
src = muon.LineIn().origin;
m_Content->GetEntryPoint(muon.LineIn(),pt1);
m_Sphere1->SetCenter(pt1(0),pt1(1),pt1(2));
m_Line1->SetPoint1(src(0),src(1),src(2));
m_Line1->SetPoint2(pt1(0),pt1(1),pt1(2));
HLine3f line_out = muon.LineOut();
src = line_out.origin;
float *direction = line_out.direction.data();
for(int i=0;i<3;++i) direction[i] *= -1;
m_Content->GetEntryPoint(line_out,pt2);
m_Sphere2->SetCenter(pt2(0),pt2(1),pt2(2));
m_Line2->SetPoint1(src(0),src(1),src(2));
m_Line2->SetPoint2(pt2(0),pt2(1),pt2(2));
m_Line3->SetPoint1(pt1(0),pt1(1),pt1(2));
m_Line3->SetPoint2(pt2(0),pt2(1),pt2(2));
// Create a vtkPoints object and store the points in it
vtkSmartPointer<vtkPoints> points =
vtkSmartPointer<vtkPoints>::New();
points->InsertNextPoint(pt1(0),pt1(1),pt1(2));
points->InsertNextPoint(pt2(0),pt2(1),pt2(2));
// Create a cell array to store the lines in and add the lines to it
vtkSmartPointer<vtkCellArray> lines =
vtkSmartPointer<vtkCellArray>::New();
vtkSmartPointer<vtkLine> line =
vtkSmartPointer<vtkLine>::New();
line->GetPointIds()->SetId(0,0);
line->GetPointIds()->SetId(1,1);
lines->InsertNextCell(line);
// Create a polydata to store everything in
vtkSmartPointer<vtkPolyData> linesPolyData =
vtkSmartPointer<vtkPolyData>::New();
// Add the points to the dataset
linesPolyData->SetPoints(points);
// Add the lines to the dataset
linesPolyData->SetLines(lines);
m_RayLine->RemoveAllInputs();
# if VTK_MAJOR_VERSION <= 5
m_RayLine->AddInputConnection(linesPolyData->GetProducerPort());
# endif
m_RayLine->AddInputConnection(m_Line1->GetOutputPort());
m_RayLine->AddInputConnection(m_Sphere1->GetOutputPort());
m_RayLine->AddInputConnection(m_Line2->GetOutputPort());
m_RayLine->AddInputConnection(m_Sphere2->GetOutputPort());
vtkSmartPointer<vtkMatrix4x4> vmat = vtkSmartPointer<vtkMatrix4x4>::New();
Matrix4f mat = m_Content->GetImage()->GetWorldMatrix();
for(int i=0; i<4; ++i)
for(int j=0; j<4; ++j)
vmat->SetElement(i,j,mat(i,j));
m_Transform->SetMatrix(vmat);
this->SetRay(pt1,pt2);
for (size_t i = 0; i < ray->Count(); ++i) {
int id = ray->Data().at(i).vox_id;
Vector3i idv = m_Content->GetImage()->UnMap(id);
vtkSmartPointer<vtkCubeSource> cube = vtkSmartPointer<vtkCubeSource>::New();
cube->SetBounds(idv(0), idv(0) + 1, idv(1), idv(1) + 1, idv(2), idv(2) + 1);
cube->Update();
#if VTK_MAJOR_VERSION <= 5
appender->AddInput(cube->GetOutput());
#endif
appender->Update();
}
}
void vtkVoxRaytracerRepresentation::SetMuon(MuonScatter &muon, HPoint3f poca)
{
HPoint3f pt1,pt2,src;
src = muon.LineIn().origin;
m_Content->GetEntryPoint(muon.LineIn(),pt1);
m_Sphere1->SetCenter(pt1(0),pt1(1),pt1(2));
m_Line1->SetPoint1(src(0),src(1),src(2));
m_Line1->SetPoint2(pt1(0),pt1(1),pt1(2));
HLine3f line_out = muon.LineOut();
src = line_out.origin;
float *direction = line_out.direction.data();
for(int i=0;i<3;++i) direction[i] *= -1;
m_Content->GetEntryPoint(line_out,pt2);
m_Sphere2->SetCenter(pt2(0),pt2(1),pt2(2));
m_Line2->SetPoint1(src(0),src(1),src(2));
m_Line2->SetPoint2(pt2(0),pt2(1),pt2(2));
m_Line3->SetPoint1(pt1(0),pt1(1),pt1(2));
m_Line3->SetPoint2(pt2(0),pt2(1),pt2(2));
// Create a vtkPoints object and store the points in it
vtkSmartPointer<vtkPoints> points =
vtkSmartPointer<vtkPoints>::New();
points->InsertNextPoint(pt1(0),pt1(1),pt1(2));
points->InsertNextPoint(poca(0),poca(1),poca(2));
points->InsertNextPoint(pt2(0),pt2(1),pt2(2));
// Create a cell array to store the lines in and add the lines to it
vtkSmartPointer<vtkCellArray> lines =
vtkSmartPointer<vtkCellArray>::New();
for(unsigned int i = 0; i < 2; i++)
{
vtkSmartPointer<vtkLine> line =
vtkSmartPointer<vtkLine>::New();
line->GetPointIds()->SetId(0,i);
line->GetPointIds()->SetId(1,i+1);
lines->InsertNextCell(line);
}
// Create a polydata to store everything in
vtkSmartPointer<vtkPolyData> linesPolyData =
vtkSmartPointer<vtkPolyData>::New();
// Add the points to the dataset
linesPolyData->SetPoints(points);
// Add the lines to the dataset
linesPolyData->SetLines(lines);
m_RayLine->RemoveAllInputs();
# if VTK_MAJOR_VERSION <= 5
m_RayLine->AddInputConnection(linesPolyData->GetProducerPort());
# endif
m_RayLine->AddInputConnection(m_Line1->GetOutputPort());
m_RayLine->AddInputConnection(m_Sphere1->GetOutputPort());
m_RayLine->AddInputConnection(m_Line2->GetOutputPort());
m_RayLine->AddInputConnection(m_Sphere2->GetOutputPort());
vtkSmartPointer<vtkSphereSource> poca_sphere =
vtkSmartPointer<vtkSphereSource>::New();
poca_sphere->SetRadius(default_radius);
poca_sphere->SetCenter(poca(0),poca(1),poca(2));
m_RayLine->AddInputConnection(poca_sphere->GetOutputPort());
vtkSmartPointer<vtkMatrix4x4> vmat = vtkSmartPointer<vtkMatrix4x4>::New();
Matrix4f mat = m_Content->GetImage()->GetWorldMatrix();
for(int i=0; i<4; ++i)
for(int j=0; j<4; ++j)
vmat->SetElement(i,j,mat(i,j));
m_Transform->SetMatrix(vmat);
if(m_Content->GetImage()->IsInsideBounds(poca))
this->SetRay(pt1,poca,pt2);
else
this->SetRay(pt1,pt2);
void vtkVoxRaytracerRepresentation::SetVoxelsColor(Vector4f rgba) {
this->SetColor(m_RayRepresentationActor, rgba);
}
void vtkVoxRaytracerRepresentation::SetMuon(vtkMuonScatter &muon)
{
HPoint3f poca = muon.GetPocaPoint();
MuonScatter &mu = muon.GetContent();
this->SetMuon(mu,poca);
void vtkVoxRaytracerRepresentation::SetRayColor(Vector4f rgba) {
this->SetColor(m_RayLineActor, rgba);
}
VoxRaytracer::RayData vtkVoxRaytracerRepresentation::GetRay()
{
return m_Ray;
void vtkVoxRaytracerRepresentation::SetColor(vtkActor *actor, Vector4f rgba) {
if (!actor)
return;
vtkProperty *pr = actor->GetProperty();
pr->SetDiffuseColor(rgba(0), rgba(1), rgba(2));
pr->SetOpacity(rgba(3));
pr->SetDiffuse(1);
}
void vtkVoxRaytracerRepresentation::SetRay(HPoint3f in, HPoint3f out)
{
m_Ray = m_Content->TraceBetweenPoints(in,out);
this->SetRay(&m_Ray);
void vtkVoxRaytracerRepresentation::InstallPipe() {
vtkSmartPointer<vtkAppendPolyData> append =
vtkSmartPointer<vtkAppendPolyData>::New();
append->AddInputConnection(m_Sphere1->GetOutputPort());
append->AddInputConnection(m_Sphere2->GetOutputPort());
append->AddInputConnection(m_Line1->GetOutputPort());
append->AddInputConnection(m_Line2->GetOutputPort());
append->Update();
vtkSmartPointer<vtkPolyDataMapper> mapper =
vtkSmartPointer<vtkPolyDataMapper>::New();
mapper->SetInputConnection(append->GetOutputPort());
mapper->Update();
vtkSmartPointer<vtkActor> actor = vtkSmartPointer<vtkActor>::New();
actor->SetMapper(mapper);
actor->GetProperty()->SetColor(0.6, 0.6, 1);
this->SetProp(actor);
mapper = vtkSmartPointer<vtkPolyDataMapper>::New();
mapper->SetInputConnection(m_RayLine->GetOutputPort());
mapper->Update();
m_RayLineActor->SetMapper(mapper);
m_RayLineActor->GetProperty()->SetColor(1, 0, 0);
this->SetProp(m_RayLineActor);
vtkSmartPointer<vtkTransformPolyDataFilter> polyfilter =
vtkSmartPointer<vtkTransformPolyDataFilter>::New();
polyfilter->SetInputConnection(m_RayRepresentation->GetOutputPort());
polyfilter->SetTransform(m_Transform);
mapper = vtkSmartPointer<vtkPolyDataMapper>::New();
mapper->SetInputConnection(polyfilter->GetOutputPort());
mapper->Update();
vtkActor *vra = m_RayRepresentationActor;
vra->SetMapper(mapper);
vra->GetProperty()->SetOpacity(0.2);
vra->GetProperty()->SetEdgeVisibility(true);
vra->GetProperty()->SetColor(0.5, 0.5, 0.5);
this->SetProp(vra);
}
void vtkVoxRaytracerRepresentation::SetRay(HPoint3f in, HPoint3f mid, HPoint3f out)
{
m_Ray = m_Content->TraceBetweenPoints(in,mid);
m_Ray.AppendRay( m_Content->TraceBetweenPoints(mid,out) );
this->SetRay(&m_Ray);
}
void vtkVoxRaytracerRepresentation::SetRay(VoxRaytracer::RayData *ray)
{
vtkAppendPolyData *appender = m_RayRepresentation;
appender->RemoveAllInputs();
for(int i=0; i<ray->Data().size(); ++i) {
int id = ray->Data().at(i).vox_id;
Vector3i idv = m_Content->GetImage()->UnMap(id);
vtkSmartPointer<vtkCubeSource> cube =
vtkSmartPointer<vtkCubeSource>::New();
cube->SetBounds(idv(0),idv(0)+1,idv(1),idv(1)+1,idv(2),idv(2)+1);
cube->Update();
# if VTK_MAJOR_VERSION <= 5
appender->AddInput(cube->GetOutput());
# endif
appender->Update();
}
}
void vtkVoxRaytracerRepresentation::SetVoxelsColor(Vector4f rgba)
{
this->SetColor(m_RayRepresentationActor,rgba);
}
void vtkVoxRaytracerRepresentation::SetRayColor(Vector4f rgba)
{
this->SetColor(m_RayLineActor,rgba);
}
void vtkVoxRaytracerRepresentation::SetColor(vtkActor *actor, Vector4f rgba)
{
if(!actor) return;
vtkProperty *pr = actor->GetProperty();
pr->SetDiffuseColor( rgba(0),
rgba(1),
rgba(2) );
pr->SetOpacity( rgba(3) );
pr->SetDiffuse(1);
}
void vtkVoxRaytracerRepresentation::InstallPipe()
{
vtkSmartPointer<vtkAppendPolyData> append =
vtkSmartPointer<vtkAppendPolyData>::New();
append->AddInputConnection(m_Sphere1->GetOutputPort());
append->AddInputConnection(m_Sphere2->GetOutputPort());
append->AddInputConnection(m_Line1->GetOutputPort());
append->AddInputConnection(m_Line2->GetOutputPort());
append->Update();
vtkSmartPointer<vtkPolyDataMapper> mapper =
vtkSmartPointer<vtkPolyDataMapper>::New();
mapper->SetInputConnection(append->GetOutputPort());
mapper->Update();
vtkSmartPointer<vtkActor> actor = vtkSmartPointer<vtkActor>::New();
actor->SetMapper(mapper);
actor->GetProperty()->SetColor(0.6,0.6,1);
this->SetProp(actor);
mapper = vtkSmartPointer<vtkPolyDataMapper>::New();
mapper->SetInputConnection(m_RayLine->GetOutputPort());
mapper->Update();
m_RayLineActor->SetMapper(mapper);
m_RayLineActor->GetProperty()->SetColor(1,0,0);
this->SetProp(m_RayLineActor);
vtkSmartPointer<vtkTransformPolyDataFilter> polyfilter =
vtkSmartPointer<vtkTransformPolyDataFilter>::New();
polyfilter->SetInputConnection(m_RayRepresentation->GetOutputPort());
polyfilter->SetTransform(m_Transform);
mapper = vtkSmartPointer<vtkPolyDataMapper>::New();
mapper->SetInputConnection(polyfilter->GetOutputPort());
mapper->Update();
vtkActor *vra = m_RayRepresentationActor;
vra->SetMapper(mapper);
vra->GetProperty()->SetOpacity(0.2);
vra->GetProperty()->SetEdgeVisibility(true);
vra->GetProperty()->SetColor(0.5,0.5,0.5);
this->SetProp(vra);
}
} // vtk
} // uLib
} // namespace Vtk
} // namespace uLib