refactor: migrate voxel data storage to DataAllocator for CUDA

CMakeLists.txt

@@ -9,6 +9,13 @@ cmake_minimum_required (VERSION 3.26)
 
 project(uLib)
 
+# CUDA Toolkit seems to be missing locally. Toggle ON if nvcc is made available.
+option(USE_CUDA "Enable CUDA support" ON)
+if(USE_CUDA)
+    enable_language(CUDA)
+    add_compile_definitions(USE_CUDA)
+endif()
+
 # The version number.
 set(PROJECT_VERSION_MAJOR 0)
 set(PROJECT_VERSION_MINOR 6)

docs/usage/usage.md

@@ -0,0 +1,60 @@
+# Usage and Installation Guide
+
+## Requirements
+
+### Compiling with CUDA Support
+
+The library supports running VoxImage filtering operations directly on CUDA cores via transparent RAM/VRAM memory transfers.
+
+By default, the `CMakeLists.txt` build system sets `USE_CUDA=ON` and will attempt to locate `nvcc` and the NVIDIA CUDA Toolkit. If the toolkit is missing, `CMake` will fail unless you explicitly configure the project with `-DUSE_CUDA=OFF`.
+
+### 1. Installing CUDA Environment via Micromamba
+
+If you are developing inside an isolated Conda/Micromamba environment (e.g., `mutom`), you can inject the CUDA compilers directly into your environment rather than relying on global system dependencies:
+
+```bash
+# Add the conda-forge channel if not already available
+micromamba config append channels conda-forge
+
+# Install nvcc and the necessary CUDA toolkit components
+micromamba install cuda-nvcc
+```
+
+Verify your installation:
+
+```bash
+nvcc --version
+```
+
+### 2. Building the Project
+
+Configure and compile the project using standard CMake flows:
+
+```bash
+mkdir -p build && cd build
+
+# Configure CMake
+# (Optional) Explicitly toggle CUDA: cmake -DUSE_CUDA=ON ..
+cmake ..
+
+# Compile the project and tests
+make -j $(nproc)
+```
+
+### 3. Validating CUDA Support
+
+You can verify that the CUDA kernels are launching correctly and allocating device memory through `DataAllocator` by running the mathematical unit tests.
+
+```bash
+# From the build directory
+./src/Math/testing/VoxImageFilterTest
+
+# Output should show:
+# "Data correctly stayed in VRAM after CUDA execution!"
+```
+
+## How It Works Under The Hood
+
+The `DataAllocator<T>` container automatically wraps memory allocations to transparently map to CPU RAM, or GPU VRAM. Standard iteration automatically pulls data backwards using implicit `MoveToRAM()` calls. 
+
+Filters using `#ifdef USE_CUDA` explicitly dictate `<buffer>.MoveToVRAM()` allocating directly on device bounds seamlessly. Fallbacks to Host compute iterations handle themselves automatically. Chaining specific filters together safely chains continuous VRAM operations avoiding costly Host copies in between iterations.

src/Core/Object.h

@@ -116,16 +116,15 @@ public:
   connect(typename FunctionPointer<Func1>::Object *sender, Func1 sigf,
           typename FunctionPointer<Func2>::Object *receiver, Func2 slof) {
     SignalBase *sigb = sender->findOrAddSignal(sigf);
-    typedef boost::signals2::signal<
+    ConnectSignal<typename FunctionPointer<Func1>::SignalSignature>(sigb, slof,
-        typename FunctionPointer<Func2>::SignalSignature>
+                                                                    receiver);
-        SigT;
-    ConnectSignal(sigb, slof, receiver);
     return true;
   }
 
   template <typename FuncT>
   static inline bool connect(SignalBase *sigb, FuncT slof, Object *receiver) {
-    ConnectSignal(sigb, slof, receiver);
+    ConnectSignal<typename FunctionPointer<FuncT>::SignalSignature>(sigb, slof,
+                                                                    receiver);
     return true;
   }
 

src/Core/Signal.h

@@ -50,8 +50,11 @@ using namespace boost::placeholders;
 #define SIGNAL(a) BOOST_STRINGIZE(a)
 
 #define _ULIB_DETAIL_SIGNAL_EMIT(_name, ...)                                   \
