#include <pybind11/eigen.h>
#include <pybind11/pybind11.h>
#include <pybind11/stl.h>
#include <pybind11/stl_bind.h>

#include <pybind11/numpy.h>

#include "Math/Accumulator.h"
#include "Math/ContainerBox.h"
#include "Math/Dense.h"
#include "Math/Geometry.h"
#include "Math/Structured2DGrid.h"
#include "Math/Structured4DGrid.h"
#include "Math/StructuredData.h"
#include "Math/StructuredGrid.h"
#include "Math/Transform.h"
#include "Math/TriangleMesh.h"
#include "Math/QuadMesh.h"
#include "Math/VoxImage.h"
#include "Math/VoxRaytracer.h"

namespace py = pybind11;
using namespace uLib;

PYBIND11_MAKE_OPAQUE(uLib::Vector<Scalari>);
PYBIND11_MAKE_OPAQUE(uLib::Vector<Scalarui>);
PYBIND11_MAKE_OPAQUE(uLib::Vector<Scalarl>);
PYBIND11_MAKE_OPAQUE(uLib::Vector<Scalarul>);
PYBIND11_MAKE_OPAQUE(uLib::Vector<Scalarf>);
PYBIND11_MAKE_OPAQUE(uLib::Vector<Scalard>);

PYBIND11_MAKE_OPAQUE(uLib::Vector<Vector3f>);
PYBIND11_MAKE_OPAQUE(uLib::Vector<Vector3i>);
PYBIND11_MAKE_OPAQUE(uLib::Vector<Vector4f>);
PYBIND11_MAKE_OPAQUE(uLib::Vector<Vector4i>);
PYBIND11_MAKE_OPAQUE(uLib::Vector<Vector3d>);
PYBIND11_MAKE_OPAQUE(uLib::Vector<Vector4d>);

PYBIND11_MAKE_OPAQUE(uLib::Vector<Voxel>);
PYBIND11_MAKE_OPAQUE(uLib::Vector<VoxRaytracer::RayData::Element>);

template <typename MatrixType>
void bind_eigen_type(py::module_ &m, const char *name) {
  using Scalar = typename MatrixType::Scalar;
  constexpr bool is_vector = MatrixType::IsVectorAtCompileTime;

  // Default constructor (zeros)
  m.def(name, []() -> MatrixType {
    if constexpr (MatrixType::RowsAtCompileTime == Eigen::Dynamic ||
                  MatrixType::ColsAtCompileTime == Eigen::Dynamic) {
      return MatrixType(); // Empty dynamic matrix
    } else {
      return MatrixType::Zero(); // Zero static matrix
    }
  });

  // Specialized constructor for dynamic matrices
  if constexpr (MatrixType::RowsAtCompileTime == Eigen::Dynamic ||
                MatrixType::ColsAtCompileTime == Eigen::Dynamic) {
    m.def(name, [](int rows, int cols) -> MatrixType {
      MatrixType mat;
      mat.setZero(rows, cols);
      return mat;
    });
  }

  // Initialize from list
  m.def(name, [](py::list l) -> MatrixType {
    MatrixType mat;
    if constexpr (is_vector) {
      mat.setZero(l.size());
      for (size_t i = 0; i < l.size(); ++i) {
        mat(i) = l[i].cast<Scalar>();
      }
    } else {
      int rows = MatrixType::RowsAtCompileTime == Eigen::Dynamic
                     ? (int)std::sqrt(l.size())
                     : MatrixType::RowsAtCompileTime;
      int cols = MatrixType::ColsAtCompileTime == Eigen::Dynamic
                     ? (int)std::sqrt(l.size())
                     : MatrixType::ColsAtCompileTime;
      mat.setZero(rows, cols);
      for (size_t i = 0; i < (size_t)l.size(); ++i) {
        mat(i / cols, i % cols) = l[i].cast<Scalar>();
      }
    }
    return mat;
  });

  // Initialize from py::array
  m.def(name,
        [](py::array_t<Scalar, py::array::c_style | py::array::forcecast> arr)
            -> MatrixType {
          auto buf = arr.request();
          MatrixType mat;
          if constexpr (is_vector) {
            mat.setZero(buf.size);
            Scalar *ptr = static_cast<Scalar *>(buf.ptr);
            for (ssize_t i = 0; i < buf.size; ++i)
              mat(i) = ptr[i];
          } else {
            int rows = buf.shape.size() > 0 ? (int)buf.shape[0] : 1;
            int cols = buf.shape.size() > 1 ? (int)buf.shape[1] : 1;
            mat.setZero(rows, cols);
            Scalar *ptr = static_cast<Scalar *>(buf.ptr);
            for (int i = 0; i < rows; ++i) {
              for (int j = 0; j < cols; ++j) {
                mat(i, j) = ptr[i * cols + j];
              }
            }
          }
          return mat;
        });
}

