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

refactor: migrate voxel data storage to DataAllocator for CUDA

Browse Source
This commit is contained in:
AndreaRigoni
2026-02-28 10:05:39 +00:00
parent 07915295cb
commit 52580d8cde
14 changed files with 1484 additions and 1022 deletions
Download Patch File Download Diff File Expand all files Collapse all files

232
src/Math/DataAllocator.h Normal file
View File

@@ -0,0 +1,232 @@
/*//////////////////////////////////////////////////////////////////////////////
// CMT Cosmic Muon Tomography project //////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////
Copyright (c) 2014, Universita' degli Studi di Padova, INFN sez. di Padova
All rights reserved
Authors: Andrea Rigoni Garola < andrea.rigoni@pd.infn.it >
------------------------------------------------------------------
This library is free software; you can redistribute it and/or
modify it under the terms of the GNU Lesser General Public
License as published by the Free Software Foundation; either
version 3.0 of the License, or (at your option) any later version.
This library is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
Lesser General Public License for more details.
You should have received a copy of the GNU Lesser General Public
License along with this library.
//////////////////////////////////////////////////////////////////////////////*/
#ifndef U_MATH_DATAALLOCATOR_H
#define U_MATH_DATAALLOCATOR_H
#include <algorithm>
#include <cstring>
#include <iostream>
#include <stdexcept>
#include <vector>
#ifdef USE_CUDA
#include <cuda_runtime.h>
#endif
namespace uLib {
enum class MemoryDevice { RAM, VRAM };
template <typename T> class DataAllocator {
public:
DataAllocator()
: m_Size(0), m_RamData(nullptr), m_VramData(nullptr),
m_Device(MemoryDevice::RAM) {}
DataAllocator(size_t size)
: m_Size(size), m_RamData(new T[size]()), m_VramData(nullptr),
m_Device(MemoryDevice::RAM) {}
DataAllocator(const DataAllocator<T> &other)
: m_Size(other.m_Size), m_RamData(nullptr), m_VramData(nullptr),
m_Device(other.m_Device) {
if (m_Size > 0) {
if (other.m_RamData) {
m_RamData = new T[m_Size];
std::memcpy(m_RamData, other.m_RamData, m_Size * sizeof(T));
}
#ifdef USE_CUDA
if (other.m_VramData) {
cudaMalloc((void **)&m_VramData, m_Size * sizeof(T));
cudaMemcpy(m_VramData, other.m_VramData, m_Size * sizeof(T),
cudaMemcpyDeviceToDevice);
}
#endif
}
}
~DataAllocator() {
if (m_RamData) {
delete[] m_RamData;
}
#ifdef USE_CUDA
if (m_VramData) {
cudaFree(m_VramData);
}
#endif
}
DataAllocator &operator=(const DataAllocator &other) {
if (this != &other) {
resize(other.m_Size);
m_Device = other.m_Device;
if (other.m_RamData) {
if (!m_RamData)
m_RamData = new T[m_Size];
std::memcpy(m_RamData, other.m_RamData, m_Size * sizeof(T));
}
#ifdef USE_CUDA
if (other.m_VramData) {
if (!m_VramData)
cudaMalloc((void **)&m_VramData, m_Size * sizeof(T));
cudaMemcpy(m_VramData, other.m_VramData, m_Size * sizeof(T),
cudaMemcpyDeviceToDevice);
}
#endif
}
return *this;
}
void MoveToRAM() {
if (m_Device == MemoryDevice::RAM)
return;
if (!m_RamData && m_Size > 0)
m_RamData = new T[m_Size]();
#ifdef USE_CUDA
if (m_VramData && m_Size > 0) {
cudaMemcpy(m_RamData, m_VramData, m_Size * sizeof(T),
cudaMemcpyDeviceToHost);
}
#endif
m_Device = MemoryDevice::RAM;
}
void MoveToVRAM() {
if (m_Device == MemoryDevice::VRAM)
return;
#ifdef USE_CUDA
if (!m_VramData && m_Size > 0) {
cudaMalloc((void **)&m_VramData, m_Size * sizeof(T));
}
if (m_RamData && m_Size > 0) {
cudaMemcpy(m_VramData, m_RamData, m_Size * sizeof(T),
cudaMemcpyHostToDevice);
}
#endif
m_Device = MemoryDevice::VRAM;
}
void resize(size_t size) {
if (m_Size == size)
return;
T *newRam = nullptr;
T *newVram = nullptr;
if (size > 0) {
newRam = new T[size]();
if (m_RamData) {
std::memcpy(newRam, m_RamData, std::min(m_Size, size) * sizeof(T));
}
#ifdef USE_CUDA
cudaMalloc((void **)&newVram, size * sizeof(T));
if (m_VramData) {
cudaMemcpy(newVram, m_VramData, std::min(m_Size, size) * sizeof(T),
cudaMemcpyDeviceToDevice);
}
#endif
}
if (m_RamData)
delete[] m_RamData;
#ifdef USE_CUDA
if (m_VramData)
cudaFree(m_VramData);
#endif
m_Size = size;
m_RamData = newRam;
m_VramData = newVram;
}
size_t size() const { return m_Size; }
T &at(size_t index) {
MoveToRAM();
if (index >= m_Size)
throw std::out_of_range("Index out of range");
return m_RamData[index];
}
const T &at(size_t index) const {
const_cast<DataAllocator *>(this)->MoveToRAM();
if (index >= m_Size)
throw std::out_of_range("Index out of range");
return m_RamData[index];
}
T &operator[](size_t index) {
MoveToRAM();
return m_RamData[index];
}
const T &operator[](size_t index) const {
const_cast<DataAllocator *>(this)->MoveToRAM();
return m_RamData[index];
}
T *data() { return (m_Device == MemoryDevice::RAM) ? m_RamData : m_VramData; }
const T *data() const {
return (m_Device == MemoryDevice::RAM) ? m_RamData : m_VramData;
}
T *GetRAMData() { return m_RamData; }
const T *GetRAMData() const { return m_RamData; }
T *GetVRAMData() { return m_VramData; }
const T *GetVRAMData() const { return m_VramData; }
MemoryDevice GetDevice() const { return m_Device; }
// Iterator support for RAM operations
T *begin() {
MoveToRAM();
return m_RamData;
}
T *end() {
MoveToRAM();
return m_RamData + m_Size;
}
const T *begin() const {
const_cast<DataAllocator *>(this)->MoveToRAM();
return m_RamData;
}
const T *end() const {
const_cast<DataAllocator *>(this)->MoveToRAM();
return m_RamData + m_Size;
}
private:
size_t m_Size;
T *m_RamData;
T *m_VramData;
MemoryDevice m_Device;
};
} // namespace uLib
#endif // U_MATH_DATAALLOCATOR_H

724
src/Math/VoxImage.h
View File