-  static BOOST_AUTO(sig, this->findOrAddSignal(&_name));                       \
+  do {                                                                         \
-  sig->operator()(__VA_ARGS__);
+    BOOST_AUTO(sig, this->findOrAddSignal(&_name));                            \
+    if (sig)                                                                   \
+      sig->operator()(__VA_ARGS__);                                            \
+  } while (0)
 
 /**
  * Utility macro to implement signal emission implementa una delle seguenti:
@@ -84,66 +87,61 @@ template <typename T> struct Signal {
 
 namespace detail {
 
-template <typename FuncT, int arity> struct ConnectSignal {};
+template <typename FuncT, typename SigSignature, int arity>
+struct ConnectSignal {};
 
-template <typename FuncT> struct ConnectSignal<FuncT, 0> {
+template <typename FuncT, typename SigSignature>
+struct ConnectSignal<FuncT, SigSignature, 0> {
   static void connect(SignalBase *sigb, FuncT slof,
                       typename FunctionPointer<FuncT>::Object *receiver) {
-    typedef
+    typedef typename Signal<SigSignature>::type SigT;
-        typename Signal<typename FunctionPointer<FuncT>::SignalSignature>::type
-            SigT;
     reinterpret_cast<SigT *>(sigb)->connect(slof);
   }
 };
 
-template <typename FuncT> struct ConnectSignal<FuncT, 1> {
+template <typename FuncT, typename SigSignature>
+struct ConnectSignal<FuncT, SigSignature, 1> {
   static void connect(SignalBase *sigb, FuncT slof,
                       typename FunctionPointer<FuncT>::Object *receiver) {
-    typedef
+    typedef typename Signal<SigSignature>::type SigT;
-        typename Signal<typename FunctionPointer<FuncT>::SignalSignature>::type
-            SigT;
     reinterpret_cast<SigT *>(sigb)->connect(boost::bind(slof, receiver));
   }
 };
 
-template <typename FuncT> struct ConnectSignal<FuncT, 2> {
+template <typename FuncT, typename SigSignature>
+struct ConnectSignal<FuncT, SigSignature, 2> {
   static void connect(SignalBase *sigb, FuncT slof,
                       typename FunctionPointer<FuncT>::Object *receiver) {
-    typedef
+    typedef typename Signal<SigSignature>::type SigT;
-        typename Signal<typename FunctionPointer<FuncT>::SignalSignature>::type
-            SigT;
     reinterpret_cast<SigT *>(sigb)->connect(boost::bind(slof, receiver, _1));
   }
 };
 
-template <typename FuncT> struct ConnectSignal<FuncT, 3> {
+template <typename FuncT, typename SigSignature>
+struct ConnectSignal<FuncT, SigSignature, 3> {
   static void connect(SignalBase *sigb, FuncT slof,
                       typename FunctionPointer<FuncT>::Object *receiver) {
-    typedef
+    typedef typename Signal<SigSignature>::type SigT;
-        typename Signal<typename FunctionPointer<FuncT>::SignalSignature>::type
-            SigT;
     reinterpret_cast<SigT *>(sigb)->connect(
         boost::bind(slof, receiver, _1, _2));
   }
 };
 
-template <typename FuncT> struct ConnectSignal<FuncT, 4> {
+template <typename FuncT, typename SigSignature>
+struct ConnectSignal<FuncT, SigSignature, 4> {
   static void connect(SignalBase *sigb, FuncT slof,
                       typename FunctionPointer<FuncT>::Object *receiver) {
-    typedef
+    typedef typename Signal<SigSignature>::type SigT;
-        typename Signal<typename FunctionPointer<FuncT>::SignalSignature>::type
-            SigT;
     reinterpret_cast<SigT *>(sigb)->connect(
         boost::bind(slof, receiver, _1, _2, _3));
   }
 };
 
-template <typename FuncT> struct ConnectSignal<FuncT, 5> {
+template <typename FuncT, typename SigSignature>
+struct ConnectSignal<FuncT, SigSignature, 5> {
   static void connect(SignalBase *sigb, FuncT slof,
                       typename FunctionPointer<FuncT>::Object *receiver) {
-    typedef
+    typedef typename Signal<SigSignature>::type SigT;
-        typename Signal<typename FunctionPointer<FuncT>::SignalSignature>::type
-            SigT;
     reinterpret_cast<SigT *>(sigb)->connect(
         boost::bind(slof, receiver, _1, _2, _3, _4));
   }
@@ -152,15 +150,16 @@ template <typename FuncT> struct ConnectSignal<FuncT, 5> {
 } // namespace detail
 
 template <typename FuncT> SignalBase *NewSignal(FuncT f) {
-  // seems to work wow !
+  return new
-  return new Signal<void()>::type;
+      typename Signal<typename FunctionPointer<FuncT>::SignalSignature>::type;
 }
 
-template <typename FuncT>
+template <typename SigSignature, typename FuncT>
 void ConnectSignal(SignalBase *sigb, FuncT slof,
                    typename FunctionPointer<FuncT>::Object *receiver) {
-  detail::ConnectSignal<FuncT, FunctionPointer<FuncT>::arity>::connect(
+  detail::ConnectSignal<FuncT, SigSignature,
-      sigb, slof, receiver);
+                        FunctionPointer<FuncT>::arity>::connect(sigb, slof,
+                                                                receiver);
 }
 
 } // namespace uLib

src/Core/testing/CMakeLists.txt

@@ -19,6 +19,7 @@ set( TESTS
         UuidTest
         TypeIntrospectionTraversal
         OptionsTest
+        PingPongTest
 )
 
 set(LIBRARIES

src/Core/testing/PingPongTest.cpp

@@ -0,0 +1,52 @@
+#include "Core/Object.h"
+#include "Core/Signal.h"
+#include "testing-prototype.h"
+#include <iostream>
+
+using namespace uLib;
+
+class Ping : public Object {
+public:
+signals:
+  void PingSignal(int count);
+public slots:
+  void OnPong(int count) {
+    std::cout << "Ping received Pong " << count << std::endl;
+    if (count > 0)
+      ULIB_SIGNAL_EMIT(Ping::PingSignal, count - 1);
+  }
+};
+
+void Ping::PingSignal(int count) { ULIB_SIGNAL_EMIT(Ping::PingSignal, count); }
+
+class Pong : public Object {
+public:
+signals:
+  void PongSignal(int count);
+public slots:
+  void OnPing(int count) {
+    std::cout << "Pong received Ping " << count << std::endl;
+    if (count > 0)
+      ULIB_SIGNAL_EMIT(Pong::PongSignal, count - 1);
+  }
+};
+
+void Pong::PongSignal(int count) { ULIB_SIGNAL_EMIT(Pong::PongSignal, count); }
+
+int main() {
+  BEGIN_TESTING(PingPong);
+  Ping ping;
+  Pong pong;
+
+  std::cout << "Connecting ping to pong" << std::endl;
+  Object::connect(&ping, &Ping::PingSignal, &pong, &Pong::OnPing);
+
+  std::cout << "Connecting pong to ping" << std::endl;
+  Object::connect(&pong, &Pong::PongSignal, &ping, &Ping::OnPong);
+
+  std::cout << "Emitting PingSignal(5)" << std::endl;
+  ping.PingSignal(5);
+
+  END_TESTING;
+  return 0;
+}

src/Core/testing/SignalTest.cpp

@@ -23,93 +23,63 @@
 
 //////////////////////////////////////////////////////////////////////////////*/
 
-
-
 #include <iostream>
 #include <typeinfo>
 
-
-#include "testing-prototype.h"
-#include "Core/Types.h"
 #include "Core/Object.h"
 #include "Core/Signal.h"
-
+#include "Core/Types.h"
+#include "testing-prototype.h"
 
 using namespace uLib;
 
-
-
-
 class Ob1 : public Object {
 public:
-
-
 signals:
-    void V0();
+  void V0();
-
-    int V1(int a);
 
+  void V1(int a);
 };
 
-
-
-
 // should be done by moc //
-void Ob1::V0() {
+void Ob1::V0() { ULIB_SIGNAL_EMIT(Ob1::V0); }
-    ULIB_SIGNAL_EMIT(Ob1::V0);
-}
-
-int Ob1::V1(int a) {
-    ULIB_SIGNAL_EMIT(Ob1::V1,a);
-}
-
-
-
 
+void Ob1::V1(int a) { ULIB_SIGNAL_EMIT(Ob1::V1, a); }
 
 class Ob2 : public Object {
 
 public slots:
-    void PrintV0() {
+  void PrintV0() { std::cout << "Ob2 prints V0\n" << std::flush; }
-        std::cout << "Ob2 prints V0\n" << std::flush;
-    }
 };
 
 class Ob3 : public Object {
 
 public slots:
-    void PrintV0() {
+  void PrintV0() { std::cout << "Ob3 prints V0\n" << std::flush; }
-        std::cout << "Ob3 prints V0\n" << std::flush;
-    }
 
-    void PrintNumber(int n) {
+  void PrintNumber(int n) {
-        std::cout << "Ob3 is printing number: " << n << "\n";
+    std::cout << "Ob3 is printing number: " << n << "\n";
-    }
+  }
 };
 