void init_math(py::module_ &m) {

  // 1. Basic Eigen Types (Vectors and Matrices)
  bind_eigen_type<Vector1f>(m, "Vector1f");
  bind_eigen_type<Vector2f>(m, "Vector2f");
  bind_eigen_type<Vector3f>(m, "Vector3f");
  bind_eigen_type<Vector4f>(m, "Vector4f");
  bind_eigen_type<Vector1i>(m, "Vector1i");
  bind_eigen_type<Vector2i>(m, "Vector2i");
  bind_eigen_type<Vector3i>(m, "Vector3i");
  bind_eigen_type<Vector4i>(m, "Vector4i");
  bind_eigen_type<Vector1d>(m, "Vector1d");
  bind_eigen_type<Vector2d>(m, "Vector2d");
  bind_eigen_type<Vector3d>(m, "Vector3d");
  bind_eigen_type<Vector4d>(m, "Vector4d");

  bind_eigen_type<Matrix2f>(m, "Matrix2f");
  bind_eigen_type<Matrix3f>(m, "Matrix3f");
  bind_eigen_type<Matrix4f>(m, "Matrix4f");
  bind_eigen_type<Matrix2i>(m, "Matrix2i");
  bind_eigen_type<Matrix3i>(m, "Matrix3i");
  bind_eigen_type<Matrix4i>(m, "Matrix4i");
  bind_eigen_type<Matrix2d>(m, "Matrix2d");
  bind_eigen_type<Matrix3d>(m, "Matrix3d");
  bind_eigen_type<Matrix4d>(m, "Matrix4d");

  bind_eigen_type<MatrixXi>(m, "MatrixXi");
  bind_eigen_type<MatrixXf>(m, "MatrixXf");
  bind_eigen_type<MatrixXd>(m, "MatrixXd");

  // 2. Homogeneous types
  py::class_<HPoint3f>(m, "HPoint3f")
      .def(py::init<>())
      .def(py::init<float, float, float>())
      .def(py::init<Vector3f &>());
  py::class_<HVector3f>(m, "HVector3f")
      .def(py::init<>())
      .def(py::init<float, float, float>())
      .def(py::init<Vector3f &>());
  py::class_<HLine3f>(m, "HLine3f")
      .def(py::init<>())
      .def_readwrite("origin", &HLine3f::origin)
      .def_readwrite("direction", &HLine3f::direction);
  py::class_<HError3f>(m, "HError3f")
      .def(py::init<>())
      .def_readwrite("position_error", &HError3f::position_error)
      .def_readwrite("direction_error", &HError3f::direction_error);

  // 3. Dynamic Vectors (uLib::Vector)
  py::bind_vector<uLib::Vector<Scalari>>(m, "Vector_i")
      .def("MoveToVRAM", &uLib::Vector<Scalari>::MoveToVRAM)
      .def("MoveToRAM", &uLib::Vector<Scalari>::MoveToRAM);
  py::bind_vector<uLib::Vector<Scalarui>>(m, "Vector_ui")
      .def("MoveToVRAM", &uLib::Vector<Scalarui>::MoveToVRAM)
      .def("MoveToRAM", &uLib::Vector<Scalarui>::MoveToRAM);
  py::bind_vector<uLib::Vector<Scalarl>>(m, "Vector_l")
      .def("MoveToVRAM", &uLib::Vector<Scalarl>::MoveToVRAM)
      .def("MoveToRAM", &uLib::Vector<Scalarl>::MoveToRAM);
  py::bind_vector<uLib::Vector<Scalarul>>(m, "Vector_ul")
      .def("MoveToVRAM", &uLib::Vector<Scalarul>::MoveToVRAM)
      .def("MoveToRAM", &uLib::Vector<Scalarul>::MoveToRAM);
  py::bind_vector<uLib::Vector<Scalarf>>(m, "Vector_f")
      .def("MoveToVRAM", &uLib::Vector<Scalarf>::MoveToVRAM)
      .def("MoveToRAM", &uLib::Vector<Scalarf>::MoveToRAM);
  py::bind_vector<uLib::Vector<Scalard>>(m, "Vector_d")
      .def("MoveToVRAM", &uLib::Vector<Scalard>::MoveToVRAM)
      .def("MoveToRAM", &uLib::Vector<Scalard>::MoveToRAM);

  py::bind_vector<uLib::Vector<Vector3f>>(m, "Vector_Vector3f")
      .def("MoveToVRAM", &uLib::Vector<Vector3f>::MoveToVRAM)
      .def("MoveToRAM", &uLib::Vector<Vector3f>::MoveToRAM);
  py::bind_vector<uLib::Vector<Vector3i>>(m, "Vector_Vector3i")
      .def("MoveToVRAM", &uLib::Vector<Vector3i>::MoveToVRAM)
      .def("MoveToRAM", &uLib::Vector<Vector3i>::MoveToRAM);
  py::bind_vector<uLib::Vector<Vector4f>>(m, "Vector_Vector4f")
      .def("MoveToVRAM", &uLib::Vector<Vector4f>::MoveToVRAM)
      .def("MoveToRAM", &uLib::Vector<Vector4f>::MoveToRAM);
  py::bind_vector<uLib::Vector<Vector4i>>(m, "Vector_Vector4i")
      .def("MoveToVRAM", &uLib::Vector<Vector4i>::MoveToVRAM)
      .def("MoveToRAM", &uLib::Vector<Vector4i>::MoveToRAM);
  py::bind_vector<uLib::Vector<Vector3d>>(m, "Vector_Vector3d")
      .def("MoveToVRAM", &uLib::Vector<Vector3d>::MoveToVRAM)
      .def("MoveToRAM", &uLib::Vector<Vector3d>::MoveToRAM);
  py::bind_vector<uLib::Vector<Vector4d>>(m, "Vector_Vector4d")
      .def("MoveToVRAM", &uLib::Vector<Vector4d>::MoveToVRAM)
      .def("MoveToRAM", &uLib::Vector<Vector4d>::MoveToRAM);
  py::bind_vector<uLib::Vector<Voxel>>(m, "Vector_Voxel")
      .def("MoveToVRAM", &uLib::Vector<Voxel>::MoveToVRAM)
      .def("MoveToRAM", &uLib::Vector<Voxel>::MoveToRAM);
  py::bind_vector<uLib::Vector<VoxRaytracer::RayData::Element>>(
      m, "Vector_VoxRaytracerRayDataElement")
      .def("MoveToVRAM",
           &uLib::Vector<VoxRaytracer::RayData::Element>::MoveToVRAM)
      .def("MoveToRAM",
           &uLib::Vector<VoxRaytracer::RayData::Element>::MoveToRAM);

  // 4. Accumulators
  py::class_<Accumulator_Mean<float>>(m, "Accumulator_Mean_f")
      .def(py::init<>())
      .def("AddPass", &Accumulator_Mean<float>::AddPass)
      .def("__call__",
           py::overload_cast<const float>(&Accumulator_Mean<float>::operator()))
      .def("__call__", py::overload_cast<>(&Accumulator_Mean<float>::operator(),
                                           py::const_));

  py::class_<Accumulator_Mean<double>>(m, "Accumulator_Mean_d")
      .def(py::init<>())
      .def("AddPass", &Accumulator_Mean<double>::AddPass)
      .def("__call__", py::overload_cast<const double>(
                           &Accumulator_Mean<double>::operator()))
      .def("__call__", py::overload_cast<>(
                           &Accumulator_Mean<double>::operator(), py::const_));

  py::class_<Accumulator_ABTrim<float>>(m, "Accumulator_ABTrim_f")
      .def(py::init<>())
      .def("SetABTrim", &Accumulator_ABTrim<float>::SetABTrim)
      .def("__iadd__",
           [](Accumulator_ABTrim<float> &self, float val) {
             self += val;
             return &self;
           })
      .def("__call__", &Accumulator_ABTrim<float>::operator());

  py::class_<Accumulator_ABTrim<double>>(m, "Accumulator_ABTrim_d")
      .def(py::init<>())
      .def("SetABTrim", &Accumulator_ABTrim<double>::SetABTrim)
      .def("__iadd__",
           [](Accumulator_ABTrim<double> &self, double val) {
             self += val;
             return &self;
           })
      .def("__call__", &Accumulator_ABTrim<double>::operator());

  py::class_<Accumulator_ABClip<float>>(m, "Accumulator_ABClip_f")
      .def(py::init<>())
      .def("SetABTrim", &Accumulator_ABClip<float>::SetABTrim)
      .def("__iadd__",
           [](Accumulator_ABClip<float> &self, float val) {
             self += val;
             return &self;
           })
      .def("__call__", &Accumulator_ABClip<float>::operator());

  py::class_<Accumulator_ABClip<double>>(m, "Accumulator_ABClip_d")
      .def(py::init<>())
      .def("SetABTrim", &Accumulator_ABClip<double>::SetABTrim)
      .def("__iadd__",
           [](Accumulator_ABClip<double> &self, double val) {
             self += val;
             return &self;
           })
      .def("__call__", &Accumulator_ABClip<double>::operator());

  // 5. Core Math Structures
  py::class_<AffineTransform, std::shared_ptr<AffineTransform>>(m,
                                                              "AffineTransform")
      .def(py::init<>())
      .def("GetWorldMatrix", &AffineTransform::GetWorldMatrix)
      .def("SetPosition", &AffineTransform::SetPosition)
      .def("GetPosition", &AffineTransform::GetPosition)
      .def("Translate", &AffineTransform::Translate)
      .def("Scale", &AffineTransform::Scale)
      .def("SetRotation", &AffineTransform::SetRotation)
      .def("GetRotation", &AffineTransform::GetRotation)
      .def("Rotate",
           py::overload_cast<const Matrix3f>(&AffineTransform::Rotate))
      .def("Rotate",
           py::overload_cast<float, Vector3f>(&AffineTransform::Rotate))
      .def("Rotate", py::overload_cast<Vector3f>(&AffineTransform::Rotate))
      .def("EulerYZYRotate", &AffineTransform::EulerYZYRotate)
      .def("FlipAxes", &AffineTransform::FlipAxes);

  py::class_<Geometry, AffineTransform, std::shared_ptr<Geometry>>(m, "Geometry")
      .def(py::init<>())
      .def("GetWorldPoint", py::overload_cast<const Vector4f>(
                                &Geometry::GetWorldPoint, py::const_))
      .def("GetWorldPoint",
           py::overload_cast<float, float, float>(&Geometry::GetWorldPoint))
      .def("GetLocalPoint", py::overload_cast<const Vector4f>(
                                &Geometry::GetLocalPoint, py::const_))
      .def("GetLocalPoint",
           py::overload_cast<float, float, float>(&Geometry::GetLocalPoint));

  py::class_<ContainerBox, AffineTransform, Object, std::shared_ptr<ContainerBox>>(
      m, "ContainerBox")
      .def(py::init<>())
      .def("SetOrigin", &ContainerBox::SetOrigin)
      .def("GetOrigin", &ContainerBox::GetOrigin)
      .def("SetSize", &ContainerBox::SetSize)
      .def("GetSize", &ContainerBox::GetSize)
      .def("GetLocalMatrix", &ContainerBox::GetLocalMatrix)
      .def("GetWorldMatrix", &ContainerBox::GetWorldMatrix)
      .def("GetWorldPoint", py::overload_cast<const Vector4f &>(
                                &ContainerBox::GetWorldPoint, py::const_))
      .def("GetWorldPoint",
           py::overload_cast<float, float, float>(&ContainerBox::GetWorldPoint))
      .def("GetLocalPoint", py::overload_cast<const Vector4f &>(
                                &ContainerBox::GetLocalPoint, py::const_))
      .def("GetLocalPoint",
           py::overload_cast<float, float, float>(&ContainerBox::GetLocalPoint))
      .def("FlipLocalAxes", &ContainerBox::FlipLocalAxes);

  py::enum_<StructuredData::_Order>(m, "StructuredDataOrder")
      .value("CustomOrder", StructuredData::CustomOrder)
      .value("XYZ", StructuredData::XYZ)
      .value("XZY", StructuredData::XZY)
      .value("YXZ", StructuredData::YXZ)
      .value("YZX", StructuredData::YZX)
      .value("ZXY", StructuredData::ZXY)
      .value("ZYX", StructuredData::ZYX)
      .export_values();

  py::class_<StructuredData, std::shared_ptr<StructuredData>>(m, "StructuredData")
      .def(py::init<const Vector3i &>())
      .def("GetDims", &StructuredData::GetDims)
      .def("SetDims", &StructuredData::SetDims)
      .def("GetIncrements", &StructuredData::GetIncrements)
      .def("SetIncrements", &StructuredData::SetIncrements)
      .def("SetDataOrder", &StructuredData::SetDataOrder)
      .def("GetDataOrder", &StructuredData::GetDataOrder)
      .def("IsInsideGrid", &StructuredData::IsInsideGrid)
      .def("Map", &StructuredData::Map)
      .def("UnMap", &StructuredData::UnMap);

  py::class_<StructuredGrid, ContainerBox, StructuredData,
             std::shared_ptr<StructuredGrid>>(m, "StructuredGrid")
      .def(py::init<const Vector3i &>())
      .def("SetSpacing", &StructuredGrid::SetSpacing)
      .def("GetSpacing", &StructuredGrid::GetSpacing)
      .def("IsInsideBounds", &StructuredGrid::IsInsideBounds)
      .def("Find", [](StructuredGrid &self, Vector3f pt) {
        return self.Find(HPoint3f(pt));
      });

  py::class_<Structured2DGrid>(m, "Structured2DGrid")
      .def(py::init<>())
      .def("SetDims", &Structured2DGrid::SetDims)
      .def("GetDims", &Structured2DGrid::GetDims)
      .def("IsInsideGrid", &Structured2DGrid::IsInsideGrid)
      .def("Map", &Structured2DGrid::Map)
      .def("UnMap", &Structured2DGrid::UnMap)
      .def("SetPhysicalSpace", &Structured2DGrid::SetPhysicalSpace)
      .def("GetSpacing", &Structured2DGrid::GetSpacing)
      .def("GetOrigin", &Structured2DGrid::GetOrigin)
      .def("IsInsideBounds", &Structured2DGrid::IsInsideBounds)
      .def("PhysicsToUnitSpace", &Structured2DGrid::PhysicsToUnitSpace)
      .def("UnitToPhysicsSpace", &Structured2DGrid::UnitToPhysicsSpace)
      .def("SetDebug", &Structured2DGrid::SetDebug);

  py::class_<Structured4DGrid>(m, "Structured4DGrid")
      .def(py::init<>())
      .def("SetDims", &Structured4DGrid::SetDims)
      .def("GetDims", &Structured4DGrid::GetDims)
      .def("IsInsideGrid", &Structured4DGrid::IsInsideGrid)
      .def("Map", &Structured4DGrid::Map)
      .def("UnMap", &Structured4DGrid::UnMap)
      .def("SetPhysicalSpace", &Structured4DGrid::SetPhysicalSpace)
      .def("GetSpacing", &Structured4DGrid::GetSpacing)
      .def("GetOrigin", &Structured4DGrid::GetOrigin)
      .def("IsInsideBounds", &Structured4DGrid::IsInsideBounds)
      .def("PhysicsToUnitSpace", &Structured4DGrid::PhysicsToUnitSpace)
      .def("UnitToPhysicsSpace", &Structured4DGrid::UnitToPhysicsSpace)
      .def("SetDebug", &Structured4DGrid::SetDebug);

  // 6. High-level Structures
  py::class_<Voxel>(m, "Voxel")
      .def(py::init<>())
      .def_readwrite("Value", &Voxel::Value)
      .def_readwrite("Count", &Voxel::Count);

  py::class_<Abstract::VoxImage, StructuredGrid,
             std::shared_ptr<Abstract::VoxImage>>(m, "AbstractVoxImage")
      .def("GetValue", py::overload_cast<const Vector3i &>(
                           &Abstract::VoxImage::GetValue, py::const_))
      .def("GetValue", py::overload_cast<const int>(
                           &Abstract::VoxImage::GetValue, py::const_))
      .def("SetValue", py::overload_cast<const Vector3i &, float>(
                           &Abstract::VoxImage::SetValue))
      .def("SetValue",
           py::overload_cast<const int, float>(&Abstract::VoxImage::SetValue))
      .def("ExportToVtk", &Abstract::VoxImage::ExportToVtk)
      .def("ExportToVti", &Abstract::VoxImage::ExportToVti)
      .def("ImportFromVtk", &Abstract::VoxImage::ImportFromVtk)
      .def("ImportFromVti", &Abstract::VoxImage::ImportFromVti);

  py::class_<VoxImage<Voxel>, Abstract::VoxImage,
             std::shared_ptr<VoxImage<Voxel>>>(m, "VoxImage")
      .def(py::init<>())
      .def(py::init<const Vector3i &>())
      .def("Data", &VoxImage<Voxel>::Data,
           py::return_value_policy::reference_internal)
      .def("InitVoxels", &VoxImage<Voxel>::InitVoxels)
      .def("Abs", &VoxImage<Voxel>::Abs)
      .def("clipImage", py::overload_cast<const Vector3i, const Vector3i>(
                            &VoxImage<Voxel>::clipImage, py::const_))
      .def("clipImage", py::overload_cast<const HPoint3f, const HPoint3f>(
                            &VoxImage<Voxel>::clipImage, py::const_))
      .def("clipImage", py::overload_cast<const float>(
                            &VoxImage<Voxel>::clipImage, py::const_))
      .def("maskImage",
           py::overload_cast<const HPoint3f, const HPoint3f, float>(
               &VoxImage<Voxel>::maskImage, py::const_))
      .def("maskImage",
           py::overload_cast<const float, float, float>(
               &VoxImage<Voxel>::maskImage, py::const_),
           py::arg("threshold"), py::arg("belowValue") = 0,
           py::arg("aboveValue") = 0)
      .def("fixVoxels", py::overload_cast<const float, float>(
                            &VoxImage<Voxel>::fixVoxels, py::const_))
      .def("__getitem__",
           py::overload_cast<unsigned int>(&VoxImage<Voxel>::operator[]))
      .def("__getitem__",
           py::overload_cast<const Vector3i &>(&VoxImage<Voxel>::operator[]));

  py::class_<TriangleMesh>(m, "TriangleMesh")
      .def(py::init<>())
      .def("AddPoint", &TriangleMesh::AddPoint)
      .def("AddTriangle",
           py::overload_cast<const Vector3i &>(&TriangleMesh::AddTriangle))
      .def("Points", py::overload_cast<>(&TriangleMesh::Points),
           py::return_value_policy::reference_internal)
      .def("Triangles", py::overload_cast<>(&TriangleMesh::Triangles),
           py::return_value_policy::reference_internal)
      .def("GetTriangle", &TriangleMesh::GetTriangle)
      .def("GetNormal", &TriangleMesh::GetNormal);

  py::class_<QuadMesh>(m, "QuadMesh")
      .def(py::init<>())
      .def("AddPoint", &QuadMesh::AddPoint)
      .def("AddQuad",
           py::overload_cast<const Vector4i &>(&QuadMesh::AddQuad))
      .def("Points", py::overload_cast<>(&QuadMesh::Points),
           py::return_value_policy::reference_internal)
      .def("Quads", py::overload_cast<>(&QuadMesh::Quads),
           py::return_value_policy::reference_internal)
      .def("GetQuad", &QuadMesh::GetQuad)
      .def("GetNormal", &QuadMesh::GetNormal);

  py::class_<VoxRaytracer::RayData::Element>(m, "VoxRaytracerRayDataElement")
      .def(py::init<>())
      .def_readwrite("vox_id", &VoxRaytracer::RayData::Element::vox_id)
      .def_readwrite("L", &VoxRaytracer::RayData::Element::L);

  py::class_<VoxRaytracer::RayData>(m, "VoxRaytracerRayData")
      .def(py::init<>())
      .def("AppendRay", &VoxRaytracer::RayData::AppendRay)
      .def("Data", py::overload_cast<>(&VoxRaytracer::RayData::Data),
           py::return_value_policy::reference_internal)
      .def("Count", &VoxRaytracer::RayData::Count)
      .def("TotalLength", &VoxRaytracer::RayData::TotalLength)
      .def("SetCount", &VoxRaytracer::RayData::SetCount)
      .def("SetTotalLength", &VoxRaytracer::RayData::SetTotalLength);

  py::class_<VoxRaytracer>(m, "VoxRaytracer")
      .def(py::init<StructuredGrid &>(), py::keep_alive<1, 2>())
      .def("GetImage", &VoxRaytracer::GetImage,
           py::return_value_policy::reference_internal)
      .def("TraceLine", &VoxRaytracer::TraceLine)
      .def("TraceBetweenPoints", &VoxRaytracer::TraceBetweenPoints);
}