diff --git a/.gitignore b/.gitignore index 4156b34..95b16f7 100644 --- a/.gitignore +++ b/.gitignore @@ -1,3 +1,7 @@ CMakeFiles/ build/ .cache/ +build_warnings*.log +final_build.log +cmake_configure.log +compile_commands.json \ No newline at end of file diff --git a/.vscode/settings.json b/.vscode/settings.json index ca906a4..ff5c426 100644 --- a/.vscode/settings.json +++ b/.vscode/settings.json @@ -13,7 +13,8 @@ "-D__global__=", "-D__constant__=", "-D__shared__=", - "-DUSE_CUDA" + "-DUSE_CUDA", + "-D__CUDACC__" ], "clangd.semanticHighlighting.enable": true, "clangd.arguments": [ diff --git a/src/Core/CMakeLists.txt b/src/Core/CMakeLists.txt index bbbe149..a1988fa 100644 --- a/src/Core/CMakeLists.txt +++ b/src/Core/CMakeLists.txt @@ -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) diff --git a/src/Math/VoxImage.h b/src/Math/VoxImage.h index 06dee03..7c98d89 100644 --- a/src/Math/VoxImage.h +++ b/src/Math/VoxImage.h @@ -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 clipImage(const Vector3i begin, const Vector3i end) const; diff --git a/src/Math/VoxRaytracer.cpp b/src/Math/VoxRaytracer.cpp index 074b8b0..4dab1f8 100644 --- a/src/Math/VoxRaytracer.cpp +++ b/src/Math/VoxRaytracer.cpp @@ -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::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); } diff --git a/src/Math/VoxRaytracer.h b/src/Math/VoxRaytracer.h index 204252c..39e8ea7 100644 --- a/src/Math/VoxRaytracer.h +++ b/src/Math/VoxRaytracer.h @@ -26,6 +26,7 @@ #ifndef VOXRAYTRACER_H #define VOXRAYTRACER_H +#include "Math/DataAllocator.h" #include #include @@ -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 &Data() const { return this->m_Data; } + inline DataAllocator &Data() { return this->m_Data; } + + inline const DataAllocator &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 m_Data; + DataAllocator m_Data; Scalarf m_TotalLength; + size_t m_Count; }; public: @@ -87,6 +97,15 @@ public: template void AccumulateLinesCUDA(const HLine3f *lines, size_t num_lines, VoxImage &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: diff --git a/src/Math/VoxRaytracerCUDA.hpp b/src/Math/VoxRaytracerCUDA.hpp index a9528d6..c0ffd19 100644 --- a/src/Math/VoxRaytracerCUDA.hpp +++ b/src/Math/VoxRaytracerCUDA.hpp @@ -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 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<<>>( + 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 h_out_counts(num_lines); + std::vector 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 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<<>>( + 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 h_out_counts(num_lines); + std::vector 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 diff --git a/src/Math/testing/CMakeLists.txt b/src/Math/testing/CMakeLists.txt index 00268eb..a88969b 100644 --- a/src/Math/testing/CMakeLists.txt +++ b/src/Math/testing/CMakeLists.txt @@ -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() diff --git a/src/Math/testing/VoxRaytracerTest.cpp b/src/Math/testing/VoxRaytracerTest.cpp index 3ba9086..ed03bba 100644 --- a/src/Math/testing/VoxRaytracerTest.cpp +++ b/src/Math/testing/VoxRaytracerTest.cpp @@ -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 img_cuda(Vector3i(4, 4, 4)); img_cuda.SetSpacing(Vector3f(2, 2, 2)); img_cuda.SetPosition(Vector3f(-4, -4, -4)); diff --git a/src/Math/testing/VoxRaytracerTestExtended.cpp b/src/Math/testing/VoxRaytracerTestExtended.cpp new file mode 100644 index 0000000..9696c97 --- /dev/null +++ b/src/Math/testing/VoxRaytracerTestExtended.cpp @@ -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 +#include +#include + +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 in_pts(NUM_RAYS); + std::vector 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 dist(-49.9f, 49.9f); + // Pick a random face for in/out to ensure rays cross the volume + std::uniform_int_distribution 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 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 cpu_ms = end_cpu - start_cpu; + std::cout << "CPU Execution Time: " << cpu_ms.count() << " ms\n"; + +#ifdef USE_CUDA + std::vector 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 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 +} diff --git a/src/Vtk/CMakeLists.txt b/src/Vtk/CMakeLists.txt index 9c193ea..8654e45 100644 --- a/src/Vtk/CMakeLists.txt +++ b/src/Vtk/CMakeLists.txt @@ -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) diff --git a/src/Vtk/vtkVoxRaytracerRepresentation.cpp b/src/Vtk/vtkVoxRaytracerRepresentation.cpp index 214b250..d4ffc8d 100644 --- a/src/Vtk/vtkVoxRaytracerRepresentation.cpp +++ b/src/Vtk/vtkVoxRaytracerRepresentation.cpp @@ -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 points = vtkSmartPointer::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 lines = vtkSmartPointer::New(); + + vtkSmartPointer line = vtkSmartPointer::New(); + line->GetPointIds()->SetId(0, 0); + line->GetPointIds()->SetId(1, 1); + lines->InsertNextCell(line); + + // Create a polydata to store everything in + vtkSmartPointer linesPolyData = + vtkSmartPointer::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 vmat = vtkSmartPointer::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 points = vtkSmartPointer::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 lines = vtkSmartPointer::New(); + + for (unsigned int i = 0; i < 2; i++) { + vtkSmartPointer line = vtkSmartPointer::New(); + line->GetPointIds()->SetId(0, i); + line->GetPointIds()->SetId(1, i + 1); + lines->InsertNextCell(line); + } + + // Create a polydata to store everything in + vtkSmartPointer linesPolyData = + vtkSmartPointer::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 poca_sphere = + vtkSmartPointer::New(); + poca_sphere->SetRadius(default_radius); + poca_sphere->SetCenter(poca(0), poca(1), poca(2)); + m_RayLine->AddInputConnection(poca_sphere->GetOutputPort()); + + vtkSmartPointer vmat = vtkSmartPointer::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 points = - vtkSmartPointer::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 lines = - vtkSmartPointer::New(); - - vtkSmartPointer line = - vtkSmartPointer::New(); - line->GetPointIds()->SetId(0,0); - line->GetPointIds()->SetId(1,1); - lines->InsertNextCell(line); - - // Create a polydata to store everything in - vtkSmartPointer linesPolyData = - vtkSmartPointer::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 vmat = vtkSmartPointer::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 cube = vtkSmartPointer::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 points = - vtkSmartPointer::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 lines = - vtkSmartPointer::New(); - - for(unsigned int i = 0; i < 2; i++) - { - vtkSmartPointer line = - vtkSmartPointer::New(); - line->GetPointIds()->SetId(0,i); - line->GetPointIds()->SetId(1,i+1); - lines->InsertNextCell(line); - } - - // Create a polydata to store everything in - vtkSmartPointer linesPolyData = - vtkSmartPointer::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 poca_sphere = - vtkSmartPointer::New(); - poca_sphere->SetRadius(default_radius); - poca_sphere->SetCenter(poca(0),poca(1),poca(2)); - m_RayLine->AddInputConnection(poca_sphere->GetOutputPort()); - - - vtkSmartPointer vmat = vtkSmartPointer::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 append = + vtkSmartPointer::New(); + append->AddInputConnection(m_Sphere1->GetOutputPort()); + append->AddInputConnection(m_Sphere2->GetOutputPort()); + append->AddInputConnection(m_Line1->GetOutputPort()); + append->AddInputConnection(m_Line2->GetOutputPort()); + + append->Update(); + vtkSmartPointer mapper = + vtkSmartPointer::New(); + + mapper->SetInputConnection(append->GetOutputPort()); + mapper->Update(); + + vtkSmartPointer actor = vtkSmartPointer::New(); + actor->SetMapper(mapper); + actor->GetProperty()->SetColor(0.6, 0.6, 1); + this->SetProp(actor); + + mapper = vtkSmartPointer::New(); + mapper->SetInputConnection(m_RayLine->GetOutputPort()); + mapper->Update(); + + m_RayLineActor->SetMapper(mapper); + m_RayLineActor->GetProperty()->SetColor(1, 0, 0); + this->SetProp(m_RayLineActor); + + vtkSmartPointer polyfilter = + vtkSmartPointer::New(); + + polyfilter->SetInputConnection(m_RayRepresentation->GetOutputPort()); + polyfilter->SetTransform(m_Transform); + + mapper = vtkSmartPointer::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; iData().size(); ++i) { - int id = ray->Data().at(i).vox_id; - Vector3i idv = m_Content->GetImage()->UnMap(id); - vtkSmartPointer cube = - vtkSmartPointer::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 append = - vtkSmartPointer::New(); - append->AddInputConnection(m_Sphere1->GetOutputPort()); - append->AddInputConnection(m_Sphere2->GetOutputPort()); - append->AddInputConnection(m_Line1->GetOutputPort()); - append->AddInputConnection(m_Line2->GetOutputPort()); - - append->Update(); - vtkSmartPointer mapper = - vtkSmartPointer::New(); - - mapper->SetInputConnection(append->GetOutputPort()); - mapper->Update(); - - vtkSmartPointer actor = vtkSmartPointer::New(); - actor->SetMapper(mapper); - actor->GetProperty()->SetColor(0.6,0.6,1); - this->SetProp(actor); - - mapper = vtkSmartPointer::New(); - mapper->SetInputConnection(m_RayLine->GetOutputPort()); - mapper->Update(); - - m_RayLineActor->SetMapper(mapper); - m_RayLineActor->GetProperty()->SetColor(1,0,0); - this->SetProp(m_RayLineActor); - - vtkSmartPointer polyfilter = - vtkSmartPointer::New(); - - polyfilter->SetInputConnection(m_RayRepresentation->GetOutputPort()); - polyfilter->SetTransform(m_Transform); - - mapper = vtkSmartPointer::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