-
-
-
 int main() {
-    BEGIN_TESTING(Signals);
+  BEGIN_TESTING(Signals);
 
-    Ob1 ob1;
+  Ob1 ob1;
-    Ob2 ob2;
+  Ob2 ob2;
-    Ob3 ob3;
+  Ob3 ob3;
 
-    Object::connect(&ob1,&Ob1::V0,&ob2,&Ob2::PrintV0);
+  Object::connect(&ob1, &Ob1::V0, &ob2, &Ob2::PrintV0);
-    Object::connect(&ob1,&Ob1::V0,&ob3,&Ob3::PrintV0);
+  Object::connect(&ob1, &Ob1::V0, &ob3, &Ob3::PrintV0);
-    Object::connect(&ob1,&Ob1::V1,&ob3,&Ob3::PrintNumber);
+  Object::connect(&ob1, &Ob1::V1, &ob3, &Ob3::PrintNumber);
 
-    //  not working yet
+  //  not working yet
-    //    Object::connect(&ob1,SIGNAL(V0(),&ob2,SLOT(PrintV0())
+  //    Object::connect(&ob1,SIGNAL(V0(),&ob2,SLOT(PrintV0())
 
-    ob1.PrintSelf(std::cout);
+  ob1.PrintSelf(std::cout);
 
-    emit ob1.V0();
+  emit ob1.V0();
-    emit ob1.V1(5552368);
+  emit ob1.V1(5552368);
 
-    END_TESTING;
+  END_TESTING;
 }
-
-

src/Math/DataAllocator.h

@@ -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

src/Math/VoxImage.h

@@ -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 Vector3i &id) const = 0;
-    virtual float GetValue(const int id) const = 0;
+  virtual float GetValue(const int id) const = 0;
-    virtual void SetValue(const Vector3i &id, float value) = 0;
+  virtual void SetValue(const Vector3i &id, float value) = 0;
-    virtual void SetValue(const int 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
+  // use this function to export to VTK binary format
-    void ExportToVti (const char *file, bool density_type = 0, bool compressed = 0);
+  void ExportToVti(const char *file, bool density_type = 0,
+                   bool compressed = 0);
 
-    // this function has been deprecated in favor of ExportToVti
+  // this function has been deprecated in favor of ExportToVti
-    // but it is kept for backward compatibility and because it 
+  // but it is kept for backward compatibility and because it
-    // does not depend on vtk library
+  // does not depend on vtk library
-    void ExportToVtkXml(const char *file, bool density_type = 0);
+  void ExportToVtkXml(const char *file, bool density_type = 0);
 
-    int ImportFromVtk(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);
+  int ImportFromVti(const char *file, bool density_type = 0);
 
 protected:
-
+  virtual ~VoxImage() {}
-    virtual ~VoxImage() {}
+  VoxImage(const Vector3i &size) : BaseClass(size) {}
-    VoxImage(const Vector3i &size) : BaseClass(size) {}
 };
 
-}
+} // namespace Abstract
 
 ////////////////////////////////////////////////////////////////////////////////
 //  VOXEL   ////////////////////////////////////////////////////////////////////