@@ -23,8 +23,6 @@
//////////////////////////////////////////////////////////////////////////////*/
#ifndef U_MATH_VOXIMAGE_H
#define U_MATH_VOXIMAGE_H
@@ -36,6 +34,8 @@
#include <stdlib.h>
#include <vector>
#include "Math/DataAllocator.h"
namespace uLib {
////////////////////////////////////////////////////////////////////////////////
@@ -46,36 +46,36 @@ namespace Abstract {
class VoxImage : public uLib::StructuredGrid {
public:
typedef uLib::StructuredGrid BaseClass;
typedef uLib::StructuredGrid BaseClass;
virtual float GetValue(const Vector3i &id) const = 0;
virtual float GetValue(const int id) const = 0;
virtual void SetValue(const Vector3i &id, float value) = 0;
virtual void SetValue(const int id, float value) = 0;
virtual float GetValue(const Vector3i &id) const = 0;
virtual float GetValue(const int id) const = 0;
virtual void SetValue(const Vector3i &id, float value) = 0;
virtual void SetValue(const int id, float value) = 0;
virtual void SetDims(const Vector3i &size) = 0;
virtual void SetDims(const Vector3i &size) = 0;
void ExportToVtk(const char *file, bool density_type = 0);
void ExportToVtk(const char *file, bool density_type = 0);
// use this function to export to VTK binary format
void ExportToVti (const char *file, bool density_type = 0, bool compressed = 0);
// this function has been deprecated in favor of ExportToVti
// but it is kept for backward compatibility and because it
// does not depend on vtk library
void ExportToVtkXml(const char *file, bool density_type = 0);
// use this function to export to VTK binary format
void ExportToVti(const char *file, bool density_type = 0,
bool compressed = 0);
int ImportFromVtk(const char *file, bool density_type = 0);
// this function has been deprecated in favor of ExportToVti
// but it is kept for backward compatibility and because it
// does not depend on vtk library
void ExportToVtkXml(const char *file, bool density_type = 0);
int ImportFromVti(const char *file, bool density_type = 0);
int ImportFromVtk(const char *file, bool density_type = 0);
int ImportFromVti(const char *file, bool density_type = 0);
protected:
virtual ~VoxImage() {}
VoxImage(const Vector3i &size) : BaseClass(size) {}
virtual ~VoxImage() {}
VoxImage(const Vector3i &size) : BaseClass(size) {}
};
}
} // namespace Abstract
////////////////////////////////////////////////////////////////////////////////
// VOXEL ////////////////////////////////////////////////////////////////////
@@ -83,421 +83,415 @@ protected:
namespace Interface {
struct Voxel {
template<class Self> void check_structural() {
uLibCheckMember(Self,Value, Scalarf);
}
template <class Self> void check_structural() {
uLibCheckMember(Self, Value, Scalarf);
}
};
}
} // namespace Interface
struct Voxel {
Scalarf Value;
Scalarf Value;
};
////////////////////////////////////////////////////////////////////////////////
// VOX IMAGE /////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////
template< typename T >
class VoxImage : public Abstract::VoxImage {
template <typename T> class VoxImage : public Abstract::VoxImage {
public:
typedef Abstract::VoxImage BaseClass;
typedef Abstract::VoxImage BaseClass;
VoxImage();
VoxImage();
VoxImage(const Vector3i &size);
VoxImage(const Vector3i &size);
VoxImage(const VoxImage<T> &copy) :
BaseClass(copy)
{
this->m_Data = copy.m_Data;
VoxImage(const VoxImage<T> &copy) : BaseClass(copy) {
this->m_Data = copy.m_Data;
}
inline DataAllocator<T> &Data() { return this->m_Data; }
inline const DataAllocator<T> &ConstData() const { return m_Data; }
inline const T &At(int i) const { return m_Data.at(i); }
inline const T &At(const Vector3i &id) const { return m_Data.at(Map(id)); }
inline T &operator[](unsigned int i) { return m_Data[i]; }
inline T &operator[](const Vector3i &id) { return m_Data[Map(id)]; }
// this implements Abstract interface //
inline Scalarf GetValue(const Vector3i &id) const {
return this->At(id).Value;
}
inline Scalarf GetValue(const int id) const { return this->At(id).Value; }
inline void SetValue(const Vector3i &id, Scalarf value) {
this->operator[](id).Value = value;
}
inline void SetValue(const int id, float value) {
this->operator[](id).Value = value;
}
inline void SetDims(const Vector3i &size) {
this->m_Data.resize(size.prod());
BaseClass::BaseClass::SetDims(size); // FIX horrible coding style !
}
inline VoxImage<T> clipImage(const Vector3i begin, const Vector3i end) const;
inline VoxImage<T> clipImage(const HPoint3f begin, const HPoint3f end) const;
inline VoxImage<T> clipImage(const float density) const;
inline VoxImage<T> clipImage(const float densityMin,
const float densityMax) const;
inline VoxImage<T> maskImage(const HPoint3f begin, const HPoint3f end,
float value) const;
inline VoxImage<T> maskImage(const float threshold, float belowValue = 0,
float aboveValue = 0) const;
inline VoxImage<T> fixVoxels(const float threshold, float tolerance) const;
inline VoxImage<T> fixVoxels(const float threshold, float tolerance,
const HPoint3f begin, const HPoint3f end) const;
inline VoxImage<T> fixVoxelsAroundPlane(const float threshold,
float tolerance, const HPoint3f begin,
const HPoint3f end,
bool aboveAir) const;
inline VoxImage<T> fixVoxels(const HPoint3f begin, const HPoint3f end) const;
inline VoxImage<T> Abs() const;
inline void SelectScalarfComponent(T &element, Scalarf *scalar);
inline void InitVoxels(T t);
// MATH OPERATORS //
inline void operator*=(Scalarf scalar) {
for (unsigned int i = 0; i < m_Data.size(); ++i)
m_Data[i].Value *= scalar;
}
inline void operator+=(Scalarf scalar) {
for (unsigned int i = 0; i < m_Data.size(); ++i)
m_Data[i].Value += scalar;
}
inline void operator/=(Scalarf scalar) {
for (unsigned int i = 0; i < m_Data.size(); ++i)
m_Data[i].Value /= scalar;
}
inline void operator-=(Scalarf scalar) {
for (unsigned int i = 0; i < m_Data.size(); ++i)
m_Data[i].Value -= scalar;
}
// MATH VoxImage Operators //
template <typename S> void operator+=(VoxImage<S> &sibling) {
if (this->GetDims() != sibling.GetDims()) {
// printf("Warning when adding VoxImages: I'm NOT doing it!\n");
return;
} // WARNING! You must Warn the user!
for (unsigned int i = 0; i < m_Data.size(); ++i) {
m_Data[i].Value += sibling.At(i).Value;
}
}
inline std::vector<T> & Data() { return this->m_Data; }
inline const std::vector<T>& ConstData() const { return m_Data; }
inline const T& At(int i) const { return m_Data.at(i); }
inline const T& At(const Vector3i &id) const { return m_Data.at(Map(id)); }
inline T& operator[](unsigned int i) { return m_Data[i]; }
inline T& operator[](const Vector3i &id) { return m_Data[Map(id)]; }
// this implements Abstract interface //
inline Scalarf GetValue(const Vector3i &id) const {
return this->At(id).Value;
}
inline Scalarf GetValue(const int id) const {
return this->At(id).Value;
template <typename S> void operator-=(VoxImage<S> &sibling) {
if (this->GetDims() != sibling.GetDims()) {
// printf("Warning when subtracting VoxImages: I'm NOT doing it!\n");
return;
} // WARNING! You must Warn the user!
for (unsigned int i = 0; i < m_Data.size(); ++i) {
m_Data[i].Value -= sibling.At(i).Value;
}
}
inline void SetValue(const Vector3i &id, Scalarf value) {
this->operator [](id).Value = value;
}
inline void SetValue(const int id, float value) {
this->operator [](id).Value = value;
template <typename S> void operator*=(VoxImage<S> &sibling) {
if (this->GetDims() != sibling.GetDims()) {
// printf("Warning when multiplying VoxImages: I'm NOT doing it!\n");
return;
} // WARNING! You must Warn the user!
for (unsigned int i = 0; i < m_Data.size(); ++i) {
m_Data[i].Value *= sibling.At(i).Value;
}
}
inline void SetDims(const Vector3i &size) {
this->m_Data.resize(size.prod());
BaseClass::BaseClass::SetDims(size); // FIX horrible coding style !
}
inline VoxImage<T> clipImage(const Vector3i begin, const Vector3i end) const;
inline VoxImage<T> clipImage(const HPoint3f begin, const HPoint3f end) const;
inline VoxImage<T> clipImage(const float density) const;
inline VoxImage<T> clipImage(const float densityMin, const float densityMax) const;
inline VoxImage<T> maskImage(const HPoint3f begin, const HPoint3f end, float value) const;
inline VoxImage<T> maskImage(const float threshold, float belowValue=0, float aboveValue=0) const;
inline VoxImage<T> fixVoxels(const float threshold, float tolerance) const;
inline VoxImage<T> fixVoxels(const float threshold, float tolerance, const HPoint3f begin, const HPoint3f end) const;
inline VoxImage<T> fixVoxelsAroundPlane(const float threshold, float tolerance, const HPoint3f begin, const HPoint3f end, bool aboveAir) const;
inline VoxImage<T> fixVoxels(const HPoint3f begin, const HPoint3f end) const;
inline VoxImage<T> Abs() const;
inline void SelectScalarfComponent(T &element, Scalarf *scalar);
inline void InitVoxels(T t);
// MATH OPERATORS //
inline void operator *=(Scalarf scalar) {
for(unsigned int i = 0; i < m_Data.size(); ++i)
m_Data[i].Value *= scalar;
}
inline void operator +=(Scalarf scalar) {
for(unsigned int i = 0; i < m_Data.size(); ++i)
m_Data[i].Value += scalar;
}
inline void operator /=(Scalarf scalar) {
for(unsigned int i = 0; i < m_Data.size(); ++i)
m_Data[i].Value /= scalar;
}
inline void operator -=(Scalarf scalar) {
for(unsigned int i = 0; i < m_Data.size(); ++i)
m_Data[i].Value -= scalar;
}
// MATH VoxImage Operators //
template <typename S>
void operator +=(VoxImage<S> &sibling) {
if (this->GetDims() != sibling.GetDims()) {
//printf("Warning when adding VoxImages: I'm NOT doing it!\n");
return;
}// WARNING! You must Warn the user!
for(unsigned int i = 0; i < m_Data.size(); ++i) {
m_Data[i].Value += sibling.At(i).Value;
}
}
template <typename S>
void operator -=(VoxImage<S> &sibling) {
if (this->GetDims() != sibling.GetDims()) {
//printf("Warning when subtracting VoxImages: I'm NOT doing it!\n");
return;
}// WARNING! You must Warn the user!
for(unsigned int i = 0; i < m_Data.size(); ++i) {
m_Data[i].Value -= sibling.At(i).Value;
}
}
template <typename S>
void operator *=(VoxImage<S> &sibling) {
if (this->GetDims() != sibling.GetDims()) {
//printf("Warning when multiplying VoxImages: I'm NOT doing it!\n");
return;
}// WARNING! You must Warn the user!
for(unsigned int i = 0; i < m_Data.size(); ++i) {
m_Data[i].Value *= sibling.At(i).Value;
}
}
template <typename S>
void operator /=(VoxImage<S> &sibling) {
if (this->GetDims() != sibling.GetDims()) {
//printf("Warning when dividing VoxImages: I'm NOT doing it!\n");
return;
}// WARNING! You must Warn the user!
for(unsigned int i = 0; i < m_Data.size(); ++i) {
m_Data[i].Value /= sibling.At(i).Value;
}
template <typename S> void operator/=(VoxImage<S> &sibling) {
if (this->GetDims() != sibling.GetDims()) {
// printf("Warning when dividing VoxImages: I'm NOT doing it!\n");
return;
} // WARNING! You must Warn the user!
for (unsigned int i = 0; i < m_Data.size(); ++i) {
m_Data[i].Value /= sibling.At(i).Value;
}
}
private:
std::vector<T> m_Data;
DataAllocator<T> m_Data;
};
template<typename T>
VoxImage<T>::VoxImage() :
m_Data(0),
BaseClass(Vector3i(0,0,0))
{ Interface::IsA <T,Interface::Voxel>(); /* structural check for T */ }
template<typename T>
VoxImage<T>::VoxImage(const Vector3i &size) :
m_Data(size.prod()),
BaseClass(size)
{ Interface::IsA <T,Interface::Voxel>(); /* structural check for T */ }
template <typename T>
VoxImage<T> VoxImage<T>::clipImage(const Vector3i begin, const Vector3i end) const
{
Vector3i dim = (end-begin)+Vector3i(1,1,1);
VoxImage<T> out(*this);
out.SetDims(dim);
out.SetPosition(this->GetPosition() + this->GetSpacing().cwiseProduct(begin.cast<float>()) );
for(uint x = 0; x<dim(0); ++x )
for(uint y = 0; y<dim(1); ++y )
for(uint z = 0; z<dim(2); ++z )
{
Vector3i id = Vector3i(x,y,z);
out[id] = this->At(begin + id);
}
return out;
VoxImage<T>::VoxImage() : m_Data(0), BaseClass(Vector3i(0, 0, 0)) {
Interface::IsA<T, Interface::Voxel>(); /* structural check for T */
}
template <typename T>
VoxImage<T> VoxImage<T>::clipImage(const HPoint3f begin, const HPoint3f end) const
{
Vector3i v1 = this->Find(begin);
Vector3i v2 = this->Find(end);
return this->clipImage(v1,v2);
VoxImage<T>::VoxImage(const Vector3i &size)
: m_Data(size.prod()), BaseClass(size) {
Interface::IsA<T, Interface::Voxel>(); /* structural check for T */
}
template <typename T>
VoxImage<T> VoxImage<T>::clipImage(const float density) const
{
Vector3i v1 = this->GetDims();
Vector3i v2 = Vector3i(0,0,0);
for(uint i=0; i< this->m_Data.size(); ++i) {
if( this->GetValue(i) >= density ) {
Vector3i id = this->UnMap(i);
v1 = v1.array().min(id.array());
v2 = v2.array().max(id.array());
}
VoxImage<T> VoxImage<T>::clipImage(const Vector3i begin,
const Vector3i end) const {
Vector3i dim = (end - begin) + Vector3i(1, 1, 1);
VoxImage<T> out(*this);
out.SetDims(dim);
out.SetPosition(this->GetPosition() +
this->GetSpacing().cwiseProduct(begin.cast<float>()));
for (uint x = 0; x < dim(0); ++x)
for (uint y = 0; y < dim(1); ++y)
for (uint z = 0; z < dim(2); ++z) {
Vector3i id = Vector3i(x, y, z);
out[id] = this->At(begin + id);
}
return out;
}
template <typename T>
VoxImage<T> VoxImage<T>::clipImage(const HPoint3f begin,
const HPoint3f end) const {
Vector3i v1 = this->Find(begin);
Vector3i v2 = this->Find(end);
return this->clipImage(v1, v2);
}
template <typename T>
VoxImage<T> VoxImage<T>::clipImage(const float density) const {
Vector3i v1 = this->GetDims();
Vector3i v2 = Vector3i(0, 0, 0);
for (uint i = 0; i < this->m_Data.size(); ++i) {
if (this->GetValue(i) >= density) {
Vector3i id = this->UnMap(i);
v1 = v1.array().min(id.array());
v2 = v2.array().max(id.array());
}
return this->clipImage(v1,v2);
}
return this->clipImage(v1, v2);
}
template <typename T>
VoxImage<T> VoxImage<T>::clipImage(const float densityMin, const float densityMax) const
{
Vector3i v1 = this->GetDims();
Vector3i v2 = Vector3i(0,0,0);
for(uint i=0; i< this->m_Data.size(); ++i) {
if( this->GetValue(i) >= densityMin && this->GetValue(i) <= densityMax) {
Vector3i id = this->UnMap(i);
v1 = v1.array().min(id.array());
v2 = v2.array().max(id.array());
}
VoxImage<T> VoxImage<T>::clipImage(const float densityMin,
const float densityMax) const {
Vector3i v1 = this->GetDims();
Vector3i v2 = Vector3i(0, 0, 0);
for (uint i = 0; i < this->m_Data.size(); ++i) {
if (this->GetValue(i) >= densityMin && this->GetValue(i) <= densityMax) {
Vector3i id = this->UnMap(i);
v1 = v1.array().min(id.array());
v2 = v2.array().max(id.array());
}
return this->clipImage(v1,v2);
}
return this->clipImage(v1, v2);
}
template <typename T>
VoxImage<T> VoxImage<T>::maskImage(const HPoint3f begin, const HPoint3f end, float value) const
{
VoxImage<T> out(*this);
out.SetDims(this->GetDims());
out.SetPosition(this->GetPosition());
VoxImage<T> VoxImage<T>::maskImage(const HPoint3f begin, const HPoint3f end,
float value) const {
VoxImage<T> out(*this);
out.SetDims(this->GetDims());
out.SetPosition(this->GetPosition());
Vector3i voxB = this->Find(begin);
Vector3i voxE = this->Find(end);
Vector3i voxB = this->Find(begin);
Vector3i voxE = this->Find(end);
Vector3i ID;
Vector3i ID;
for(int ix=voxB(0); ix<voxE(0); ix++)
for(int iy=voxB(1); iy<voxE(1); iy++)
for(int iz=voxB(2); iz<voxE(2); iz++){
ID << ix,iy,iz;
out.SetValue(ID,value*1.E-6);
}
for (int ix = voxB(0); ix < voxE(0); ix++)
for (int iy = voxB(1); iy < voxE(1); iy++)
for (int iz = voxB(2); iz < voxE(2); iz++) {
ID << ix, iy, iz;
out.SetValue(ID, value * 1.E-6);
}
return out;
return out;
}
template <typename T>
VoxImage<T> VoxImage<T>::maskImage(const float threshold, float belowValue, float aboveValue) const
{
std::cout << "VoxImage: maskImage, fixing voxels under threshold " << threshold;
if(belowValue)
std::cout << " at value " << belowValue;
else
std::cout << " at -value";
std::cout << ", voxels above threshold at value ";
if(aboveValue)
std::cout << aboveValue;
else
std::cout << "found";
VoxImage<T> VoxImage<T>::maskImage(const float threshold, float belowValue,
float aboveValue) const {
std::cout << "VoxImage: maskImage, fixing voxels under threshold "
<< threshold;
if (belowValue)
std::cout << " at value " << belowValue;
else
std::cout << " at -value";
std::cout << ", voxels above threshold at value ";
if (aboveValue)
std::cout << aboveValue;
else
std::cout << "found";
VoxImage<T> out(*this);
out.SetDims(this->GetDims());
out.SetPosition(this->GetPosition());
VoxImage<T> out(*this);
out.SetDims(this->GetDims());
out.SetPosition(this->GetPosition());
for(uint i=0; i< this->m_Data.size(); ++i) {
// skip negative voxels: they are already frozen
if( this->GetValue(i) >= 0 ){
// voxels under threshold
if( this->GetValue(i) <= threshold*1.E-6 ){
if(belowValue){
// std::cout << "vox " << i << ", " << this->GetValue(i);
// std::cout << " ----> set to " << -1.*belowValue*1.E-6 << std::endl;
out.SetValue(i,-1.*belowValue*1.E-6);}
else
out.SetValue(i,-1.*this->GetValue(i));
}
// voxels over threshold
else{
if(aboveValue)
out.SetValue(i,aboveValue*1.E-6);
else
out.SetValue(i,this->GetValue(i));
}
}
for (uint i = 0; i < this->m_Data.size(); ++i) {
// skip negative voxels: they are already frozen
if (this->GetValue(i) >= 0) {
// voxels under threshold
if (this->GetValue(i) <= threshold * 1.E-6) {
if (belowValue) {
// std::cout << "vox " << i << ", " <<
// this->GetValue(i); std::cout << " ----> set to " <<
// -1.*belowValue*1.E-6 << std::endl;
out.SetValue(i, -1. * belowValue * 1.E-6);
} else
out.SetValue(i, -1. * this->GetValue(i));
}
// voxels over threshold
else {
if (aboveValue)
out.SetValue(i, aboveValue * 1.E-6);
else
out.SetValue(i, this->GetValue(i));
}
}
return out;
}
return out;
}
template <typename T>
VoxImage<T> VoxImage<T>::fixVoxels(const float threshold, float tolerance) const
{
std::cout << "VoxImage: fixing voxels with value " << threshold << std::endl;
VoxImage<T> VoxImage<T>::fixVoxels(const float threshold,
float tolerance) const {
std::cout << "VoxImage: fixing voxels with value " << threshold << std::endl;
VoxImage<T> out(*this);
out.SetDims(this->GetDims());
out.SetPosition(this->GetPosition());
VoxImage<T> out(*this);
out.SetDims(this->GetDims());
out.SetPosition(this->GetPosition());
for(uint i=0; i< this->m_Data.size(); ++i) {
for (uint i = 0; i < this->m_Data.size(); ++i) {
// voxels around threshold
if (fabs(this->GetValue(i) - threshold * 1.E-6) < tolerance * 1.E-6) {
// std::cout << "vox " << i << ", " << this->GetValue(i);
// std::cout << " ----> set to " << -1.*this->GetValue(i) <<
// std::endl;
out.SetValue(i, -1. * this->GetValue(i));
}
}
return out;
}
template <typename T> VoxImage<T> VoxImage<T>::Abs() const {
std::cout << "VoxImage: set abs voxels value " << std::endl;
VoxImage<T> out(*this);
out.SetDims(this->GetDims());
out.SetPosition(this->GetPosition());
for (uint i = 0; i < this->m_Data.size(); ++i)
out.SetValue(i, fabs(this->GetValue(i)));
return out;
}
template <typename T>
VoxImage<T> VoxImage<T>::fixVoxels(const float threshold, float tolerance,
const HPoint3f begin,
const HPoint3f end) const {
VoxImage<T> out(*this);
out.SetDims(this->GetDims());
out.SetPosition(this->GetPosition());
Vector3i voxB = this->Find(begin);
Vector3i voxE = this->Find(end);
Vector3i ID;
for (int ix = voxB(0); ix < voxE(0); ix++)
for (int iy = voxB(1); iy < voxE(1); iy++)
for (int iz = voxB(2); iz < voxE(2); iz++) {
ID << ix, iy, iz;
// voxels around threshold
if( fabs(this->GetValue(i) - threshold*1.E-6) < tolerance* 1.E-6 ){
// std::cout << "vox " << i << ", " << this->GetValue(i);
// std::cout << " ----> set to " << -1.*this->GetValue(i) << std::endl;
out.SetValue(i,-1.*this->GetValue(i));
if (fabs(this->GetValue(ID) - threshold * 1.E-6) < tolerance * 1.E-6) {
out.SetValue(ID, -1. * this->GetValue(ID));
}
}
return out;
}
return out;
}
template <typename T>
VoxImage<T> VoxImage<T>::Abs() const
{
std::cout << "VoxImage: set abs voxels value " << std::endl;
VoxImage<T> VoxImage<T>::fixVoxels(const HPoint3f begin,
const HPoint3f end) const {
VoxImage<T> out(*this);
out.SetDims(this->GetDims());
out.SetPosition(this->GetPosition());
VoxImage<T> out(*this);
out.SetDims(this->GetDims());
out.SetPosition(this->GetPosition());
Vector3i voxB = this->Find(begin);
Vector3i voxE = this->Find(end);
for(uint i=0; i< this->m_Data.size(); ++i)
out.SetValue(i,fabs(this->GetValue(i)));
Vector3i ID;
return out;
for (int ix = voxB(0); ix < voxE(0); ix++)
for (int iy = voxB(1); iy < voxE(1); iy++)
for (int iz = voxB(2); iz < voxE(2); iz++) {
ID << ix, iy, iz;
// voxels around threshold
out.SetValue(ID, -1. * this->GetValue(ID));
}
return out;
}
template <typename T>
VoxImage<T> VoxImage<T>::fixVoxels( const float threshold, float tolerance, const HPoint3f begin, const HPoint3f end) const
{
VoxImage<T> out(*this);
out.SetDims(this->GetDims());
out.SetPosition(this->GetPosition());
VoxImage<T> VoxImage<T>::fixVoxelsAroundPlane(const float threshold,
float tolerance, const HPoint3f B,
const HPoint3f E,
bool aboveAir) const {
VoxImage<T> out(*this);
Vector3i dim = this->GetDims();
out.SetDims(dim);
out.SetPosition(this->GetPosition());
Vector3i voxB = this->Find(begin);
Vector3i voxE = this->Find(end);
HPoint3f Bcoll = this->GetPosition().homogeneous();
Vector3i ID;
Vector3i ID;
for (int ix = 0; ix < dim(0); ix++)
for (int iy = 0; iy < dim(1); iy++)
for (int iz = 0; iz < dim(2); iz++) {
ID << ix, iy, iz;
for(int ix=voxB(0); ix<voxE(0); ix++)
for(int iy=voxB(1); iy<voxE(1); iy++)
for(int iz=voxB(2); iz<voxE(2); iz++){
ID << ix,iy,iz;
// voxels around threshold
if( fabs(this->GetValue(ID) - threshold*1.E-6) < tolerance*1.E-6 ){
out.SetValue(ID,-1.*this->GetValue(ID));
}
}
// B, E voxel position
Vector3i iv(ix, iy, iz);
Vector3f v =
Vector3f(iv.cast<float>().cwiseProduct(this->GetSpacing()));
HPoint3f Bvox = Bcoll + HPoint3f(v);
HPoint3f Evox = Bvox + this->GetSpacing().homogeneous();
HPoint3f V = Bvox + 0.5 * (this->GetSpacing().homogeneous());
return out;
// if distance point (x0,y0) from line by points (x1,y1) and (x2,y2) is
// less than tolerance
float x1 = B[1];
float y1 = B[2];
float x2 = E[1];
float y2 = E[2];
float x0 = V[1];
float y0 = V[2];
float dist = fabs((x2 - x1) * (y1 - y0) - ((x1 - x0) * (y2 - y1))) /
sqrt((x2 - x1) * (x2 - x1) + ((y2 - y1) * (y2 - y1)));
float distSign = (x2 - x1) * (y1 - y0) - ((x1 - x0) * (y2 - y1));
// set voxel air value
if (dist < tolerance) {
// std::cout << "voxel " << iv << ", line " << dist << ", tolerance "
// << tolerance << std::endl;
out.SetValue(ID, threshold * 1.E-6);
} else
out.SetValue(ID, this->GetValue(ID));
if ((distSign > 0 && aboveAir) || (distSign < 0 && !aboveAir))
out.SetValue(ID, threshold * 1.E-6);
}
return out;
}
template <typename T>
VoxImage<T> VoxImage<T>::fixVoxels(const HPoint3f begin, const HPoint3f end) const
{
VoxImage<T> out(*this);
out.SetDims(this->GetDims());
out.SetPosition(this->GetPosition());
Vector3i voxB = this->Find(begin);
Vector3i voxE = this->Find(end);
Vector3i ID;
for(int ix=voxB(0); ix<voxE(0); ix++)
for(int iy=voxB(1); iy<voxE(1); iy++)
for(int iz=voxB(2); iz<voxE(2); iz++){
ID << ix,iy,iz;
// voxels around threshold
out.SetValue(ID,-1.*this->GetValue(ID));
}
return out;
template <typename T> void VoxImage<T>::InitVoxels(T t) {
std::fill(m_Data.begin(), m_Data.end(), t); // warning... stl function //
}
template <typename T>
VoxImage<T> VoxImage<T>::fixVoxelsAroundPlane( const float threshold, float tolerance, const HPoint3f B, const HPoint3f E, bool aboveAir) const
{
VoxImage<T> out(*this);
Vector3i dim = this->GetDims();
out.SetDims(dim);
out.SetPosition(this->GetPosition());
HPoint3f Bcoll = this->GetPosition().homogeneous();
Vector3i ID;
for(int ix=0; ix<dim(0); ix++)
for(int iy=0; iy<dim(1); iy++)
for(int iz=0; iz<dim(2); iz++){
ID << ix,iy,iz;
// B, E voxel position
Vector3i iv(ix,iy,iz);
Vector3f v = Vector3f(iv.cast<float>().cwiseProduct(this->GetSpacing()));
HPoint3f Bvox = Bcoll + HPoint3f(v);
HPoint3f Evox = Bvox + this->GetSpacing().homogeneous();
HPoint3f V = Bvox + 0.5*(this->GetSpacing().homogeneous());
// if distance point (x0,y0) from line by points (x1,y1) and (x2,y2) is less than tolerance
float x1 = B[1];
float y1 = B[2];
float x2 = E[1];
float y2 = E[2];
float x0 = V[1];
float y0 = V[2];
float dist = fabs( (x2-x1)*(y1-y0) - ((x1-x0)*(y2-y1))) / sqrt( (x2-x1)*(x2-x1)+((y2-y1)*(y2-y1)));
float distSign = (x2-x1)*(y1-y0) - ((x1-x0)*(y2-y1));
// set voxel air value
if(dist < tolerance){
//std::cout << "voxel " << iv << ", line " << dist << ", tolerance " << tolerance << std::endl;
out.SetValue(ID,threshold*1.E-6);
}
else
out.SetValue(ID,this->GetValue(ID));
if((distSign>0 && aboveAir) || (distSign<0 && !aboveAir) )
out.SetValue(ID,threshold*1.E-6);
}
return out;
}
template<typename T>
void VoxImage<T>::InitVoxels(T t)
{
std::fill( m_Data.begin(), m_Data.end(), t ); // warning... stl function //
}
}
} // namespace uLib
#endif // VOXIMAGE_H

89
src/Math/VoxImageFilter.h
View File

@@ -23,8 +23,6 @@
//////////////////////////////////////////////////////////////////////////////*/
#ifndef VOXIMAGEFILTER_H
#define VOXIMAGEFILTER_H
@@ -33,96 +31,83 @@
#include "Math/VoxImage.h"
namespace uLib {
namespace Interface {
struct VoxImageFilterShape {
template <class Self> void check_structural() {
uLibCheckFunction(Self,operator(),float,float);
uLibCheckFunction(Self,operator(),float,const Vector3f&);
}
template <class Self> void check_structural() {
uLibCheckFunction(Self, operator(), float, float);
uLibCheckFunction(Self, operator(), float, const Vector3f &);
}
};
}
template < typename VoxelT > class Kernel;
} // namespace Interface
template <typename VoxelT> class Kernel;
namespace Abstract {
class VoxImageFilter {
public:
virtual void Run() = 0;
virtual void Run() = 0;
virtual void SetImage(Abstract::VoxImage *image) = 0;
virtual void SetImage(Abstract::VoxImage *image) = 0;
protected:
virtual ~VoxImageFilter() {}
virtual ~VoxImageFilter() {}
};
}
} // namespace Abstract
template < typename VoxelT, typename AlgorithmT >
class VoxImageFilter : public Abstract::VoxImageFilter
{
template <typename VoxelT, typename AlgorithmT>
class VoxImageFilter : public Abstract::VoxImageFilter {
public:
VoxImageFilter(const Vector3i &size);
VoxImageFilter(const Vector3i &size);
void Run();
void Run();
void SetKernelNumericXZY(const std::vector<float> &numeric);
void SetKernelNumericXZY(const std::vector<float> &numeric);
void SetKernelSpherical(float (*shape)(float));
void SetKernelSpherical(float (*shape)(float));
template < class ShapeT >
void SetKernelSpherical( ShapeT shape );
template <class ShapeT> void SetKernelSpherical(ShapeT shape);
void SetKernelWeightFunction(float (*shape)(const Vector3f &));
void SetKernelWeightFunction(float (*shape)(const Vector3f &));
template < class ShapeT >
void SetKernelWeightFunction( ShapeT shape );
template <class ShapeT> void SetKernelWeightFunction(ShapeT shape);
inline Kernel<VoxelT> GetKernelData() const { return this->m_KernelData; }
inline const Kernel<VoxelT> &GetKernelData() const {
return this->m_KernelData;
}
inline Kernel<VoxelT> &GetKernelData() { return this->m_KernelData; }
inline VoxImage<VoxelT>* GetImage() const { return this->m_Image; }
inline VoxImage<VoxelT> *GetImage() const { return this->m_Image; }
void SetImage(Abstract::VoxImage *image);
void SetImage(Abstract::VoxImage *image);
protected:
float Convolve(const VoxImage<VoxelT> &buffer, int index); // remove //
float Convolve(const VoxImage<VoxelT> &buffer, int index); // remove //
void SetKernelOffset();
void SetKernelOffset();
float Distance2(const Vector3i &v);
float Distance2(const Vector3i &v);
// protected members for algorithm access //
Kernel<VoxelT> m_KernelData;
VoxImage<VoxelT> *m_Image;
// protected members for algorithm access //
Kernel<VoxelT> m_KernelData;
VoxImage<VoxelT> *m_Image;
private:
AlgorithmT *t_Algoritm;
AlgorithmT *t_Algoritm;
};
}
} // namespace uLib
#endif // VOXIMAGEFILTER_H
#include "VoxImageFilter.hpp"
#include "VoxImageFilterLinear.hpp"
#include "VoxImageFilterThreshold.hpp"
#include "VoxImageFilterMedian.hpp"
#include "VoxImageFilter2ndStat.hpp"
#include "VoxImageFilterABTrim.hpp"
#include "VoxImageFilterBilateral.hpp"
#include "VoxImageFilter2ndStat.hpp"
#include "VoxImageFilterCustom.hpp"
#include "VoxImageFilterLinear.hpp"
#include "VoxImageFilterMedian.hpp"
#include "VoxImageFilterThreshold.hpp"

368
src/Math/VoxImageFilter.hpp
View File

@@ -23,280 +23,238 @@
//////////////////////////////////////////////////////////////////////////////*/
#ifndef VOXIMAGEFILTER_HPP
#define VOXIMAGEFILTER_HPP
#include <Math/Dense.h>
#include "Math/StructuredData.h"
#include "Math/VoxImage.h"
#include "VoxImageFilter.h"
#include <Math/Dense.h>
namespace uLib {
// KERNEL //////////////////////////////////////////////////////////////////////
template < typename T >
class Kernel : public StructuredData {
typedef StructuredData BaseClass;
template <typename T> class Kernel : public StructuredData {
typedef StructuredData BaseClass;
public:
Kernel(const Vector3i &size);
Kernel(const Vector3i &size);
inline T& operator[](const Vector3i &id) { return m_Data[Map(id)]; }
inline T& operator[](const int &id) { return m_Data[id]; }
inline int GetCenterData() const;
inline T &operator[](const Vector3i &id) { return m_Data[Map(id)]; }
inline T &operator[](const int &id) { return m_Data[id]; }
inline int GetCenterData() const;
inline std::vector<T> & Data() { return this->m_Data; }
inline DataAllocator<T> &Data() { return this->m_Data; }
inline const std::vector<T>& ConstData() const { return this->m_Data; }
inline const DataAllocator<T> &ConstData() const { return this->m_Data; }
void PrintSelf(std::ostream &o) const;
void PrintSelf(std::ostream &o) const;
private:
std::vector<T> m_Data;
DataAllocator<T> m_Data;
};
template < typename T >
Kernel<T>::Kernel(const Vector3i &size) :
BaseClass(size),
m_Data(size.prod())
{
Interface::IsA<T,Interface::Voxel>();
template <typename T>
Kernel<T>::Kernel(const Vector3i &size) : BaseClass(size), m_Data(size.prod()) {
Interface::IsA<T, Interface::Voxel>();
}
template < typename T >
inline int Kernel<T>::GetCenterData() const
{
static int center = Map(this->GetDims() / 2);
return center;
template <typename T> inline int Kernel<T>::GetCenterData() const {
static int center = Map(this->GetDims() / 2);
return center;
}
template < typename T >
void Kernel<T>::PrintSelf(std::ostream &o) const
{
o << " Filter Kernel Dump [XZ_Y]: \n";
Vector3i index;
o << "\n Value: \n\n"
<< "------------------------------------------------- \n";
for (int y = 0 ; y < this->GetDims()(1); ++y ) {
o << "[y=" << y << "]\n";
for (int z = 0 ; z < this->GetDims()(2); ++z ) {
for (int x = 0 ; x < this->GetDims()(0); ++x ) {
index << x,y,z;
o << m_Data[Map(index)].Value << " ";
} o << "\n";
} o << " --------------------------------------------------- \n";
template <typename T> void Kernel<T>::PrintSelf(std::ostream &o) const {
o << " Filter Kernel Dump [XZ_Y]: \n";
Vector3i index;
o << "\n Value: \n\n"
<< "------------------------------------------------- \n";
for (int y = 0; y < this->GetDims()(1); ++y) {
o << "[y=" << y << "]\n";
for (int z = 0; z < this->GetDims()(2); ++z) {
for (int x = 0; x < this->GetDims()(0); ++x) {
index << x, y, z;
o << m_Data[Map(index)].Value << " ";
}
o << "\n";
}
o << "\n Offset: \n"
<< "------------------------------------------------- \n";
for (int y = 0 ; y < this->GetDims()(1); ++y ) {
o << "[y=" << y << "]\n";
for (int z = 0 ; z < this->GetDims()(2); ++z ) {
for (int x = 0 ; x < this->GetDims()(0); ++x ) {
index << x,y,z;
o << m_Data[Map(index)].Count << " ";
} o << "\n";
} o << " --------------------------------------------------- \n";
o << " --------------------------------------------------- \n";
}
o << "\n Offset: \n"
<< "------------------------------------------------- \n";
for (int y = 0; y < this->GetDims()(1); ++y) {
o << "[y=" << y << "]\n";
for (int z = 0; z < this->GetDims()(2); ++z) {
for (int x = 0; x < this->GetDims()(0); ++x) {
index << x, y, z;
o << m_Data[Map(index)].Count << " ";
}
o << "\n";
}
o << " --------------------------------------------------- \n";
}
}
////////////////////////////////////////////////////////////////////////////////
#define _TPL_ template < typename VoxelT , typename AlgorithmT >
#define _TPLT_ VoxelT,AlgorithmT
#define _TPL_ template <typename VoxelT, typename AlgorithmT>
#define _TPLT_ VoxelT, AlgorithmT
_TPL_
VoxImageFilter<_TPLT_>::VoxImageFilter(const Vector3i &size) :
m_KernelData(size),
t_Algoritm(static_cast<AlgorithmT *>(this))
{
VoxImageFilter<_TPLT_>::VoxImageFilter(const Vector3i &size)
: m_KernelData(size), t_Algoritm(static_cast<AlgorithmT *>(this)) {}
_TPL_
void VoxImageFilter<_TPLT_>::Run() {
VoxImage<VoxelT> buffer = *m_Image;
#pragma omp parallel for
for (int i = 0; i < m_Image->Data().size(); ++i)
m_Image->operator[](i).Value = this->t_Algoritm->Evaluate(buffer, i);
#pragma omp barrier
}
_TPL_
void VoxImageFilter<_TPLT_>::Run()
{
VoxImage<VoxelT> buffer = *m_Image;
#pragma omp parallel for
for(int i=0 ; i < m_Image->Data().size() ; ++i)
m_Image->operator [](i).Value = this->t_Algoritm->Evaluate(buffer,i);
#pragma omp barrier
}
_TPL_
void VoxImageFilter<_TPLT_>::SetKernelOffset()
{
Vector3i id(0,0,0);
for( int z=0 ; z < m_KernelData.GetDims()(2); ++z ) {
for( int x=0 ; x < m_KernelData.GetDims()(0); ++x ) {
for( int y=0 ; y < m_KernelData.GetDims()(1); ++y ) {
id << x,y,z;
m_KernelData[id].Count = id.transpose() * m_Image->GetIncrements();
}
}
void VoxImageFilter<_TPLT_>::SetKernelOffset() {
Vector3i id(0, 0, 0);
for (int z = 0; z < m_KernelData.GetDims()(2); ++z) {
for (int x = 0; x < m_KernelData.GetDims()(0); ++x) {
for (int y = 0; y < m_KernelData.GetDims()(1); ++y) {
id << x, y, z;
m_KernelData[id].Count = id.transpose() * m_Image->GetIncrements();
}
}
}
}
_TPL_
float VoxImageFilter<_TPLT_>::Distance2(const Vector3i &v)
{
Vector3i tmp = v;
const Vector3i &dim = this->m_KernelData.GetDims();
Vector3i center = dim / 2;
tmp = tmp - center;
center = center.cwiseProduct(center);
tmp = tmp.cwiseProduct(tmp);
return (float)(tmp.sum()) / (float)( center.sum() + 0.25 *
(3 - (dim(0) % 2) - (dim(1) % 2) - (dim(2) % 2)));
float VoxImageFilter<_TPLT_>::Distance2(const Vector3i &v) {
Vector3i tmp = v;
const Vector3i &dim = this->m_KernelData.GetDims();
Vector3i center = dim / 2;
tmp = tmp - center;
center = center.cwiseProduct(center);
tmp = tmp.cwiseProduct(tmp);
return (float)(tmp.sum()) /
(float)(center.sum() +
0.25 * (3 - (dim(0) % 2) - (dim(1) % 2) - (dim(2) % 2)));
}
_TPL_
void VoxImageFilter<_TPLT_>::SetKernelNumericXZY(const std::vector<float> &numeric)
{
// set data order //
StructuredData::Order order = m_KernelData.GetDataOrder();
//m_KernelData.SetDataOrder(StructuredData::XZY);
Vector3i id;
int index = 0;
for( int y=0 ; y < m_KernelData.GetDims()(1); ++y ) {
for( int z=0 ; z < m_KernelData.GetDims()(2); ++z ) {
for( int x=0 ; x < m_KernelData.GetDims()(0); ++x ) {
id << x,y,z;
m_KernelData[id].Value = numeric[index++];
}
}
void VoxImageFilter<_TPLT_>::SetKernelNumericXZY(
const std::vector<float> &numeric) {
// set data order //
StructuredData::Order order = m_KernelData.GetDataOrder();
// m_KernelData.SetDataOrder(StructuredData::XZY);
Vector3i id;
int index = 0;
for (int y = 0; y < m_KernelData.GetDims()(1); ++y) {
for (int z = 0; z < m_KernelData.GetDims()(2); ++z) {
for (int x = 0; x < m_KernelData.GetDims()(0); ++x) {
id << x, y, z;
m_KernelData[id].Value = numeric[index++];
}
}
//m_KernelData.SetDataOrder(order);
}
// m_KernelData.SetDataOrder(order);
}
_TPL_
void VoxImageFilter<_TPLT_>::SetKernelSpherical(float(* shape)(float))
{
Vector3i id;
for( int y=0 ; y < m_KernelData.GetDims()(1); ++y ) {
for( int z=0 ; z < m_KernelData.GetDims()(2); ++z ) {
for( int x=0 ; x < m_KernelData.GetDims()(0); ++x ) {
id << x,y,z;
m_KernelData[id].Value = shape(this->Distance2(id));
}
}
void VoxImageFilter<_TPLT_>::SetKernelSpherical(float (*shape)(float)) {
Vector3i id;
for (int y = 0; y < m_KernelData.GetDims()(1); ++y) {
for (int z = 0; z < m_KernelData.GetDims()(2); ++z) {
for (int x = 0; x < m_KernelData.GetDims()(0); ++x) {
id << x, y, z;
m_KernelData[id].Value = shape(this->Distance2(id));
}
}
}
}
_TPL_ template <class ShapeT>
void VoxImageFilter<_TPLT_>::SetKernelSpherical(ShapeT shape)
{
Interface::IsA<ShapeT,Interface::VoxImageFilterShape>();
Vector3i id;
for( int y=0 ; y < m_KernelData.GetDims()(1); ++y ) {
for( int z=0 ; z < m_KernelData.GetDims()(2); ++z ) {
for( int x=0 ; x < m_KernelData.GetDims()(0); ++x ) {
id << x,y,z;
m_KernelData[id].Value = shape(this->Distance2(id));
}
}
void VoxImageFilter<_TPLT_>::SetKernelSpherical(ShapeT shape) {
Interface::IsA<ShapeT, Interface::VoxImageFilterShape>();
Vector3i id;
for (int y = 0; y < m_KernelData.GetDims()(1); ++y) {
for (int z = 0; z < m_KernelData.GetDims()(2); ++z) {
for (int x = 0; x < m_KernelData.GetDims()(0); ++x) {
id << x, y, z;
m_KernelData[id].Value = shape(this->Distance2(id));
}
}
}
}
_TPL_
void VoxImageFilter<_TPLT_>::SetKernelWeightFunction(float (*shape)(const Vector3f &))
{
const Vector3i &dim = m_KernelData.GetDims();
Vector3i id;
Vector3f pt;
for( int y=0 ; y < dim(1); ++y ) {
for( int z=0 ; z < dim(2); ++z ) {
for( int x=0 ; x < dim(0); ++x ) {
// get voxels centroid coords from kernel center //
id << x,y,z;
pt << id(0) - dim(0)/2 + 0.5 * !(dim(0) % 2),
id(1) - dim(1)/2 + 0.5 * !(dim(1) % 2),
id(2) - dim(2)/2 + 0.5 * !(dim(2) % 2);
// compute function using given shape //
m_KernelData[id].Value = shape(pt);
}
}
void VoxImageFilter<_TPLT_>::SetKernelWeightFunction(
float (*shape)(const Vector3f &)) {
const Vector3i &dim = m_KernelData.GetDims();
Vector3i id;
Vector3f pt;
for (int y = 0; y < dim(1); ++y) {
for (int z = 0; z < dim(2); ++z) {
for (int x = 0; x < dim(0); ++x) {
// get voxels centroid coords from kernel center //
id << x, y, z;
pt << id(0) - dim(0) / 2 + 0.5 * !(dim(0) % 2),
id(1) - dim(1) / 2 + 0.5 * !(dim(1) % 2),
id(2) - dim(2) / 2 + 0.5 * !(dim(2) % 2);
// compute function using given shape //
m_KernelData[id].Value = shape(pt);
}
}
}
}
_TPL_ template < class ShapeT >
void VoxImageFilter<_TPLT_>::SetKernelWeightFunction(ShapeT shape)
{
Interface::IsA<ShapeT,Interface::VoxImageFilterShape>();
const Vector3i &dim = m_KernelData.GetDims();
Vector3i id;
Vector3f pt;
for( int y=0 ; y < dim(1); ++y ) {
for( int z=0 ; z < dim(2); ++z ) {
for( int x=0 ; x < dim(0); ++x ) {
// get voxels centroid coords from kernel center //
id << x,y,z;
pt << id(0) - dim(0)/2 + 0.5 * !(dim(0) % 2),
id(1) - dim(1)/2 + 0.5 * !(dim(1) % 2),
id(2) - dim(2)/2 + 0.5 * !(dim(2) % 2);
// compute function using given shape //
m_KernelData[id].Value = shape(pt);
}
}
_TPL_ template <class ShapeT>
void VoxImageFilter<_TPLT_>::SetKernelWeightFunction(ShapeT shape) {
Interface::IsA<ShapeT, Interface::VoxImageFilterShape>();
const Vector3i &dim = m_KernelData.GetDims();
Vector3i id;
Vector3f pt;
for (int y = 0; y < dim(1); ++y) {
for (int z = 0; z < dim(2); ++z) {
for (int x = 0; x < dim(0); ++x) {
// get voxels centroid coords from kernel center //
id << x, y, z;
pt << id(0) - dim(0) / 2 + 0.5 * !(dim(0) % 2),
id(1) - dim(1) / 2 + 0.5 * !(dim(1) % 2),
id(2) - dim(2) / 2 + 0.5 * !(dim(2) % 2);
// compute function using given shape //
m_KernelData[id].Value = shape(pt);
}
}
}
}
_TPL_
void VoxImageFilter<_TPLT_>::SetImage(Abstract::VoxImage *image)
{
this->m_Image = reinterpret_cast<VoxImage<VoxelT> *> (image);
this->SetKernelOffset();
void VoxImageFilter<_TPLT_>::SetImage(Abstract::VoxImage *image) {
this->m_Image = reinterpret_cast<VoxImage<VoxelT> *>(image);
this->SetKernelOffset();
}
_TPL_
float VoxImageFilter<_TPLT_>::Convolve(const VoxImage<VoxelT> &buffer, int index)
{
const std::vector<VoxelT> &vbuf = buffer.ConstData();
const std::vector<VoxelT> &vker = m_KernelData.ConstData();
int vox_size = vbuf.size();
int ker_size = vker.size();
int pos;
float conv = 0, ksum = 0;
for (int ik = 0; ik < ker_size; ++ik) {
pos = index + vker[ik].Count - vker[m_KernelData.GetCenterData()].Count;
pos = (pos + vox_size) % vox_size;
conv += vbuf[pos].Value * vker[ik].Value;
ksum += vker[ik].Value;
}
return conv / ksum;
float VoxImageFilter<_TPLT_>::Convolve(const VoxImage<VoxelT> &buffer,
int index) {
const DataAllocator<VoxelT> &vbuf = buffer.ConstData();
const DataAllocator<VoxelT> &vker = m_KernelData.ConstData();
int vox_size = vbuf.size();
int ker_size = vker.size();
int pos;
float conv = 0, ksum = 0;
for (int ik = 0; ik < ker_size; ++ik) {
pos = index + vker[ik].Count - vker[m_KernelData.GetCenterData()].Count;
pos = (pos + vox_size) % vox_size;
conv += vbuf[pos].Value * vker[ik].Value;
ksum += vker[ik].Value;
}
return conv / ksum;
}
#undef _TPLT_
#undef _TPL_
}
} // namespace uLib
#endif // VOXIMAGEFILTER_HPP

67
src/Math/VoxImageFilter2ndStat.hpp
View File

@@ -23,14 +23,12 @@
//////////////////////////////////////////////////////////////////////////////*/
#ifndef VOXIMAGEFILTER2NDSTAT_HPP
#define VOXIMAGEFILTER2NDSTAT_HPP
#include <Math/Dense.h>
#include "Math/VoxImage.h"
#include "VoxImageFilter.h"
#include <Math/Dense.h>
////////////////////////////////////////////////////////////////////////////////
///// VOXIMAGE FILTER ABTRIM /////////////////////////////////////////////////
@@ -39,45 +37,42 @@
namespace uLib {
template <typename VoxelT>
class VoxFilterAlgorithm2ndStat :
public VoxImageFilter<VoxelT, VoxFilterAlgorithm2ndStat<VoxelT> > {
class VoxFilterAlgorithm2ndStat
: public VoxImageFilter<VoxelT, VoxFilterAlgorithm2ndStat<VoxelT>> {
public:
typedef VoxImageFilter<VoxelT, VoxFilterAlgorithm2ndStat<VoxelT> > BaseClass;
VoxFilterAlgorithm2ndStat(const Vector3i &size) :
BaseClass(size)
{ }
typedef VoxImageFilter<VoxelT, VoxFilterAlgorithm2ndStat<VoxelT>> BaseClass;
VoxFilterAlgorithm2ndStat(const Vector3i &size) : BaseClass(size) {}
float Evaluate(const VoxImage<VoxelT> &buffer, int index)
{
const std::vector<VoxelT> &vbuf = buffer.ConstData();
const std::vector<VoxelT> &vker = this->m_KernelData.ConstData();
int vox_size = vbuf.size();
int ker_size = vker.size();
int pos;
float Evaluate(const VoxImage<VoxelT> &buffer, int index) {
const DataAllocator<VoxelT> &vbuf = buffer.ConstData();
const DataAllocator<VoxelT> &vker = this->m_KernelData.ConstData();
int vox_size = vbuf.size();
int ker_size = vker.size();
int pos;
// mean //
float conv = 0, ksum = 0;
for (int ik = 0; ik < ker_size; ++ik) {
pos = index + vker[ik].Count - vker[this->m_KernelData.GetCenterData()].Count;
pos = (pos + vox_size) % vox_size;
conv += vbuf[pos].Value * vker[ik].Value;
ksum += vker[ik].Value;
}
float mean = conv / ksum;
// rms //
conv = 0;
for (int ik = 0; ik < ker_size; ++ik) {
pos = index + vker[ik].Count - vker[this->m_KernelData.GetCenterData()].Count;
pos = (pos + vox_size) % vox_size;
conv += pow((vbuf[pos].Value * vker[ik].Value) - mean , 2);
}
return conv / (vker.size() - 1) ;
// mean //
float conv = 0, ksum = 0;
for (int ik = 0; ik < ker_size; ++ik) {
pos = index + vker[ik].Count -
vker[this->m_KernelData.GetCenterData()].Count;
pos = (pos + vox_size) % vox_size;
conv += vbuf[pos].Value * vker[ik].Value;
ksum += vker[ik].Value;
}
float mean = conv / ksum;
// rms //
conv = 0;
for (int ik = 0; ik < ker_size; ++ik) {
pos = index + vker[ik].Count -
vker[this->m_KernelData.GetCenterData()].Count;
pos = (pos + vox_size) % vox_size;
conv += pow((vbuf[pos].Value * vker[ik].Value) - mean, 2);
}
return conv / (vker.size() - 1);
}
};
}
} // namespace uLib
#endif // VOXIMAGEFILTER2NDSTAT_HPP

335
src/Math/VoxImageFilterABTrim.hpp
View File

@@ -23,14 +23,12 @@
//////////////////////////////////////////////////////////////////////////////*/
#ifndef VOXIMAGEFILTERABTRIM_HPP
#define VOXIMAGEFILTERABTRIM_HPP
#include <Math/Dense.h>
#include "Math/VoxImage.h"
#include "VoxImageFilter.h"
#include <Math/Dense.h>
////////////////////////////////////////////////////////////////////////////////
///// VOXIMAGE FILTER ABTRIM /////////////////////////////////////////////////
@@ -38,142 +36,257 @@
namespace uLib {
#ifdef USE_CUDA
template <typename VoxelT>
class VoxFilterAlgorithmAbtrim :
public VoxImageFilter<VoxelT, VoxFilterAlgorithmAbtrim<VoxelT> > {
__global__ void ABTrimFilterKernel(const VoxelT *in, VoxelT *out,
const VoxelT *kernel, int vox_size,
int ker_size, int center_count, int mAtrim,
int mBtrim) {
int index = blockIdx.x * blockDim.x + threadIdx.x;
if (index < vox_size) {
// Allocate space for sorting
extern __shared__ char shared_mem[];
VoxelT *mfh =
(VoxelT *)&shared_mem[threadIdx.x * ker_size * sizeof(VoxelT)];
struct KernelSortAscending
{
bool operator()(const VoxelT& e1, const VoxelT& e2)
{ return e1.Value < e2.Value; }
};
for (int i = 0; i < ker_size; ++i) {
mfh[i].Count = i;
}
for (int ik = 0; ik < ker_size; ik++) {
int pos = index + kernel[ik].Count - center_count;
if (pos < 0) {
pos += vox_size * ((-pos / vox_size) + 1);
}
pos = pos % vox_size;
mfh[ik].Value = in[pos].Value;
}
// Simple bubble sort for small arrays
for (int i = 0; i < ker_size - 1; i++) {
for (int j = 0; j < ker_size - i - 1; j++) {
if (mfh[j].Value > mfh[j + 1].Value) {
VoxelT temp = mfh[j];
mfh[j] = mfh[j + 1];
mfh[j + 1] = temp;
}
}
}
float ker_sum = 0;
float fconv = 0;
for (int ik = 0; ik < mAtrim; ik++) {
ker_sum += kernel[mfh[ik].Count].Value;
}
for (int ik = mAtrim; ik < ker_size - mBtrim; ik++) {
fconv += mfh[ik].Value * kernel[mfh[ik].Count].Value;
ker_sum += kernel[mfh[ik].Count].Value;
}
for (int ik = ker_size - mBtrim; ik < ker_size; ik++) {
ker_sum += kernel[mfh[ik].Count].Value;
}
out[index].Value = fconv / ker_sum;
}
}
#endif
template <typename VoxelT>
class VoxFilterAlgorithmAbtrim
: public VoxImageFilter<VoxelT, VoxFilterAlgorithmAbtrim<VoxelT>> {
struct KernelSortAscending {
bool operator()(const VoxelT &e1, const VoxelT &e2) {
return e1.Value < e2.Value;
}
};
public:
typedef VoxImageFilter<VoxelT, VoxFilterAlgorithmAbtrim<VoxelT> > BaseClass;
VoxFilterAlgorithmAbtrim(const Vector3i &size) :
BaseClass(size)
{
mAtrim = 0;
mBtrim = 0;
typedef VoxImageFilter<VoxelT, VoxFilterAlgorithmAbtrim<VoxelT>> BaseClass;
VoxFilterAlgorithmAbtrim(const Vector3i &size) : BaseClass(size) {
mAtrim = 0;
mBtrim = 0;
}
#ifdef USE_CUDA
void Run() {
if (this->m_Image->Data().GetDevice() == MemoryDevice::VRAM ||
this->m_KernelData.Data().GetDevice() == MemoryDevice::VRAM) {
this->m_Image->Data().MoveToVRAM();
this->m_KernelData.Data().MoveToVRAM();
VoxImage<VoxelT> buffer = *(this->m_Image);
buffer.Data().MoveToVRAM();
int vox_size = buffer.Data().size();
int ker_size = this->m_KernelData.Data().size();
VoxelT *d_img_out = this->m_Image->Data().GetVRAMData();
const VoxelT *d_img_in = buffer.Data().GetVRAMData();
const VoxelT *d_kernel = this->m_KernelData.Data().GetVRAMData();
int center_count =
this->m_KernelData[this->m_KernelData.GetCenterData()].Count;
int threadsPerBlock = 256;
int blocksPerGrid = (vox_size + threadsPerBlock - 1) / threadsPerBlock;
size_t shared_mem_size = threadsPerBlock * ker_size * sizeof(VoxelT);
ABTrimFilterKernel<<<blocksPerGrid, threadsPerBlock, shared_mem_size>>>(
d_img_in, d_img_out, d_kernel, vox_size, ker_size, center_count,
mAtrim, mBtrim);
cudaDeviceSynchronize();
} else {
BaseClass::Run();
}
}
#endif
float Evaluate(const VoxImage<VoxelT> &buffer, int index) {
const DataAllocator<VoxelT> &vbuf = buffer.ConstData();
const DataAllocator<VoxelT> &vker = this->m_KernelData.ConstData();
int vox_size = vbuf.size();
int ker_size = vker.size();
int pos;
std::vector<VoxelT> mfh(ker_size);
for (int i = 0; i < ker_size; ++i)
mfh[i].Count = i; // index key for ordering function
for (int ik = 0; ik < ker_size; ik++) {
pos = index + vker[ik].Count -
vker[this->m_KernelData.GetCenterData()].Count;
pos = (pos + vox_size) % vox_size;
mfh[ik].Value = vbuf[pos].Value;
}
float Evaluate(const VoxImage<VoxelT> &buffer, int index)
{
const std::vector<VoxelT> &vbuf = buffer.ConstData();
const std::vector<VoxelT> &vker = this->m_KernelData.ConstData();
int vox_size = vbuf.size();
int ker_size = vker.size();
int pos;
std::vector<VoxelT> mfh(ker_size);
for (int i = 0; i < ker_size; ++i)
mfh[i].Count = i; //index key for ordering function
for (int ik = 0; ik < ker_size; ik++) {
pos = index + vker[ik].Count - vker[this->m_KernelData.GetCenterData()].Count;
pos = (pos + vox_size) % vox_size;
mfh[ik].Value = vbuf[pos].Value;
}
std::sort(mfh.begin(), mfh.end(), KernelSortAscending());
float ker_sum = 0;
float fconv = 0;
for (int ik = 0; ik < mAtrim; ik++)
ker_sum += vker[ mfh[ik].Count ].Value;
for (int ik = mAtrim; ik < ker_size - mBtrim; ik++) {
fconv += mfh[ik].Value * vker[ mfh[ik].Count ].Value; // convloution //
ker_sum += vker[ mfh[ik].Count ].Value;
}
for (int ik = ker_size - mBtrim; ik < ker_size; ik++)
ker_sum += vker[ mfh[ik].Count ].Value;
return fconv / ker_sum;
std::sort(mfh.begin(), mfh.end(), KernelSortAscending());
float ker_sum = 0;
float fconv = 0;
for (int ik = 0; ik < mAtrim; ik++)
ker_sum += vker[mfh[ik].Count].Value;
for (int ik = mAtrim; ik < ker_size - mBtrim; ik++) {
fconv += mfh[ik].Value * vker[mfh[ik].Count].Value; // convloution //
ker_sum += vker[mfh[ik].Count].Value;
}
for (int ik = ker_size - mBtrim; ik < ker_size; ik++)
ker_sum += vker[mfh[ik].Count].Value;
inline void SetABTrim(int a, int b) { mAtrim = a; mBtrim = b; }
return fconv / ker_sum;
}
inline void SetABTrim(int a, int b) {
mAtrim = a;
mBtrim = b;
}
private:
int mAtrim;
int mBtrim;
int mAtrim;
int mBtrim;
};
////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////
// Roberspierre Filter //
template <typename VoxelT>
class VoxFilterAlgorithmSPR :
public VoxImageFilter<VoxelT, VoxFilterAlgorithmSPR<VoxelT> > {
class VoxFilterAlgorithmSPR
: public VoxImageFilter<VoxelT, VoxFilterAlgorithmSPR<VoxelT>> {
struct KernelSortAscending
{
bool operator()(const VoxelT& e1, const VoxelT& e2)
{ return e1.Value < e2.Value; }
};
struct KernelSortAscending {
bool operator()(const VoxelT &e1, const VoxelT &e2) {
return e1.Value < e2.Value;
}
};
public:
typedef VoxImageFilter<VoxelT, VoxFilterAlgorithmSPR<VoxelT> > BaseClass;
VoxFilterAlgorithmSPR(const Vector3i &size) :
BaseClass(size)
{
mAtrim = 0;
mBtrim = 0;
typedef VoxImageFilter<VoxelT, VoxFilterAlgorithmSPR<VoxelT>> BaseClass;
VoxFilterAlgorithmSPR(const Vector3i &size) : BaseClass(size) {
mAtrim = 0;
mBtrim = 0;
}
#ifdef USE_CUDA
void Run() {
if (this->m_Image->Data().GetDevice() == MemoryDevice::VRAM ||
this->m_KernelData.Data().GetDevice() == MemoryDevice::VRAM) {
this->m_Image->Data().MoveToVRAM();
this->m_KernelData.Data().MoveToVRAM();
VoxImage<VoxelT> buffer = *(this->m_Image);
buffer.Data().MoveToVRAM();
int vox_size = buffer.Data().size();
int ker_size = this->m_KernelData.Data().size();
VoxelT *d_img_out = this->m_Image->Data().GetVRAMData();
const VoxelT *d_img_in = buffer.Data().GetVRAMData();
const VoxelT *d_kernel = this->m_KernelData.Data().GetVRAMData();
int center_count =
this->m_KernelData[this->m_KernelData.GetCenterData()].Count;
int threadsPerBlock = 256;
int blocksPerGrid = (vox_size + threadsPerBlock - 1) / threadsPerBlock;
size_t shared_mem_size = threadsPerBlock * ker_size * sizeof(VoxelT);
ABTrimFilterKernel<<<blocksPerGrid, threadsPerBlock, shared_mem_size>>>(
d_img_in, d_img_out, d_kernel, vox_size, ker_size, center_count,
mAtrim, mBtrim);
cudaDeviceSynchronize();
} else {
BaseClass::Run();
}
}
#endif
float Evaluate(const VoxImage<VoxelT> &buffer, int index) {
const DataAllocator<VoxelT> &vbuf = buffer.ConstData();
const DataAllocator<VoxelT> &vker = this->m_KernelData.ConstData();
int vox_size = vbuf.size();
int ker_size = vker.size();
int pos;
std::vector<VoxelT> mfh(ker_size);
for (int i = 0; i < ker_size; ++i)
mfh[i].Count = i; // index key for ordering function
for (int ik = 0; ik < ker_size; ik++) {
pos = index + vker[ik].Count -
vker[this->m_KernelData.GetCenterData()].Count;
pos = (pos + vox_size) % vox_size;
mfh[ik].Value = vbuf[pos].Value;
}
float Evaluate(const VoxImage<VoxelT> &buffer, int index)
{
const std::vector<VoxelT> &vbuf = buffer.ConstData();
const std::vector<VoxelT> &vker = this->m_KernelData.ConstData();
int vox_size = vbuf.size();
int ker_size = vker.size();
int pos;
std::sort(mfh.begin(), mfh.end(), KernelSortAscending());
float spr = vbuf[index].Value;
if ((mAtrim > 0 && spr <= mfh[mAtrim - 1].Value) ||
(mBtrim > 0 && spr >= mfh[ker_size - mBtrim].Value)) {
float ker_sum = 0;
float fconv = 0;
for (int ik = 0; ik < mAtrim; ik++)
ker_sum += vker[mfh[ik].Count].Value;
for (int ik = mAtrim; ik < ker_size - mBtrim; ik++) {
fconv += mfh[ik].Value * vker[mfh[ik].Count].Value;
ker_sum += vker[mfh[ik].Count].Value;
}
for (int ik = ker_size - mBtrim; ik < ker_size; ik++)
ker_sum += vker[mfh[ik].Count].Value;
std::vector<VoxelT> mfh(ker_size);
for (int i = 0; i < ker_size; ++i)
mfh[i].Count = i; //index key for ordering function
for (int ik = 0; ik < ker_size; ik++) {
pos = index + vker[ik].Count -
vker[this->m_KernelData.GetCenterData()].Count;
pos = (pos + vox_size) % vox_size;
mfh[ik].Value = vbuf[pos].Value;
}
return fconv / ker_sum;
} else
return spr;
}
std::sort(mfh.begin(), mfh.end(), KernelSortAscending());
float spr = vbuf[index].Value;
if( (mAtrim > 0 && spr <= mfh[mAtrim-1].Value) ||
(mBtrim > 0 && spr >= mfh[ker_size - mBtrim].Value) )
{
float ker_sum = 0;
float fconv = 0;
for (int ik = 0; ik < mAtrim; ik++)
ker_sum += vker[ mfh[ik].Count ].Value;
for (int ik = mAtrim; ik < ker_size - mBtrim; ik++) {
fconv += mfh[ik].Value * vker[ mfh[ik].Count ].Value;
ker_sum += vker[ mfh[ik].Count ].Value;
}
for (int ik = ker_size - mBtrim; ik < ker_size; ik++)
ker_sum += vker[ mfh[ik].Count ].Value;
return fconv / ker_sum;
}
else
return spr;
}
inline void SetABTrim(int a, int b) { mAtrim = a; mBtrim = b; }
inline void SetABTrim(int a, int b) {
mAtrim = a;
mBtrim = b;
}
private:
int mAtrim;
int mBtrim;
int mAtrim;
int mBtrim;
};
}
} // namespace uLib
#endif // VOXIMAGEFILTERABTRIM_HPP

180
src/Math/VoxImageFilterBilateral.hpp
View File

@@ -23,14 +23,12 @@
//////////////////////////////////////////////////////////////////////////////*/
#ifndef VOXIMAGEFILTERBILATERAL_HPP
#define VOXIMAGEFILTERBILATERAL_HPP
#include <Math/Dense.h>
#include "Math/VoxImage.h"
#include "VoxImageFilter.h"
#include <Math/Dense.h>
////////////////////////////////////////////////////////////////////////////////
///// VOXIMAGE FILTER LINEAR /////////////////////////////////////////////////
@@ -38,115 +36,119 @@
namespace uLib {
template <typename VoxelT>
class VoxFilterAlgorithmBilateral :
public VoxImageFilter<VoxelT, VoxFilterAlgorithmBilateral<VoxelT> > {
class VoxFilterAlgorithmBilateral
: public VoxImageFilter<VoxelT, VoxFilterAlgorithmBilateral<VoxelT>> {
public:
typedef VoxImageFilter<VoxelT, VoxFilterAlgorithmBilateral<VoxelT> > BaseClass;
VoxFilterAlgorithmBilateral(const Vector3i &size) : BaseClass(size) {
m_sigma = 1;
}
typedef VoxImageFilter<VoxelT, VoxFilterAlgorithmBilateral<VoxelT>> BaseClass;
VoxFilterAlgorithmBilateral(const Vector3i &size) : BaseClass(size) {
m_sigma = 1;
}
float Evaluate(const VoxImage<VoxelT> &buffer, int index)
{
const std::vector<VoxelT> &vbuf = buffer.ConstData();
const std::vector<VoxelT> &vker = this->m_KernelData.ConstData();
int vox_size = vbuf.size();
int ker_size = vker.size();
int pos;
float conv = 0, ksum = 0;
float gamma_smooth;
for (int ik = 0; ik < ker_size; ++ik) {
// if (ik==this->m_KernelData.GetCenterData()) continue;
pos = index + vker[ik].Count - vker[this->m_KernelData.GetCenterData()].Count;
pos = (pos + vox_size) % vox_size;
gamma_smooth = compute_gauss( fabs(vbuf[index].Value - vbuf[pos].Value) * 1.E6 );
conv += vbuf[pos].Value * vker[ik].Value * gamma_smooth;
ksum += vker[ik].Value * gamma_smooth;
}
return conv / ksum;
float Evaluate(const VoxImage<VoxelT> &buffer, int index) {
const DataAllocator<VoxelT> &vbuf = buffer.ConstData();
const DataAllocator<VoxelT> &vker = this->m_KernelData.ConstData();
int vox_size = vbuf.size();
int ker_size = vker.size();
int pos;
float conv = 0, ksum = 0;
float gamma_smooth;
for (int ik = 0; ik < ker_size; ++ik) {
// if (ik==this->m_KernelData.GetCenterData()) continue;
pos = index + vker[ik].Count -
vker[this->m_KernelData.GetCenterData()].Count;
pos = (pos + vox_size) % vox_size;
gamma_smooth =
compute_gauss(fabs(vbuf[index].Value - vbuf[pos].Value) * 1.E6);
conv += vbuf[pos].Value * vker[ik].Value * gamma_smooth;
ksum += vker[ik].Value * gamma_smooth;
}
return conv / ksum;
}
inline void SetIntensitySigma(const float s) { m_sigma = s; }
inline void SetIntensitySigma(const float s) { m_sigma = s; }
private:
inline float compute_gauss(const float x) {
return 1/(sqrt(2*M_PI)* m_sigma) * exp(-0.5*(x*x)/(m_sigma*m_sigma));
}
inline float compute_gauss(const float x) {
return 1 / (sqrt(2 * M_PI) * m_sigma) *
exp(-0.5 * (x * x) / (m_sigma * m_sigma));
}
Scalarf m_sigma;
Scalarf m_sigma;
};
template <typename VoxelT>
class VoxFilterAlgorithmBilateralTrim :
public VoxImageFilter<VoxelT, VoxFilterAlgorithmBilateralTrim<VoxelT> > {
class VoxFilterAlgorithmBilateralTrim
: public VoxImageFilter<VoxelT, VoxFilterAlgorithmBilateralTrim<VoxelT>> {
typedef std::pair<float,float> FPair;
typedef std::pair<float, float> FPair;
struct KernelSortAscending
{
bool operator()(const FPair& e1, const FPair& e2)
{ return e1.second < e2.second; }
};
struct KernelSortAscending {
bool operator()(const FPair &e1, const FPair &e2) {
return e1.second < e2.second;
}
};
public:
typedef VoxImageFilter<VoxelT, VoxFilterAlgorithmBilateralTrim<VoxelT> > BaseClass;
VoxFilterAlgorithmBilateralTrim(const Vector3i &size) : BaseClass(size) {
m_sigma = 1;
mAtrim = 0;
mBtrim = 0;
typedef VoxImageFilter<VoxelT, VoxFilterAlgorithmBilateralTrim<VoxelT>>
BaseClass;
VoxFilterAlgorithmBilateralTrim(const Vector3i &size) : BaseClass(size) {
m_sigma = 1;
mAtrim = 0;
mBtrim = 0;
}
float Evaluate(const VoxImage<VoxelT> &buffer, int index) {
const DataAllocator<VoxelT> &vbuf = buffer.ConstData();
const DataAllocator<VoxelT> &vker = this->m_KernelData.ConstData();
int img_size = vbuf.size();
int ker_size = vker.size();
int pos;
std::vector<FPair> mfh(ker_size);
for (int i = 0; i < ker_size; ++i)
mfh[i].first = vker[i].Value; // kernel value in first
for (int ik = 0; ik < ker_size; ik++) {
pos = index + vker[ik].Count -
vker[this->m_KernelData.GetCenterData()].Count;
pos = (pos + img_size) % img_size;
mfh[ik].second = vbuf[pos].Value; // image value in second
}
std::sort(mfh.begin(), mfh.end(), KernelSortAscending());
float Evaluate(const VoxImage<VoxelT> &buffer, int index)
{
const std::vector<VoxelT> &vbuf = buffer.ConstData();
const std::vector<VoxelT> &vker = this->m_KernelData.ConstData();
int img_size = vbuf.size();
int ker_size = vker.size();
int pos;
std::vector<FPair> mfh(ker_size);
for (int i = 0; i < ker_size; ++i)
mfh[i].first = vker[i].Value; // kernel value in first
for (int ik = 0; ik < ker_size; ik++) {
pos = index + vker[ik].Count - vker[this->m_KernelData.GetCenterData()].Count;
pos = (pos + img_size) % img_size;
mfh[ik].second = vbuf[pos].Value; // image value in second
}
std::sort(mfh.begin(), mfh.end(), KernelSortAscending());
float conv = 0, ksum = 0;
float gamma_smooth;
// for (int ik = 0; ik < mAtrim; ik++)
// ksum += mfh[ik].first;
for (int ik = mAtrim; ik < ker_size - mBtrim; ik++) {
gamma_smooth = compute_gauss( fabs(vbuf[index].Value - mfh[ik].second) * 1.E6 );
conv += mfh[ik].first * mfh[ik].second * gamma_smooth;
ksum += mfh[ik].first * gamma_smooth;
}
// for (int ik = ker_size - mBtrim; ik < ker_size; ik++)
// ksum += mfh[ik].first;
return conv / ksum;
float conv = 0, ksum = 0;
float gamma_smooth;
// for (int ik = 0; ik < mAtrim; ik++)
// ksum += mfh[ik].first;
for (int ik = mAtrim; ik < ker_size - mBtrim; ik++) {
gamma_smooth =
compute_gauss(fabs(vbuf[index].Value - mfh[ik].second) * 1.E6);
conv += mfh[ik].first * mfh[ik].second * gamma_smooth;
ksum += mfh[ik].first * gamma_smooth;
}
// for (int ik = ker_size - mBtrim; ik < ker_size; ik++)
// ksum += mfh[ik].first;
inline void SetIntensitySigma(const float s) { m_sigma = s; }
inline void SetABTrim(int a, int b) { mAtrim = a; mBtrim = b; }
return conv / ksum;
}
inline void SetIntensitySigma(const float s) { m_sigma = s; }
inline void SetABTrim(int a, int b) {
mAtrim = a;
mBtrim = b;
}
private:
inline float compute_gauss(const float x) {
return 1/(sqrt(2*M_PI)* m_sigma) * exp(-0.5*(x*x)/(m_sigma*m_sigma));
}
inline float compute_gauss(const float x) {
return 1 / (sqrt(2 * M_PI) * m_sigma) *
exp(-0.5 * (x * x) / (m_sigma * m_sigma));
}
Scalarf m_sigma;
int mAtrim;
int mBtrim;
Scalarf m_sigma;
int mAtrim;
int mBtrim;
};
}
} // namespace uLib
#endif // VOXIMAGEFILTERBILATERAL_HPP

70
src/Math/VoxImageFilterCustom.hpp
View File

@@ -23,14 +23,12 @@
//////////////////////////////////////////////////////////////////////////////*/
#ifndef VOXIMAGEFILTERCUSTOM_HPP
#define VOXIMAGEFILTERCUSTOM_HPP
#include <Math/Dense.h>
#include "Math/VoxImage.h"
#include "VoxImageFilter.h"
#include <Math/Dense.h>
#define likely(expr) __builtin_expect(!!(expr), 1)
@@ -41,50 +39,50 @@
namespace uLib {
template <typename VoxelT>
class VoxFilterAlgorithmCustom :
public VoxImageFilter<VoxelT, VoxFilterAlgorithmCustom<VoxelT> > {
class VoxFilterAlgorithmCustom
: public VoxImageFilter<VoxelT, VoxFilterAlgorithmCustom<VoxelT>> {
typedef float (*FunctionPt)(const std::vector<Scalarf> &);
typedef float (* FunctionPt)(const std::vector<Scalarf> &);
public:
typedef VoxImageFilter<VoxelT, VoxFilterAlgorithmCustom<VoxelT> > BaseClass;
VoxFilterAlgorithmCustom(const Vector3i &size) :
BaseClass(size), m_CustomEvaluate(NULL)
{}
typedef VoxImageFilter<VoxelT, VoxFilterAlgorithmCustom<VoxelT>> BaseClass;
VoxFilterAlgorithmCustom(const Vector3i &size)
: BaseClass(size), m_CustomEvaluate(NULL) {}
float Evaluate(const VoxImage<VoxelT> &buffer, int index)
{
if(likely(m_CustomEvaluate)) {
const std::vector<VoxelT> &vbuf = buffer.ConstData();
const std::vector<VoxelT> &vker = this->m_KernelData.ConstData();
int vox_size = vbuf.size();
int ker_size = vker.size();
int pos;
float Evaluate(const VoxImage<VoxelT> &buffer, int index) {
if (likely(m_CustomEvaluate)) {
const DataAllocator<VoxelT> &vbuf = buffer.ConstData();
const DataAllocator<VoxelT> &vker = this->m_KernelData.ConstData();
int vox_size = vbuf.size();
int ker_size = vker.size();
int pos;
float ker_sum = 0;
std::vector<Scalarf> mfh(ker_size);
for (int ik = 0; ik < ker_size; ik++) {
pos = index + vker[ik].Count - vker[this->m_KernelData.GetCenterData()].Count;
pos = (pos + vox_size) % vox_size;
mfh[ik] = vbuf[pos].Value * vker[ik].Value;
ker_sum += vker[ik].Value;
}
return this->m_CustomEvaluate(mfh);
}
else
std::cerr << "Custom evaluate function is NULL \n" <<
"No operation performed by filter.\n";
float ker_sum = 0;
std::vector<Scalarf> mfh(ker_size);
for (int ik = 0; ik < ker_size; ik++) {
pos = index + vker[ik].Count -
vker[this->m_KernelData.GetCenterData()].Count;
pos = (pos + vox_size) % vox_size;
mfh[ik] = vbuf[pos].Value * vker[ik].Value;
ker_sum += vker[ik].Value;
}
return this->m_CustomEvaluate(mfh);
} else {
std::cerr << "Custom evaluate function is NULL \n"
<< "No operation performed by filter.\n";
return 0;
}
}
inline void SetCustomEvaluate(FunctionPt funPt) { this->m_CustomEvaluate = funPt; }
inline void SetCustomEvaluate(FunctionPt funPt) {
this->m_CustomEvaluate = funPt;
}
private:
FunctionPt m_CustomEvaluate;
FunctionPt m_CustomEvaluate;
};
}
} // namespace uLib
#endif // VOXIMAGEFILTERCUSTOM_HPP

98
src/Math/VoxImageFilterLinear.hpp
View File

@@ -23,14 +23,12 @@
//////////////////////////////////////////////////////////////////////////////*/
#ifndef VOXIMAGEFILTERLINEAR_HPP
#define VOXIMAGEFILTERLINEAR_HPP
#include <Math/Dense.h>
#include "Math/VoxImage.h"
#include "VoxImageFilter.h"
#include <Math/Dense.h>
////////////////////////////////////////////////////////////////////////////////
///// VOXIMAGE FILTER LINEAR /////////////////////////////////////////////////
@@ -38,32 +36,86 @@
namespace uLib {
#ifdef USE_CUDA
template <typename VoxelT>
__global__ void LinearFilterKernel(const VoxelT *in, VoxelT *out,
const VoxelT *kernel, int vox_size,
int ker_size, int center_count) {
int index = blockIdx.x * blockDim.x + threadIdx.x;
if (index < vox_size) {
float conv = 0;
float ksum = 0;
for (int ik = 0; ik < ker_size; ++ik) {
int pos = index + kernel[ik].Count - center_count;
if (pos < 0) {
pos += vox_size * ((-pos / vox_size) + 1);
}
pos = pos % vox_size;
conv += in[pos].Value * kernel[ik].Value;
ksum += kernel[ik].Value;
}
out[index].Value = conv / ksum;
}
}
#endif
template <typename VoxelT>
class VoxFilterAlgorithmLinear :
public VoxImageFilter<VoxelT, VoxFilterAlgorithmLinear<VoxelT> > {
class VoxFilterAlgorithmLinear
: public VoxImageFilter<VoxelT, VoxFilterAlgorithmLinear<VoxelT>> {
public:
typedef VoxImageFilter<VoxelT, VoxFilterAlgorithmLinear<VoxelT> > BaseClass;
VoxFilterAlgorithmLinear(const Vector3i &size) : BaseClass(size) {}
typedef VoxImageFilter<VoxelT, VoxFilterAlgorithmLinear<VoxelT>> BaseClass;
VoxFilterAlgorithmLinear(const Vector3i &size) : BaseClass(size) {}
float Evaluate(const VoxImage<VoxelT> &buffer, int index)
{
const std::vector<VoxelT> &vbuf = buffer.ConstData();
const std::vector<VoxelT> &vker = this->m_KernelData.ConstData();
int vox_size = vbuf.size();
int ker_size = vker.size();
int pos;
float conv = 0, ksum = 0;
for (int ik = 0; ik < ker_size; ++ik) {
pos = index + vker[ik].Count - vker[this->m_KernelData.GetCenterData()].Count;
pos = (pos + vox_size) % vox_size;
conv += vbuf[pos].Value * vker[ik].Value;
ksum += vker[ik].Value;
}
return conv / ksum;
#ifdef USE_CUDA
void Run() {
if (this->m_Image->Data().GetDevice() == MemoryDevice::VRAM ||
this->m_KernelData.Data().GetDevice() == MemoryDevice::VRAM) {
this->m_Image->Data().MoveToVRAM();
this->m_KernelData.Data().MoveToVRAM();
VoxImage<VoxelT> buffer = *(this->m_Image);
buffer.Data().MoveToVRAM();
int vox_size = buffer.Data().size();
int ker_size = this->m_KernelData.Data().size();
VoxelT *d_img_out = this->m_Image->Data().GetVRAMData();
const VoxelT *d_img_in = buffer.Data().GetVRAMData();
const VoxelT *d_kernel = this->m_KernelData.Data().GetVRAMData();
int center_count =
this->m_KernelData[this->m_KernelData.GetCenterData()].Count;
int threadsPerBlock = 256;
int blocksPerGrid = (vox_size + threadsPerBlock - 1) / threadsPerBlock;
LinearFilterKernel<<<blocksPerGrid, threadsPerBlock>>>(
d_img_in, d_img_out, d_kernel, vox_size, ker_size, center_count);
cudaDeviceSynchronize();
} else {
BaseClass::Run();
}
}
#endif
float Evaluate(const VoxImage<VoxelT> &buffer, int index) {
const DataAllocator<VoxelT> &vbuf = buffer.ConstData();
const DataAllocator<VoxelT> &vker = this->m_KernelData.ConstData();
int vox_size = vbuf.size();
int ker_size = vker.size();
int pos;
float conv = 0, ksum = 0;
for (int ik = 0; ik < ker_size; ++ik) {
pos = index + vker[ik].Count -
vker[this->m_KernelData.GetCenterData()].Count;
pos = (pos + vox_size) % vox_size;
conv += vbuf[pos].Value * vker[ik].Value;
ksum += vker[ik].Value;
}
return conv / ksum;
}
};
}
} // namespace uLib
#endif // VOXIMAGEFILTERLINEAR_HPP

57
src/Math/VoxImageFilterMedian.hpp
View File

@@ -23,14 +23,12 @@
//////////////////////////////////////////////////////////////////////////////*/
#ifndef VOXIMAGEFILTERMEDIAN_HPP
#define VOXIMAGEFILTERMEDIAN_HPP
#include <Math/Dense.h>
#include "Math/VoxImage.h"
#include "VoxImageFilter.h"
#include <Math/Dense.h>
////////////////////////////////////////////////////////////////////////////////
///// VOXIMAGE FILTER MEDIAN /////////////////////////////////////////////////
@@ -39,37 +37,38 @@
namespace uLib {
template <typename VoxelT>
class VoxFilterAlgorithmMedian :
public VoxImageFilter<VoxelT, VoxFilterAlgorithmMedian<VoxelT> > {
class VoxFilterAlgorithmMedian
: public VoxImageFilter<VoxelT, VoxFilterAlgorithmMedian<VoxelT>> {
public:
typedef VoxImageFilter<VoxelT, VoxFilterAlgorithmMedian<VoxelT> > BaseClass;
VoxFilterAlgorithmMedian(const Vector3i &size) : BaseClass(size) {}
typedef VoxImageFilter<VoxelT, VoxFilterAlgorithmMedian<VoxelT>> BaseClass;
VoxFilterAlgorithmMedian(const Vector3i &size) : BaseClass(size) {}
float Evaluate(const VoxImage<VoxelT> &buffer, int index)
{
const std::vector<VoxelT> &vbuf = buffer.ConstData();
const std::vector<VoxelT> &vker = this->m_KernelData.ConstData();
int vox_size = vbuf.size();
int ker_size = vker.size();
int pos;
float Evaluate(const VoxImage<VoxelT> &buffer, int index) {
const DataAllocator<VoxelT> &vbuf = buffer.ConstData();
const DataAllocator<VoxelT> &vker = this->m_KernelData.ConstData();
int vox_size = vbuf.size();
int ker_size = vker.size();
int pos;
std::vector<float> mfh(ker_size);
for (int ik = 0; ik < ker_size; ik++) {
pos = index + vker[ik].Count - vker[this->m_KernelData.GetCenterData()].Count;
pos = (pos + vox_size) % vox_size;
mfh[ik] = vbuf[pos].Value * vker[ik].Value;
}
std::sort(mfh.begin(), mfh.end());
pos = 0;
// count zeroes in filter kernel to move it out of median //
for (int i = 0; i < ker_size; ++i)
if (vker[i].Value == 0.0) pos++;
// median //
pos += (ker_size - pos) / 2;
return mfh[pos];
std::vector<float> mfh(ker_size);
for (int ik = 0; ik < ker_size; ik++) {
pos = index + vker[ik].Count -
vker[this->m_KernelData.GetCenterData()].Count;
pos = (pos + vox_size) % vox_size;
mfh[ik] = vbuf[pos].Value * vker[ik].Value;
}
std::sort(mfh.begin(), mfh.end());
pos = 0;
// count zeroes in filter kernel to move it out of median //
for (int i = 0; i < ker_size; ++i)
if (vker[i].Value == 0.0)
pos++;
// median //
pos += (ker_size - pos) / 2;
return mfh[pos];
}
};
}
} // namespace uLib
#endif // VOXIMAGEFILTERMEDIAN_HPP

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

@@ -22,3 +22,7 @@ set(LIBRARIES
)
uLib_add_tests(Math)
if(USE_CUDA)
set_source_files_properties(VoxImageFilterTest.cpp PROPERTIES LANGUAGE CUDA)
endif()

215
src/Math/testing/VoxImageFilterTest.cpp
View File

@@ -23,128 +23,191 @@
//////////////////////////////////////////////////////////////////////////////*/
#include "testing-prototype.h"
#include "Math/StructuredGrid.h"
#include "testing-prototype.h"
#include "Math/VoxImage.h"
#include "Math/VoxImageFilter.h"
using namespace uLib;
struct TestVoxel {
Scalarf Value;
unsigned int Count;
Scalarf Value;
unsigned int Count;
};
float GaussianShape(float d)
{
// normalized manually .. fix //
return 4.5 * exp(-d * 4.5);
float GaussianShape(float d) {
// normalized manually .. fix //
return 4.5 * exp(-d * 4.5);
}
class GaussianShapeClass : public Interface::VoxImageFilterShape {
public:
GaussianShapeClass(float sigma) :
m_sigma(sigma)
{}
GaussianShapeClass(float sigma) : m_sigma(sigma) {}
float operator ()(float d) {
return (1/m_sigma) * exp(-d/m_sigma);
}
float operator()(float d) { return (1 / m_sigma) * exp(-d / m_sigma); }
private:
float m_sigma;
float m_sigma;
};
static float MaxInVector(const std::vector<float> &v)
{
float max = 0;
for(int i=0; i<v.size(); ++i)
if(v.at(i) > max) max = v.at(i);
return max;
static float MaxInVector(const std::vector<float> &v) {
float max = 0;
for (int i = 0; i < v.size(); ++i)
if (v.at(i) > max)
max = v.at(i);
return max;
}
int main() {
BEGIN_TESTING(VoxImageFilters);
int main()
{
BEGIN_TESTING(VoxImageFilters);
VoxImage<TestVoxel> image(Vector3i(20, 30, 40));
image[Vector3i(10, 10, 10)].Value = 1;
// image[Vector3i(10,10,8)].Value = 1;
image.ExportToVtk("test_filter_original.vtk", 0);
VoxImage<TestVoxel> image(Vector3i(20,30,40));
image[Vector3i(10,10,10)].Value = 1;
//image[Vector3i(10,10,8)].Value = 1;
image.ExportToVtk("test_filter_original.vtk",0);
////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////
// RPS //
{
VoxFilterAlgorithmSPR<TestVoxel> filter(Vector3i(2, 3, 4));
VoxImage<TestVoxel> filtered = image;
////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////
// RPS //
std::vector<float> values;
for (int i = 0; i < filter.GetKernelData().GetDims().prod(); ++i) {
values.push_back(1.);
std::cout << values[i] << " ";
}
std::cout << "\n";
{
VoxFilterAlgorithmSPR<TestVoxel> filter(Vector3i(2,3,4));
filter.SetImage(&filtered);
VoxImage<TestVoxel> filtered = image;
filter.SetKernelNumericXZY(values);
std::vector<float> values;
for(int i=0; i < filter.GetKernelData().GetDims().prod(); ++i) {
values.push_back(1.);
std::cout << values[i] << " ";
}
std::cout << "\n";
filter.SetABTrim(0, 2);
filter.SetImage(&filtered);
filter.GetKernelData().PrintSelf(std::cout);
filter.SetKernelNumericXZY(values);
filter.Run();
filter.SetABTrim(0,2);
filtered.ExportToVtk("filter_RPS_out.vtk", 0);
}
filter.GetKernelData().PrintSelf(std::cout);
{
filter.Run();
VoxImage<TestVoxel> image(Vector3i(20, 30, 40));
image[Vector3i(10, 10, 10)].Value = 1;
image[Vector3i(9, 10, 8)].Value = 2;
image.ExportToVtk("test_filter_max_original.vtk", 0);
filtered.ExportToVtk("filter_RPS_out.vtk",0);
VoxFilterAlgorithmCustom<TestVoxel> filter(Vector3i(3, 3, 4));
std::vector<float> values;
for (int i = 0; i < filter.GetKernelData().GetDims().prod(); ++i) {
values.push_back(static_cast<float>(1));
}
filter.SetImage(&image);
filter.SetKernelNumericXZY(values);
filter.SetCustomEvaluate(MaxInVector);
filter.Run();
{
image.ExportToVtk("test_filter_max.vtk", 0);
}
VoxImage<TestVoxel> image(Vector3i(20,30,40));
image[Vector3i(10,10,10)].Value = 1;
image[Vector3i(9,10,8)].Value = 2;
image.ExportToVtk("test_filter_max_original.vtk",0);
////////////////////////////////////////////////////////////////////////////
// CUDA Allocator Transfer Test //
{
VoxImage<TestVoxel> image(Vector3i(10, 10, 10));
image[Vector3i(5, 5, 5)].Value = 1;
VoxFilterAlgorithmLinear<TestVoxel> filter(Vector3i(3, 3, 3));
VoxFilterAlgorithmCustom<TestVoxel> filter(Vector3i(3,3,4));
std::vector<float> values;
for(int i=0; i < filter.GetKernelData().GetDims().prod(); ++i) {
values.push_back(static_cast<float>(1));
}
filter.SetImage(&image);
filter.SetKernelNumericXZY(values);
filter.SetCustomEvaluate(MaxInVector);
filter.Run();
image.ExportToVtk("test_filter_max.vtk",0);
std::vector<float> values;
for (int i = 0; i < filter.GetKernelData().GetDims().prod(); ++i) {
values.push_back(1.0f);
}
filter.SetImage(&image);
filter.SetKernelNumericXZY(values);
// Move the kernel data and image data to VRAM to simulate CUDA transfer
filter.GetKernelData().Data().MoveToVRAM();
image.Data().MoveToVRAM();
END_TESTING;
// Validate devices
if (filter.GetKernelData().Data().GetDevice() != MemoryDevice::VRAM ||
image.Data().GetDevice() != MemoryDevice::VRAM) {
#ifdef USE_CUDA
std::cerr << "Failed to move memory to VRAM." << std::endl;
#else
std::cout << "DataAllocator correctly simulates VRAM without crashing."
<< std::endl;
#endif
}
// Run the filter; The fallback CPU filter will trigger MoveToRAM
// behind the scenes inside Convolve / Evaluate.
filter.Run();
// Assert it came back to RAM if evaluated on CPU
if (image.Data().GetDevice() != MemoryDevice::RAM) {
#ifdef USE_CUDA
std::cout << "Data correctly stayed in VRAM after CUDA execution!"
<< std::endl;
#else
std::cout << "Data correctly stayed in RAM simulation." << std::endl;
#endif
}
image.ExportToVtk("test_filter_cuda_transfer.vtk", 0);
}
////////////////////////////////////////////////////////////////////////////
// CUDA ABTrim Allocator Transfer Test //
{
VoxImage<TestVoxel> image(Vector3i(10, 10, 10));
image[Vector3i(5, 5, 5)].Value = 10;
image[Vector3i(5, 5, 6)].Value = 2; // Test trimming
VoxFilterAlgorithmAbtrim<TestVoxel> filter(Vector3i(3, 3, 3));
std::vector<float> values;
for (int i = 0; i < filter.GetKernelData().GetDims().prod(); ++i) {
values.push_back(1.0f);
}
filter.SetImage(&image);
filter.SetKernelNumericXZY(values);
filter.SetABTrim(1, 1); // trim highest and lowest
// Move the kernel data and image data to VRAM to simulate CUDA transfer
filter.GetKernelData().Data().MoveToVRAM();
image.Data().MoveToVRAM();
// Run the filter
filter.Run();
// Ensure data stays on device if CUDA was toggled
if (image.Data().GetDevice() != MemoryDevice::RAM) {
#ifdef USE_CUDA
std::cout << "ABTrim correctly stayed in VRAM after CUDA execution!"
<< std::endl;
#else
std::cout << "ABTrim Data correctly stayed in RAM simulation."
<< std::endl;
#endif
}
image.ExportToVtk("test_filter_abtrim_cuda_transfer.vtk", 0);
}
END_TESTING;
}