@@ -83,421 +83,415 @@ protected:
 
 namespace Interface {
 struct Voxel {
-    template<class Self> void check_structural() {
+  template <class Self> void check_structural() {
-        uLibCheckMember(Self,Value, Scalarf);
+    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) :
+  VoxImage(const VoxImage<T> &copy) : BaseClass(copy) {
-        BaseClass(copy)
+    this->m_Data = copy.m_Data;
-    {
+  }
-        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; }
+  template <typename S> void operator-=(VoxImage<S> &sibling) {
-    inline const std::vector<T>& ConstData() const { return m_Data; }
+    if (this->GetDims() != sibling.GetDims()) {
-
+      // printf("Warning when subtracting VoxImages: I'm NOT doing it!\n");
-    inline const T& At(int i)              const { return m_Data.at(i); }
+      return;
-    inline const T& At(const Vector3i &id) const { return m_Data.at(Map(id)); }
+    } // WARNING! You must Warn the user!
-    inline T& operator[](unsigned int i)     { return m_Data[i]; }
+    for (unsigned int i = 0; i < m_Data.size(); ++i) {
-    inline T& operator[](const Vector3i &id) { return m_Data[Map(id)]; }
+      m_Data[i].Value -= sibling.At(i).Value;
-
-    // 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) {
+  template <typename S> void operator*=(VoxImage<S> &sibling) {
-        this->operator [](id).Value = value;
+    if (this->GetDims() != sibling.GetDims()) {
-    }
+      // printf("Warning when multiplying VoxImages: I'm NOT doing it!\n");
-    inline void SetValue(const int id, float value) {
+      return;
-        this->operator [](id).Value = value;
+    } // 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) {
+  template <typename S> void operator/=(VoxImage<S> &sibling) {
-        this->m_Data.resize(size.prod());
+    if (this->GetDims() != sibling.GetDims()) {
-        BaseClass::BaseClass::SetDims(size); // FIX horrible coding style !
+      // printf("Warning when dividing VoxImages: I'm NOT doing it!\n");
-    }
+      return;
-
+    } // WARNING! You must Warn the user!
-    inline VoxImage<T> clipImage(const Vector3i begin, const Vector3i end) const;
+    for (unsigned int i = 0; i < m_Data.size(); ++i) {
-    inline VoxImage<T> clipImage(const HPoint3f begin, const HPoint3f end) const;
+      m_Data[i].Value /= sibling.At(i).Value;
-    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;
-        }
     }
+  }
 
 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
+VoxImage<T>::VoxImage() : m_Data(0), BaseClass(Vector3i(0, 0, 0)) {
-{
+  Interface::IsA<T, Interface::Voxel>(); /* structural check for T */
-        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
+VoxImage<T>::VoxImage(const Vector3i &size)
-{
+    : m_Data(size.prod()), BaseClass(size) {
-    Vector3i v1 = this->Find(begin);
+  Interface::IsA<T, Interface::Voxel>(); /* structural check for T */
-    Vector3i v2 = this->Find(end);
-    return this->clipImage(v1,v2);
 }
 
 template <typename T>
-VoxImage<T> VoxImage<T>::clipImage(const float density) const
+VoxImage<T> VoxImage<T>::clipImage(const Vector3i begin,
-{
+                                   const Vector3i end) const {
-    Vector3i v1 = this->GetDims();
+  Vector3i dim = (end - begin) + Vector3i(1, 1, 1);
-    Vector3i v2 = Vector3i(0,0,0);
+  VoxImage<T> out(*this);
-    for(uint i=0; i< this->m_Data.size(); ++i) {
+  out.SetDims(dim);
-        if( this->GetValue(i) >= density ) {
+  out.SetPosition(this->GetPosition() +
-            Vector3i id = this->UnMap(i);
+                  this->GetSpacing().cwiseProduct(begin.cast<float>()));
-            v1 = v1.array().min(id.array());
+
-            v2 = v2.array().max(id.array());
+  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
+VoxImage<T> VoxImage<T>::clipImage(const float densityMin,
-{
+                                   const float densityMax) const {
-    Vector3i v1 = this->GetDims();
+  Vector3i v1 = this->GetDims();
-    Vector3i v2 = Vector3i(0,0,0);
+  Vector3i v2 = Vector3i(0, 0, 0);
-    for(uint i=0; i< this->m_Data.size(); ++i) {
+  for (uint i = 0; i < this->m_Data.size(); ++i) {
-        if( this->GetValue(i) >= densityMin && this->GetValue(i) <= densityMax) {
+    if (this->GetValue(i) >= densityMin && this->GetValue(i) <= densityMax) {
-            Vector3i id = this->UnMap(i);
+      Vector3i id = this->UnMap(i);
-            v1 = v1.array().min(id.array());
+      v1 = v1.array().min(id.array());
-            v2 = v2.array().max(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> VoxImage<T>::maskImage(const HPoint3f begin, const HPoint3f end,
-{
+                                   float value) const {
-    VoxImage<T> out(*this);
+  VoxImage<T> out(*this);
-    out.SetDims(this->GetDims());
+  out.SetDims(this->GetDims());
-    out.SetPosition(this->GetPosition());
+  out.SetPosition(this->GetPosition());
 
-    Vector3i voxB = this->Find(begin);
+  Vector3i voxB = this->Find(begin);
-    Vector3i voxE = this->Find(end);
+  Vector3i voxE = this->Find(end);
 
-    Vector3i ID;
+  Vector3i ID;
 
-    for(int ix=voxB(0); ix<voxE(0); ix++)
+  for (int ix = voxB(0); ix < voxE(0); ix++)
-        for(int iy=voxB(1); iy<voxE(1); iy++)
+    for (int iy = voxB(1); iy < voxE(1); iy++)
-            for(int iz=voxB(2); iz<voxE(2); iz++){
+      for (int iz = voxB(2); iz < voxE(2); iz++) {
-                ID << ix,iy,iz;
+        ID << ix, iy, iz;
-                out.SetValue(ID,value*1.E-6);
+        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
+VoxImage<T> VoxImage<T>::maskImage(const float threshold, float belowValue,
-{
+                                   float aboveValue) const {
-    std::cout << "VoxImage: maskImage, fixing voxels under threshold " << threshold;
+  std::cout << "VoxImage: maskImage, fixing voxels under threshold "
-    if(belowValue)
+            << threshold;
-        std::cout << " at value " << belowValue;
+  if (belowValue)
-    else
+    std::cout << " at value " << belowValue;
-        std::cout << " at -value";
+  else
-    std::cout << ", voxels above threshold at value ";
+    std::cout << " at -value";
-    if(aboveValue)
+  std::cout << ", voxels above threshold at value ";
-        std::cout << aboveValue;
+  if (aboveValue)
-    else
+    std::cout << aboveValue;
-        std::cout << "found";
+  else
+    std::cout << "found";
 
+  VoxImage<T> out(*this);
+  out.SetDims(this->GetDims());
+  out.SetPosition(this->GetPosition());
 
-    VoxImage<T> out(*this);
+  for (uint i = 0; i < this->m_Data.size(); ++i) {
-    out.SetDims(this->GetDims());
+    // skip negative voxels: they are already frozen
-    out.SetPosition(this->GetPosition());
+    if (this->GetValue(i) >= 0) {
-
+      // voxels under threshold
-    for(uint i=0; i< this->m_Data.size(); ++i) {
+      if (this->GetValue(i) <= threshold * 1.E-6) {
-        // skip negative voxels: they are already frozen
+        if (belowValue) {
-        if( this->GetValue(i) >= 0 ){
+          //                std::cout << "vox " << i << ", " <<
-            // voxels under threshold
+          //                this->GetValue(i); std::cout << " ----> set to " <<
-            if( this->GetValue(i) <= threshold*1.E-6 ){
+          //                -1.*belowValue*1.E-6 << std::endl;
-                if(belowValue){
+          out.SetValue(i, -1. * belowValue * 1.E-6);
-                    //                std::cout << "vox " << i << ", " << this->GetValue(i);
+        } else
-                    //                std::cout << " ----> set to " << -1.*belowValue*1.E-6 << std::endl;
+          out.SetValue(i, -1. * this->GetValue(i));
-                    out.SetValue(i,-1.*belowValue*1.E-6);}
+      }
-                else
+      // voxels over threshold
-                    out.SetValue(i,-1.*this->GetValue(i));
+      else {
-            }
+        if (aboveValue)
-            // voxels over threshold
+          out.SetValue(i, aboveValue * 1.E-6);
-            else{
+        else
-                if(aboveValue)
+          out.SetValue(i, this->GetValue(i));
-                    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
+VoxImage<T> VoxImage<T>::fixVoxels(const float threshold,
-{
+                                   float tolerance) const {
-    std::cout << "VoxImage: fixing voxels with value " << threshold << std::endl;
+  std::cout << "VoxImage: fixing voxels with value " << threshold << std::endl;
 
-    VoxImage<T> out(*this);
+  VoxImage<T> out(*this);
-    out.SetDims(this->GetDims());
+  out.SetDims(this->GetDims());
-    out.SetPosition(this->GetPosition());
+  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 ){
+        if (fabs(this->GetValue(ID) - threshold * 1.E-6) < tolerance * 1.E-6) {
-//            std::cout << "vox " << i << ", " << this->GetValue(i);
+          out.SetValue(ID, -1. * this->GetValue(ID));
-//            std::cout << " ----> set to " << -1.*this->GetValue(i) << std::endl;
-                out.SetValue(i,-1.*this->GetValue(i));
         }
-    }
+      }
-    return out;
+
+  return out;
 }
 
 template <typename T>
-VoxImage<T> VoxImage<T>::Abs() const
+VoxImage<T> VoxImage<T>::fixVoxels(const HPoint3f begin,
-{
+                                   const HPoint3f end) const {
-    std::cout << "VoxImage: set abs voxels value " << std::endl;
+  VoxImage<T> out(*this);
+  out.SetDims(this->GetDims());
+  out.SetPosition(this->GetPosition());
 
-    VoxImage<T> out(*this);
+  Vector3i voxB = this->Find(begin);
-    out.SetDims(this->GetDims());
+  Vector3i voxE = this->Find(end);
-    out.SetPosition(this->GetPosition());
 
-    for(uint i=0; i< this->m_Data.size(); ++i)
+  Vector3i ID;
-        out.SetValue(i,fabs(this->GetValue(i)));
 
-    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> VoxImage<T>::fixVoxelsAroundPlane(const float threshold,
-{
+                                              float tolerance, const HPoint3f B,
-    VoxImage<T> out(*this);
+                                              const HPoint3f E,
-    out.SetDims(this->GetDims());
+                                              bool aboveAir) const {
-    out.SetPosition(this->GetPosition());
+  VoxImage<T> out(*this);
+  Vector3i dim = this->GetDims();
+  out.SetDims(dim);
+  out.SetPosition(this->GetPosition());
 
-    Vector3i voxB = this->Find(begin);
+  HPoint3f Bcoll = this->GetPosition().homogeneous();
-    Vector3i voxE = this->Find(end);
 
-    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++)
+        // B, E voxel position
-        for(int iy=voxB(1); iy<voxE(1); iy++)
+        Vector3i iv(ix, iy, iz);
-            for(int iz=voxB(2); iz<voxE(2); iz++){
+        Vector3f v =
-                ID << ix,iy,iz;
+            Vector3f(iv.cast<float>().cwiseProduct(this->GetSpacing()));
-                // voxels around threshold
+        HPoint3f Bvox = Bcoll + HPoint3f(v);
-                if( fabs(this->GetValue(ID) - threshold*1.E-6) < tolerance*1.E-6 ){
+        HPoint3f Evox = Bvox + this->GetSpacing().homogeneous();
-                    out.SetValue(ID,-1.*this->GetValue(ID));
+        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>
+template <typename T> void VoxImage<T>::InitVoxels(T t) {
-VoxImage<T> VoxImage<T>::fixVoxels(const HPoint3f begin, const HPoint3f end) const
+  std::fill(m_Data.begin(), m_Data.end(), t); // warning... stl function //
-{
-    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;
 }
 
-
+} // namespace uLib
-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 //
-}
-
-}
 
 #endif // VOXIMAGE_H

src/Math/VoxImageFilter.h

@@ -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() {
+  template <class Self> void check_structural() {
-        uLibCheckFunction(Self,operator(),float,float);
+    uLibCheckFunction(Self, operator(), float, float);
-        uLibCheckFunction(Self,operator(),float,const Vector3f&);
+    uLibCheckFunction(Self, operator(), float, const Vector3f &);
-    }
+  }
 };
-}
+} // namespace Interface
-
-template < typename VoxelT > class Kernel;
 
+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>
-template < typename VoxelT, typename AlgorithmT >
+class VoxImageFilter : public Abstract::VoxImageFilter {
-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 >
+  template <class ShapeT> void SetKernelSpherical(ShapeT shape);
-    void SetKernelSpherical( ShapeT shape );
 
-    void SetKernelWeightFunction(float (*shape)(const Vector3f &));
+  void SetKernelWeightFunction(float (*shape)(const Vector3f &));
 
-    template < class ShapeT >
+  template <class ShapeT> void SetKernelWeightFunction(ShapeT shape);
-    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;
-    // protected members for algorithm access //
+  VoxImage<VoxelT> *m_Image;
-    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 "VoxImageFilter2ndStat.hpp"
-#include "VoxImageFilterThreshold.hpp"
-#include "VoxImageFilterMedian.hpp"
 #include "VoxImageFilterABTrim.hpp"
 #include "VoxImageFilterBilateral.hpp"
-#include "VoxImageFilter2ndStat.hpp"
 #include "VoxImageFilterCustom.hpp"
-
+#include "VoxImageFilterLinear.hpp"
+#include "VoxImageFilterMedian.hpp"
+#include "VoxImageFilterThreshold.hpp"

src/Math/VoxImageFilter.hpp

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

src/Math/VoxImageFilter2ndStat.hpp

@@ -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 :
+class VoxFilterAlgorithm2ndStat
-        public VoxImageFilter<VoxelT, VoxFilterAlgorithm2ndStat<VoxelT> > {
+    : public VoxImageFilter<VoxelT, VoxFilterAlgorithm2ndStat<VoxelT>> {
 
 public:
-    typedef VoxImageFilter<VoxelT, VoxFilterAlgorithm2ndStat<VoxelT> > BaseClass;
+  typedef VoxImageFilter<VoxelT, VoxFilterAlgorithm2ndStat<VoxelT>> BaseClass;
-    VoxFilterAlgorithm2ndStat(const Vector3i &size) :
+  VoxFilterAlgorithm2ndStat(const Vector3i &size) : BaseClass(size) {}
-        BaseClass(size)
-    { }
 
-    float Evaluate(const VoxImage<VoxelT> &buffer, int index)
+  float Evaluate(const VoxImage<VoxelT> &buffer, int index) {
-    {
+    const DataAllocator<VoxelT> &vbuf = buffer.ConstData();
-        const std::vector<VoxelT> &vbuf = buffer.ConstData();
+    const DataAllocator<VoxelT> &vker = this->m_KernelData.ConstData();
-        const std::vector<VoxelT> &vker = this->m_KernelData.ConstData();
+    int vox_size = vbuf.size();
-        int vox_size = vbuf.size();
+    int ker_size = vker.size();
-        int ker_size = vker.size();
+    int pos;
-        int pos;
 
-        // mean //
+    // mean //
-        float conv = 0, ksum = 0;
+    float conv = 0, ksum = 0;
-        for (int ik = 0; ik < ker_size; ++ik) {
+    for (int ik = 0; ik < ker_size; ++ik) {
-            pos = index + vker[ik].Count - vker[this->m_KernelData.GetCenterData()].Count;
+      pos = index + vker[ik].Count -
-            pos = (pos + vox_size) % vox_size;
+            vker[this->m_KernelData.GetCenterData()].Count;
-            conv += vbuf[pos].Value * vker[ik].Value;
+      pos = (pos + vox_size) % vox_size;
-            ksum += vker[ik].Value;
+      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) ;
     }
+    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

src/Math/VoxImageFilterABTrim.hpp

@@ -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 :
+__global__ void ABTrimFilterKernel(const VoxelT *in, VoxelT *out,
-        public VoxImageFilter<VoxelT, VoxFilterAlgorithmAbtrim<VoxelT> > {
+                                   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
+    for (int i = 0; i < ker_size; ++i) {
-    {
+      mfh[i].Count = i;
-        bool operator()(const VoxelT& e1, const VoxelT& e2)
+    }
-        { return e1.Value < e2.Value; }
+
-    };
+    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;
+  typedef VoxImageFilter<VoxelT, VoxFilterAlgorithmAbtrim<VoxelT>> BaseClass;
-    VoxFilterAlgorithmAbtrim(const Vector3i &size) :
+  VoxFilterAlgorithmAbtrim(const Vector3i &size) : BaseClass(size) {
-        BaseClass(size)
+    mAtrim = 0;
-    {
+    mBtrim = 0;
-        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)
+    std::sort(mfh.begin(), mfh.end(), KernelSortAscending());
-    {
+    float ker_sum = 0;
-        const std::vector<VoxelT> &vbuf = buffer.ConstData();
+    float fconv = 0;
-        const std::vector<VoxelT> &vker = this->m_KernelData.ConstData();
+    for (int ik = 0; ik < mAtrim; ik++)
-        int vox_size = vbuf.size();
+      ker_sum += vker[mfh[ik].Count].Value;
-        int ker_size = vker.size();
+    for (int ik = mAtrim; ik < ker_size - mBtrim; ik++) {
-        int pos;
+      fconv += mfh[ik].Value * vker[mfh[ik].Count].Value; // convloution //
-
+      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;
-        }
-
-        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;
     }
+    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 mAtrim;
-    int mBtrim;
+  int mBtrim;
 };
 
-
-
 ////////////////////////////////////////////////////////////////////////////////
 ////////////////////////////////////////////////////////////////////////////////
 ////////////////////////////////////////////////////////////////////////////////
 // Roberspierre Filter //
 
-
-
 template <typename VoxelT>
-class VoxFilterAlgorithmSPR :
+class VoxFilterAlgorithmSPR
-        public VoxImageFilter<VoxelT, VoxFilterAlgorithmSPR<VoxelT> > {
+    : public VoxImageFilter<VoxelT, VoxFilterAlgorithmSPR<VoxelT>> {
 
-    struct KernelSortAscending
+  struct KernelSortAscending {
-    {
+    bool operator()(const VoxelT &e1, const VoxelT &e2) {
-        bool operator()(const VoxelT& e1, const VoxelT& e2)
+      return e1.Value < e2.Value;
-        { return e1.Value < e2.Value; }
+    }
-    };
+  };
 
 public:
-    typedef VoxImageFilter<VoxelT, VoxFilterAlgorithmSPR<VoxelT> > BaseClass;
+  typedef VoxImageFilter<VoxelT, VoxFilterAlgorithmSPR<VoxelT>> BaseClass;
-    VoxFilterAlgorithmSPR(const Vector3i &size) :
+  VoxFilterAlgorithmSPR(const Vector3i &size) : BaseClass(size) {
-        BaseClass(size)
+    mAtrim = 0;
-    {
+    mBtrim = 0;
-        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)
+    std::sort(mfh.begin(), mfh.end(), KernelSortAscending());
-    {
+    float spr = vbuf[index].Value;
-        const std::vector<VoxelT> &vbuf = buffer.ConstData();
+    if ((mAtrim > 0 && spr <= mfh[mAtrim - 1].Value) ||
-        const std::vector<VoxelT> &vker = this->m_KernelData.ConstData();
+        (mBtrim > 0 && spr >= mfh[ker_size - mBtrim].Value)) {
-        int vox_size = vbuf.size();
+      float ker_sum = 0;
-        int ker_size = vker.size();
+      float fconv = 0;
-        int pos;
+      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);
+      return fconv / ker_sum;
-        for (int i = 0; i < ker_size; ++i)
+    } else
-            mfh[i].Count = i; //index key for ordering function
+      return spr;
-        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());
+  inline void SetABTrim(int a, int b) {
-        float spr = vbuf[index].Value;
+    mAtrim = a;
-        if( (mAtrim > 0 && spr <= mfh[mAtrim-1].Value) ||
+    mBtrim = b;
-                (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; }
 
 private:
-    int mAtrim;
+  int mAtrim;
-    int mBtrim;
+  int mBtrim;
 };
 
-
+} // namespace uLib
-
-
-}
 
 #endif // VOXIMAGEFILTERABTRIM_HPP

src/Math/VoxImageFilterBilateral.hpp

@@ -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 :
+class VoxFilterAlgorithmBilateral
-        public VoxImageFilter<VoxelT, VoxFilterAlgorithmBilateral<VoxelT> > {
+    : public VoxImageFilter<VoxelT, VoxFilterAlgorithmBilateral<VoxelT>> {
 public:
-    typedef VoxImageFilter<VoxelT, VoxFilterAlgorithmBilateral<VoxelT> > BaseClass;
+  typedef VoxImageFilter<VoxelT, VoxFilterAlgorithmBilateral<VoxelT>> BaseClass;
-    VoxFilterAlgorithmBilateral(const Vector3i &size) : BaseClass(size) {
+  VoxFilterAlgorithmBilateral(const Vector3i &size) : BaseClass(size) {
-        m_sigma = 1;
+    m_sigma = 1;
-    }
+  }
 
-    float Evaluate(const VoxImage<VoxelT> &buffer, int index)
+  float Evaluate(const VoxImage<VoxelT> &buffer, int index) {
-    {
+    const DataAllocator<VoxelT> &vbuf = buffer.ConstData();
-        const std::vector<VoxelT> &vbuf = buffer.ConstData();
+    const DataAllocator<VoxelT> &vker = this->m_KernelData.ConstData();
-        const std::vector<VoxelT> &vker = this->m_KernelData.ConstData();
+    int vox_size = vbuf.size();
-        int vox_size = vbuf.size();
+    int ker_size = vker.size();
-        int ker_size = vker.size();
+    int pos;
-        int pos;
+    float conv = 0, ksum = 0;
-        float conv = 0, ksum = 0;
+    float gamma_smooth;
-        float gamma_smooth;
+    for (int ik = 0; ik < ker_size; ++ik) {
-        for (int ik = 0; ik < ker_size; ++ik) {
+      //            if (ik==this->m_KernelData.GetCenterData()) continue;
-            //            if (ik==this->m_KernelData.GetCenterData()) continue;
+      pos = index + vker[ik].Count -
-            pos = index + vker[ik].Count - vker[this->m_KernelData.GetCenterData()].Count;
+            vker[this->m_KernelData.GetCenterData()].Count;
-            pos = (pos + vox_size) % vox_size;
+      pos = (pos + vox_size) % vox_size;
-            gamma_smooth = compute_gauss( fabs(vbuf[index].Value - vbuf[pos].Value) * 1.E6 );
+      gamma_smooth =
-            conv += vbuf[pos].Value * vker[ik].Value * gamma_smooth;
+          compute_gauss(fabs(vbuf[index].Value - vbuf[pos].Value) * 1.E6);
-            ksum += vker[ik].Value * gamma_smooth;
+      conv += vbuf[pos].Value * vker[ik].Value * gamma_smooth;
-        }
+      ksum += vker[ik].Value * gamma_smooth;
-        return conv / ksum;
     }
+    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) {
+  inline float compute_gauss(const float x) {
-        return 1/(sqrt(2*M_PI)* m_sigma) * exp(-0.5*(x*x)/(m_sigma*m_sigma));
+    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 :
+class VoxFilterAlgorithmBilateralTrim
-        public VoxImageFilter<VoxelT, VoxFilterAlgorithmBilateralTrim<VoxelT> > {
+    : public VoxImageFilter<VoxelT, VoxFilterAlgorithmBilateralTrim<VoxelT>> {
 
-    typedef std::pair<float,float> FPair;
+  typedef std::pair<float, float> FPair;
 
-    struct KernelSortAscending
+  struct KernelSortAscending {
-    {
+    bool operator()(const FPair &e1, const FPair &e2) {
-        bool operator()(const FPair& e1, const FPair& e2)
+      return e1.second < e2.second;
-        { return e1.second < e2.second; }
+    }
-    };
+  };
 
 public:
-    typedef VoxImageFilter<VoxelT, VoxFilterAlgorithmBilateralTrim<VoxelT> > BaseClass;
+  typedef VoxImageFilter<VoxelT, VoxFilterAlgorithmBilateralTrim<VoxelT>>
-    VoxFilterAlgorithmBilateralTrim(const Vector3i &size) : BaseClass(size) {
+      BaseClass;
-        m_sigma = 1;
+  VoxFilterAlgorithmBilateralTrim(const Vector3i &size) : BaseClass(size) {
-        mAtrim = 0;
+    m_sigma = 1;
-        mBtrim = 0;
+    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)
+    float conv = 0, ksum = 0;
-    {
+    float gamma_smooth;
-        const std::vector<VoxelT> &vbuf = buffer.ConstData();
+    //        for (int ik = 0; ik < mAtrim; ik++)
-        const std::vector<VoxelT> &vker = this->m_KernelData.ConstData();
+    //            ksum += mfh[ik].first;
-        int img_size = vbuf.size();
+    for (int ik = mAtrim; ik < ker_size - mBtrim; ik++) {
-        int ker_size = vker.size();
+      gamma_smooth =
-        int pos;
+          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;
-
-        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;
     }
+    //        for (int ik = ker_size - mBtrim; ik < ker_size; ik++)
+    //            ksum += mfh[ik].first;
 
-    inline void SetIntensitySigma(const float s) { m_sigma = s; }
+    return conv / ksum;
-    inline void SetABTrim(int a, int b) { mAtrim = a; mBtrim = b; }
+  }
+
+  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) {
+  inline float compute_gauss(const float x) {
-        return 1/(sqrt(2*M_PI)* m_sigma) * exp(-0.5*(x*x)/(m_sigma*m_sigma));
+    return 1 / (sqrt(2 * M_PI) * m_sigma) *
-    }
+           exp(-0.5 * (x * x) / (m_sigma * m_sigma));
+  }
 
-    Scalarf m_sigma;
+  Scalarf m_sigma;
-    int mAtrim;
+  int mAtrim;
-    int mBtrim;
+  int mBtrim;
 };
 
-}
+} // namespace uLib
 
 #endif // VOXIMAGEFILTERBILATERAL_HPP

src/Math/VoxImageFilterCustom.hpp

@@ -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 :
+class VoxFilterAlgorithmCustom
-        public VoxImageFilter<VoxelT, VoxFilterAlgorithmCustom<VoxelT> > {
+    : 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;
+  typedef VoxImageFilter<VoxelT, VoxFilterAlgorithmCustom<VoxelT>> BaseClass;
-    VoxFilterAlgorithmCustom(const Vector3i &size) :
+  VoxFilterAlgorithmCustom(const Vector3i &size)
-        BaseClass(size), m_CustomEvaluate(NULL)
+      : BaseClass(size), m_CustomEvaluate(NULL) {}
-    {}
 
-    float Evaluate(const VoxImage<VoxelT> &buffer, int index)
+  float Evaluate(const VoxImage<VoxelT> &buffer, int index) {
-    {
+    if (likely(m_CustomEvaluate)) {
-        if(likely(m_CustomEvaluate)) {
+      const DataAllocator<VoxelT> &vbuf = buffer.ConstData();
-            const std::vector<VoxelT> &vbuf = buffer.ConstData();
+      const DataAllocator<VoxelT> &vker = this->m_KernelData.ConstData();
-            const std::vector<VoxelT> &vker = this->m_KernelData.ConstData();
+      int vox_size = vbuf.size();
-            int vox_size = vbuf.size();
+      int ker_size = vker.size();
-            int ker_size = vker.size();
+      int pos;
-            int pos;
 
-            float ker_sum = 0;
+      float ker_sum = 0;
-            std::vector<Scalarf> mfh(ker_size);
+      std::vector<Scalarf> mfh(ker_size);
-            for (int ik = 0; ik < ker_size; ik++) {
+      for (int ik = 0; ik < ker_size; ik++) {
-                pos = index + vker[ik].Count - vker[this->m_KernelData.GetCenterData()].Count;
+        pos = index + vker[ik].Count -
-                pos = (pos + vox_size) % vox_size;
+              vker[this->m_KernelData.GetCenterData()].Count;
-                mfh[ik] = vbuf[pos].Value * vker[ik].Value;
+        pos = (pos + vox_size) % vox_size;
-                ker_sum += vker[ik].Value;
+        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 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

src/Math/VoxImageFilterLinear.hpp

@@ -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 :
+class VoxFilterAlgorithmLinear
-        public VoxImageFilter<VoxelT, VoxFilterAlgorithmLinear<VoxelT> > {
+    : public VoxImageFilter<VoxelT, VoxFilterAlgorithmLinear<VoxelT>> {
 public:
-    typedef VoxImageFilter<VoxelT, VoxFilterAlgorithmLinear<VoxelT> > BaseClass;
+  typedef VoxImageFilter<VoxelT, VoxFilterAlgorithmLinear<VoxelT>> BaseClass;
-    VoxFilterAlgorithmLinear(const Vector3i &size) : BaseClass(size) {}
+  VoxFilterAlgorithmLinear(const Vector3i &size) : BaseClass(size) {}
 
-    float Evaluate(const VoxImage<VoxelT> &buffer, int index)
+#ifdef USE_CUDA
-    {
+  void Run() {
-        const std::vector<VoxelT> &vbuf = buffer.ConstData();
+    if (this->m_Image->Data().GetDevice() == MemoryDevice::VRAM ||
-        const std::vector<VoxelT> &vker = this->m_KernelData.ConstData();
+        this->m_KernelData.Data().GetDevice() == MemoryDevice::VRAM) {
-        int vox_size = vbuf.size();
+
-        int ker_size = vker.size();
+      this->m_Image->Data().MoveToVRAM();
-        int pos;
+      this->m_KernelData.Data().MoveToVRAM();
-        float conv = 0, ksum = 0;
+
-        for (int ik = 0; ik < ker_size; ++ik) {
+      VoxImage<VoxelT> buffer = *(this->m_Image);
-            pos = index + vker[ik].Count - vker[this->m_KernelData.GetCenterData()].Count;
+      buffer.Data().MoveToVRAM();
-            pos = (pos + vox_size) % vox_size;
+
-            conv += vbuf[pos].Value * vker[ik].Value;
+      int vox_size = buffer.Data().size();
-            ksum += vker[ik].Value;
+      int ker_size = this->m_KernelData.Data().size();
-        }
+
-        return conv / ksum;
+      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

src/Math/VoxImageFilterMedian.hpp

@@ -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 :
+class VoxFilterAlgorithmMedian
-        public VoxImageFilter<VoxelT, VoxFilterAlgorithmMedian<VoxelT> > {
+    : public VoxImageFilter<VoxelT, VoxFilterAlgorithmMedian<VoxelT>> {
 public:
-    typedef VoxImageFilter<VoxelT, VoxFilterAlgorithmMedian<VoxelT> > BaseClass;
+  typedef VoxImageFilter<VoxelT, VoxFilterAlgorithmMedian<VoxelT>> BaseClass;
-    VoxFilterAlgorithmMedian(const Vector3i &size) : BaseClass(size) {}
+  VoxFilterAlgorithmMedian(const Vector3i &size) : BaseClass(size) {}
 
-    float Evaluate(const VoxImage<VoxelT> &buffer, int index)
+  float Evaluate(const VoxImage<VoxelT> &buffer, int index) {
-    {
+    const DataAllocator<VoxelT> &vbuf = buffer.ConstData();
-        const std::vector<VoxelT> &vbuf = buffer.ConstData();
+    const DataAllocator<VoxelT> &vker = this->m_KernelData.ConstData();
-        const std::vector<VoxelT> &vker = this->m_KernelData.ConstData();
+    int vox_size = vbuf.size();
-        int vox_size = vbuf.size();
+    int ker_size = vker.size();
-        int ker_size = vker.size();
+    int pos;
-        int pos;
 
-        std::vector<float> mfh(ker_size);
+    std::vector<float> mfh(ker_size);
-        for (int ik = 0; ik < ker_size; ik++) {
+    for (int ik = 0; ik < ker_size; ik++) {
-            pos = index + vker[ik].Count - vker[this->m_KernelData.GetCenterData()].Count;
+      pos = index + vker[ik].Count -
-            pos = (pos + vox_size) % vox_size;
+            vker[this->m_KernelData.GetCenterData()].Count;
-            mfh[ik] = vbuf[pos].Value * vker[ik].Value;
+      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::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

src/Math/testing/CMakeLists.txt

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

src/Math/testing/VoxImageFilterTest.cpp

@@ -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;
+  Scalarf Value;
-    unsigned int Count;
+  unsigned int Count;
 };
 
-float GaussianShape(float d)
+float GaussianShape(float d) {
-{
+  // normalized manually .. fix //
-    // normalized manually .. fix //
+  return 4.5 * exp(-d * 4.5);
-    return 4.5 * exp(-d * 4.5);
 }
 
-
 class GaussianShapeClass : public Interface::VoxImageFilterShape {
 public:
-    GaussianShapeClass(float sigma) :
+  GaussianShapeClass(float sigma) : m_sigma(sigma) {}
-        m_sigma(sigma)
-    {}
 
-    float operator ()(float d) {
+  float operator()(float d) { return (1 / m_sigma) * exp(-d / m_sigma); }
-        return (1/m_sigma) * exp(-d/m_sigma);
-    }
 
 private:
-    float m_sigma;
+  float m_sigma;
 };
 
-
+static float MaxInVector(const std::vector<float> &v) {
-static float MaxInVector(const std::vector<float> &v)
+  float max = 0;
-{
+  for (int i = 0; i < v.size(); ++i)
-    float max = 0;
+    if (v.at(i) > max)
-    for(int i=0; i<v.size(); ++i)
+      max = v.at(i);
-        if(v.at(i) > max) max = v.at(i);
+  return max;
-    return max;
 }
 
+int main() {
+  BEGIN_TESTING(VoxImageFilters);
 
-int main()
+  VoxImage<TestVoxel> image(Vector3i(20, 30, 40));
-{
+  image[Vector3i(10, 10, 10)].Value = 1;
-    BEGIN_TESTING(VoxImageFilters);
+  // 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;
 
-    ////////////////////////////////////////////////////////////////////////////
+    std::vector<float> values;
-    ////////////////////////////////////////////////////////////////////////////
+    for (int i = 0; i < filter.GetKernelData().GetDims().prod(); ++i) {
-    ////////////////////////////////////////////////////////////////////////////
+      values.push_back(1.);
-    // RPS //
+      std::cout << values[i] << " ";
+    }
+    std::cout << "\n";
 
-    {
+    filter.SetImage(&filtered);
-        VoxFilterAlgorithmSPR<TestVoxel> filter(Vector3i(2,3,4));
 
-        VoxImage<TestVoxel> filtered = image;
+    filter.SetKernelNumericXZY(values);
 
-        std::vector<float> values;
+    filter.SetABTrim(0, 2);
-        for(int i=0; i < filter.GetKernelData().GetDims().prod(); ++i) {
-            values.push_back(1.);
-            std::cout << values[i] << " ";
-        }
-        std::cout << "\n";
 
-        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;
+  // CUDA Allocator Transfer Test //
-        image[Vector3i(9,10,8)].Value = 2;
+  {
-        image.ExportToVtk("test_filter_max_original.vtk",0);
+    VoxImage<TestVoxel> image(Vector3i(10, 10, 10));
+    image[Vector3i(5, 5, 5)].Value = 1;
 
+    VoxFilterAlgorithmLinear<TestVoxel> filter(Vector3i(3, 3, 3));
 
-
+    std::vector<float> values;
-        VoxFilterAlgorithmCustom<TestVoxel> filter(Vector3i(3,3,4));
+    for (int i = 0; i < filter.GetKernelData().GetDims().prod(); ++i) {
-
+      values.push_back(1.0f);
-        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);
     }
 
+    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;
 }