crater/csg/

fields.rs

//! Scalar fields and isosurface definitions.
//!
//! This module provides the core abstractions for scalar fields and isosurfaces.
//! For theoretical background on scalar fields and CSG operations, see the
//! [Scalar Fields](../../book/theory/scalar-fields.md) chapter.

use burn::prelude::*;
use dyn_clone::DynClone;
use std::fmt::Debug;

use crate::{
    analysis::prelude::*,
    primitives::nvector::{NVector, length, normalize, to_tensor},
};

#[derive(Debug, Clone, Copy, PartialEq)]
/// Represents the classification of a point relative to an [`Isosurface`].
///
/// [`Classification::Inside`] always owns a negative value, and
/// [`Classification::Outside`] always owns a positive value.
/// [`Classification::On`] is a special case that occurs when the field value is within
/// [`EPSILON`] of the surface.
pub enum Classification {
    Inside(f32),
    Outside(f32),
    On,
}

impl Classification {
    pub fn is_inside(&self) -> bool {
        matches!(self, Classification::Inside(_))
    }

    pub fn is_outside(&self) -> bool {
        matches!(self, Classification::Outside(_))
    }

    pub fn is_on(&self) -> bool {
        matches!(self, Classification::On)
    }
}

impl From<f32> for Classification {
    fn from(value: f32) -> Self {
        if value < -EPSILON {
            Classification::Inside(value)
        } else if value > EPSILON {
            Classification::Outside(value)
        } else {
            Classification::On
        }
    }
}

pub trait Classify<B: Backend, const D: usize> {
    /// Mask for points inside the surface.
    fn is_inside_mask(&self) -> Tensor<B, D, Bool>;

    /// Mask for points outside the surface.
    fn is_outside_mask(&self) -> Tensor<B, D, Bool>;

    /// Mask for points on the surface.
    fn is_on_mask(&self) -> Tensor<B, D, Bool>;

    /// Get the classification of a point at a given index.
    fn classification_of_index(&self, index: usize) -> Classification;
}

impl<B: Backend, const D: usize> Classify<B, D> for Tensor<B, D, Float> {
    fn is_inside_mask(&self) -> Tensor<B, D, Bool> {
        self.clone().lower_elem(-EPSILON)
    }

    fn is_outside_mask(&self) -> Tensor<B, D, Bool> {
        self.clone().greater_elem(EPSILON)
    }

    fn is_on_mask(&self) -> Tensor<B, D, Bool> {
        self.clone().abs().lower_elem(EPSILON)
    }

    fn classification_of_index(&self, index: usize) -> Classification {
        match self.clone().into_data().as_slice().unwrap()[index] {
            x if x < -EPSILON => Classification::Inside(x),
            x if x > EPSILON => Classification::Outside(x),
            _ => Classification::On,
        }
    }
}

/// A collection of points in N-dimensional space
#[allow(type_alias_bounds)]
pub type Origins<B: Backend, const N: usize> = Tensor<B, 2, Float>;
/// A collection of directions in N-dimensional space
#[allow(type_alias_bounds)]
pub type Directions<B: Backend, const N: usize> = Tensor<B, 2, Float>;
/// A collection of scalars
#[allow(type_alias_bounds)]
pub type Scalars<B: Backend> = Tensor<B, 1, Float>;

/// A [`ScalarField`] is a function that maps a collection of points to a collection of scalars.
pub trait ScalarField<const N: usize, B: Backend>: DynClone + Send + Sync {
    /// Evaluates the scalar field at a collection of points.
    fn evaluate(&self, origins: Origins<B, N>) -> Scalars<B>;

    /// Get the device on which the scalar field is allocated.
    fn device(&self) -> &B::Device;
}

dyn_clone::clone_trait_object!(<const N: usize, B: Backend> ScalarField<N, B>);

/// Absolute value for a floating point value of a scalar field to be considered zero.
// pub const EPSILON: f32 = 5e-3;
pub const EPSILON: f32 = 1e-5;

/// Represents a surface in N-dimensional space. The surface is defined by a mathematical
/// expression that takes N variables as input.
#[derive(Clone)]
pub struct Isosurface<const N: usize, B: Backend> {
    /// The scalar field
    pub field: Box<dyn ScalarField<N, B>>,
    /// The ray field (optional)
    pub ray_field: Option<Box<dyn RayField<N, B>>>,
    /// The constant value of the surface.
    pub constant: Tensor<B, 1, Float>,
}

impl<const N: usize, B: Backend> PartialEq for Isosurface<N, B> {
    fn eq(&self, _other: &Self) -> bool {
        // All surfaces are unequal
        false
    }
}

impl<const N: usize, B: Backend> Debug for Isosurface<N, B> {
    fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
        write!(f, "Isosurface<{N}>")
    }
}

impl<const N: usize, B: Backend> ScalarField<N, B> for Isosurface<N, B> {
    fn evaluate(&self, points: Origins<B, N>) -> Scalars<B> {
        self.field.evaluate(points) - self.constant.clone()
    }

    fn device(&self) -> &B::Device {
        self.field.device()
    }
}

impl<const N: usize, B: Backend> Isosurface<N, B> {
    /// Classifies a single point relative to the [`Isosurface`].
    pub fn classify_point(&self, point: &NVector<N>) -> Classification {
        let points = Tensor::from_data([*point], self.device());
        self.evaluate(points).classification_of_index(0)
    }
}

/// Conversion of a type into a [`Isosurface`] object.
pub trait IntoIsosurface<const N: usize, B: Backend> {
    /// Converts the object into a [`Isosurface<N, B>`].
    fn into_isosurface(self, constant: f32) -> Isosurface<N, B>;
}

/// Represents a 2D geometric shape that can be converted into a [`Isosurface`].
#[derive(Debug, Clone)]
pub enum Field2D<B: Backend> {
    /// A circle defined by its radius
    Circle {
        r: Tensor<B, 2, Float>,
        device: B::Device,
    },
    /// An ellipse defined by its semi-major and semi-minor axes and center.
    Ellipse {
        a: Tensor<B, 2, Float>,
        b: Tensor<B, 2, Float>,
        device: B::Device,
    },
    /// Line in 2D space
    Line {
        normal: Tensor<B, 2, Float>,
        device: B::Device,
    },
    /// A cylinder in 2D space (infinite circle)
    Cylinder {
        r: Tensor<B, 2, Float>,
        device: B::Device,
    },
}

impl<B: Backend> Field2D<B> {
    /// Allocate a circle [`ScalarField`] in the given device's memory.
    pub fn circle(r: f32, device: B::Device) -> Self {
        Self::Circle {
            r: to_tensor([r], &device.clone()),
            device: device.clone(),
        }
    }

    /// Allocate an ellipse [`ScalarField`] in the given device's memory.
    pub fn ellipse(a: f32, b: f32, device: B::Device) -> Self {
        Self::Ellipse {
            a: to_tensor([a], &device.clone()),
            b: to_tensor([b], &device.clone()),
            device: device.clone(),
        }
    }

    /// Allocate a line [`ScalarField`] in the given device's memory.
    pub fn line(normal: NVector<2>, device: B::Device) -> Self {
        Self::Line {
            normal: to_tensor(normal, &device.clone()),
            device: device.clone(),
        }
    }

    /// Allocate a cylinder [`ScalarField`] in the given device's memory.
    /// This is equivalent to a hypercylinder in 2D space.
    pub fn cylinder(r: f32, device: B::Device) -> Self {
        Self::Cylinder {
            r: to_tensor([r], &device.clone()),
            device: device.clone(),
        }
    }
}

impl<B: Backend> ScalarField<2, B> for Field2D<B> {
    fn evaluate(&self, points: Origins<B, 2>) -> Scalars<B> {
        assert!(
            points.shape().dims::<2>()[1] == 2,
            "Points must be of shape (num_points, 2)"
        );

        match self {
            Field2D::Circle { r, .. } => FieldND::<2, B>::Hypersphere {
                r: r.clone(),
                device: self.device().clone(),
            }
            .evaluate(points),
            Field2D::Ellipse { a, b, .. } => (points.clone().slice([0..1, 0..1])
                / a.clone().unsqueeze().powi_scalar(2)
                + points.clone().slice([0..1, 1..2]) / b.clone().unsqueeze().powi_scalar(2)
                - 1.0)
                .squeeze(1),

            Field2D::Line { normal, .. } => FieldND::<2, B>::Hyperplane {
                normal: normal.clone(),
                device: self.device().clone(),
            }
            .evaluate(points),
            Field2D::Cylinder { r, .. } => FieldND::<2, B>::Hypercylinder {
                r: r.clone(),
                device: self.device().clone(),
            }
            .evaluate(points),
        }
    }

    fn device(&self) -> &B::Device {
        match self {
            Field2D::Circle { device, .. } => device,
            Field2D::Ellipse { device, .. } => device,
            Field2D::Line { device, .. } => device,
            Field2D::Cylinder { device, .. } => device,
        }
    }
}

impl<B: Backend> IntoIsosurface<2, B> for Field2D<B> {
    fn into_isosurface(self, constant: f32) -> Isosurface<2, B> {
        Isosurface {
            field: Box::new(self.clone()),
            ray_field: Some(Box::new(self.clone())),
            constant: Tensor::from_data([constant], ScalarField::device(&self)),
        }
    }
}

#[derive(Debug, Clone)]
pub enum Field3D<B: Backend> {
    /// A cone defined by its direction and opening angle
    Cone {
        axis: Tensor<B, 2, Float>,
        theta: Tensor<B, 2, Float>,
        device: B::Device,
    },
    /// A cylinder defined by its radius
    Cylinder {
        r: Tensor<B, 2, Float>,
        device: B::Device,
    },
    /// A sphere defined by its radius.
    Sphere {
        r: Tensor<B, 2, Float>,
        device: B::Device,
    },
    /// A torus defined by its two radii.
    Torus {
        r1: Tensor<B, 2, Float>,
        r2: Tensor<B, 2, Float>,
        device: B::Device,
    },
    /// Planes in 3D space
    Plane {
        normal: Tensor<B, 2, Float>,
        device: B::Device,
    },
}

impl<B: Backend> Field3D<B> {
    /// Creates a 3D sphere scalar field with specified radius.
    ///
    /// Generates a sphere centered at the origin with the given radius. The scalar field
    /// represents the signed distance from any point to the sphere surface, where:
    /// - Negative values indicate points inside the sphere
    /// - Zero values indicate points exactly on the sphere surface  
    /// - Positive values indicate points outside the sphere
    ///
    /// # Mathematical Definition
    ///
    /// For a point **x** = (x, y, z) and sphere radius r, the field function is:
    ///
    /// f(**x**) = ||**x**|| - r = √(x² + y² + z²) - r
    ///
    /// This is the signed distance function (SDF) for a sphere, providing exact
    /// geometric distances that enable precise CSG operations.
    ///
    /// # Parameters
    ///
    /// * `r` - Sphere radius in world units. Must be positive.
    /// * `device` - Backend device for tensor allocation and computation
    ///
    /// # Returns
    ///
    /// A [`Field3D`] representing the sphere scalar field, ready for use in
    /// CSG operations, isosurface extraction, or ray casting.
    ///
    /// # Panics
    ///
    /// This function does not panic, but negative radii will produce inverted spheres
    /// where the interior and exterior are swapped.
    ///
    /// # Examples
    ///
    /// ```rust
    /// use crater::csg::prelude::*;
    /// use crater::primitives::prelude::*;
    /// use burn::prelude::*;
    /// use burn::backend::ndarray::NdArrayDevice;
    ///
    /// // Create a unit sphere
    /// let sphere = Field3D::<burn::backend::ndarray::NdArray>::sphere(1.0, NdArrayDevice::Cpu);
    ///
    /// // Convert to isosurface for CSG operations
    /// let isosurface = sphere.into_isosurface(0.0);
    /// let region = Region::HalfSpace(isosurface, Side::Negative);
    ///
    /// // Test point classification
    /// let algebra = Algebra::default();
    /// let center = Tensor::from_data([[0.0, 0.0, 0.0]], &NdArrayDevice::Cpu);
    /// let surface = Tensor::from_data([[1.0, 0.0, 0.0]], &NdArrayDevice::Cpu);
    /// let exterior = Tensor::from_data([[2.0, 0.0, 0.0]], &NdArrayDevice::Cpu);
    ///
    /// let center_val = region.evaluate(center, &algebra);    // ≈ -1.0 (inside)
    /// let surface_val = region.evaluate(surface, &algebra);  // ≈  0.0 (on surface)  
    /// let exterior_val = region.evaluate(exterior, &algebra); // ≈  1.0 (outside)
    /// ```
    ///
    /// ```rust
    /// use crater::csg::prelude::*;
    /// use crater::primitives::prelude::*;
    /// use burn::prelude::*;
    /// use burn::backend::ndarray::NdArrayDevice;
    ///
    /// // Create a large sphere for mesh generation
    /// let large_sphere = Field3D::<burn::backend::ndarray::NdArray>::sphere(5.0, NdArrayDevice::Cpu).into_isosurface(0.0);
    /// let region = Region::HalfSpace(large_sphere, Side::Negative);
    ///
    /// // Generate high-resolution mesh
    /// let params = MarchingCubesParams {
    ///     region,
    ///     bounds: BoundingBox::new([-6.0, -6.0, -6.0], [6.0, 6.0, 6.0]),
    ///     resolution: (64, 64, 64),
    ///     algebra: Algebra::default(),
    /// };
    ///
    /// let mesh = marching_cubes(&params, &NdArrayDevice::Cpu);
    /// println!("Generated sphere mesh with {} triangles", mesh.triangles.len());
    /// ```
    pub fn sphere(r: f32, device: B::Device) -> Self {
        Self::Sphere {
            r: to_tensor([r], &device.clone()),
            device,
        }
    }

    /// Allocate a cylinder [`ScalarField`] in the given device's memory.
    pub fn cylinder(r: f32, device: B::Device) -> Self {
        Self::Cylinder {
            r: to_tensor([r], &device.clone()),
            device,
        }
    }

    /// Allocate a cone [`ScalarField`] in the given device's memory.
    pub fn cone(direction: NVector<3>, theta: f32, device: B::Device) -> Self {
        assert!(theta >= 0.0, "Theta must be non-negative");
        assert!(length(&direction) > 0.0, "direction must be non-zero");
        Self::Cone {
            axis: to_tensor(normalize(&direction), &device.clone()),
            theta: to_tensor([theta], &device.clone()),
            device: device.clone(),
        }
    }

    /// Allocate a torus [`ScalarField`] in the given device's memory.
    pub fn torus(r1: f32, r2: f32, device: B::Device) -> Self {
        Self::Torus {
            // Radii must be a column vectors
            r1: to_tensor([r1], &device.clone()),
            r2: to_tensor([r2], &device.clone()),
            device: device.clone(),
        }
    }

    /// Allocate a plane [`ScalarField`] in the given device's memory.      
    pub fn plane(normal: NVector<3>, device: B::Device) -> Self {
        Self::Plane {
            // Normal vector must be a column vector
            normal: to_tensor(normal, &device.clone()),
            device: device.clone(),
        }
    }

    /// Allocate a hypercylinder [`ScalarField`] in the given device's memory.
    /// This is equivalent to a hypercylinder in 3D space.
    pub fn hypercylinder(r: f32, device: B::Device) -> Self {
        Self::Cylinder {
            r: to_tensor([r], &device.clone()),
            device: device.clone(),
        }
    }

    /// Allocate a hypertorus [`ScalarField`] in the given device's memory.
    /// This is equivalent to a hypertorus in 3D space.
    pub fn hypertorus(r1: f32, r2: f32, device: B::Device) -> Self {
        Self::Torus {
            r1: to_tensor([r1], &device.clone()),
            r2: to_tensor([r2], &device.clone()),
            device: device.clone(),
        }
    }
}

impl<B: Backend> ScalarField<3, B> for Field3D<B> {
    fn evaluate(&self, points: Origins<B, 3>) -> Scalars<B> {
        assert!(
            points.shape().dims::<2>()[1] == 3,
            "Points must be of shape (num_points, 3)"
        );

        match self {
            Field3D::Sphere { r, .. } => FieldND::<3, B>::Hypersphere {
                r: r.clone(),
                device: self.device().clone(),
            }
            .evaluate(points),
            Field3D::Cone { axis, theta, .. } => FieldND::<3, B>::Hypercone {
                axis: axis.clone(),
                theta: theta.clone(),
                device: self.device().clone(),
            }
            .evaluate(points),
            Field3D::Cylinder { r, .. } => FieldND::<3, B>::Hypercylinder {
                r: r.clone(),
                device: self.device().clone(),
            }
            .evaluate(points),
            Field3D::Torus { r1, r2, .. } => FieldND::<3, B>::Hypertorus {
                r1: r1.clone(),
                r2: r2.clone(),
                device: self.device().clone(),
            }
            .evaluate(points),
            Field3D::Plane { normal, .. } => FieldND::<3, B>::Hyperplane {
                normal: normal.clone(),
                device: self.device().clone(),
            }
            .evaluate(points),
        }
    }

    fn device(&self) -> &B::Device {
        match self {
            Field3D::Sphere { device, .. } => device,
            Field3D::Cone { device, .. } => device,
            Field3D::Cylinder { device, .. } => device,
            Field3D::Torus { device, .. } => device,
            Field3D::Plane { device, .. } => device,
        }
    }
}

impl<B: Backend> IntoIsosurface<3, B> for Field3D<B> {
    fn into_isosurface(self, constant: f32) -> Isosurface<3, B> {
        Isosurface {
            field: Box::new(self.clone()),
            ray_field: Some(Box::new(self.clone())),
            constant: Tensor::from_data([constant], ScalarField::device(&self)),
        }
    }
}

/// Represents a surface in N-dimensional space. The surface is defined by a mathematical
/// expression that takes N variables as input.
#[derive(Clone)]
pub enum FieldND<const N: usize, B: Backend> {
    /// A hyperplane defined by its normal vector and distance from the origin.
    Hyperplane {
        normal: Tensor<B, 2, Float>,
        device: B::Device,
    },

    /// A hypersphere defined by its radius.
    Hypersphere {
        r: Tensor<B, 2, Float>,
        device: B::Device,
    },
    /// A hypercone defined by its axis and opening angle.
    Hypercone {
        axis: Tensor<B, 2, Float>,
        theta: Tensor<B, 2, Float>,
        device: B::Device,
    },
    /// A hypercylinder defined by its radius, extending along the last dimension.
    Hypercylinder {
        r: Tensor<B, 2, Float>,
        device: B::Device,
    },

    /// A hypertorus defined by its two radii.
    Hypertorus {
        r1: Tensor<B, 2, Float>,
        r2: Tensor<B, 2, Float>,
        device: B::Device,
    },
}

impl<const N: usize, B: Backend> FieldND<N, B> {
    /// Allocate a hyperplane [`ScalarField`] in the given device's memory.
    pub fn hyperplane(normal: NVector<N>, device: B::Device) -> Self {
        Self::Hyperplane {
            normal: to_tensor(normalize(&normal), &device.clone()),
            device: device.clone(),
        }
    }

    /// Allocate a hypersphere [`ScalarField`] in the given device's memory.
    pub fn hypersphere(r: f32, device: B::Device) -> Self {
        Self::Hypersphere {
            r: to_tensor([r], &device.clone()),
            device: device.clone(),
        }
    }

    /// Allocate a hypercylinder [`ScalarField`] in the given device's memory.
    /// The cylinder extends infinitely along the last dimension.
    pub fn hypercylinder(r: f32, device: B::Device) -> Self {
        Self::Hypercylinder {
            r: to_tensor([r], &device.clone()),
            device: device.clone(),
        }
    }

    /// Allocate a hypertorus [`ScalarField`] in the given device's memory.
    /// r1 is the major radius, r2 is the minor radius.
    /// The torus lies in the first N-1 dimensions.
    pub fn hypertorus(r1: f32, r2: f32, device: B::Device) -> Self {
        Self::Hypertorus {
            r1: to_tensor([r1], &device.clone()),
            r2: to_tensor([r2], &device.clone()),
            device: device.clone(),
        }
    }
}

impl<const N: usize, B: Backend> ScalarField<N, B> for FieldND<N, B> {
    fn evaluate(&self, points: Origins<B, N>) -> Scalars<B> {
        assert!(
            points.shape().dims::<2>()[1] == N,
            "Points must be of shape (num_points, {N})"
        );

        match self {
            FieldND::Hyperplane { normal, .. } => points.matmul(normal.clone()).squeeze(1),
            FieldND::Hypersphere { r, .. } => (points.clone().powf_scalar(2.0).sum_dim(1)
                - r.clone().powf_scalar(2.0).unsqueeze())
            .squeeze(1),
            FieldND::Hypercone { axis, theta, .. } => {
                let axis_dot_x = points.clone().matmul(axis.clone()).squeeze(1);
                let cos_theta_sq = theta.clone().cos().powf_scalar(2.0).squeeze(1);
                let x_squared = points.clone().mul(points.clone()).sum_dim(1).squeeze(1);

                // Double cone: cos²(θ) * |x|² - (axis·x)² = 0
                // This defines a cone that extends in both directions from the origin
                cos_theta_sq.clone() * x_squared - axis_dot_x.clone().powf_scalar(2.0)
            }
            FieldND::Hypercylinder { r, .. } => {
                // Hypercylinder: sum of squares of first N-1 dimensions minus r²
                // This creates a cylinder that extends infinitely along the last dimension
                if N <= 1 {
                    panic!("Hypercylinder requires at least 2 dimensions");
                }
                let first_n_minus_1 = points.clone().slice([None, Some((0i64, (N - 1) as i64))]);
                (first_n_minus_1.powf_scalar(2.0).sum_dim(1)
                    - r.clone().powf_scalar(2.0).unsqueeze())
                .squeeze(1)
            }
            FieldND::Hypertorus { r1, r2, .. } => {
                // Hypertorus: generalization of the 3D torus to N dimensions
                // The torus lies in the first N-1 dimensions
                if N <= 2 {
                    panic!("HyperTorus requires at least 3 dimensions");
                }

                // Extract the last coordinate (N-1 dimension)
                let device = self.device();
                let last_coord = points
                    .clone()
                    .select(1, Tensor::<B, 1, Int>::from_data([N - 1], device));

                // Distance in the first N-1 dimensions
                let first_n_minus_1 = points.clone().slice([None, Some((0i64, (N - 1) as i64))]);
                let radial_distance_sq = first_n_minus_1.powf_scalar(2.0).sum_dim(1);
                let radial_distance = radial_distance_sq.clone().sqrt();

                // Hypertorus equation: (sqrt(x₁² + x₂² + ... + x_{N-1}²) - r1)² + x_N² - r2²
                let major_term = (radial_distance - r1.clone().unsqueeze()).powf_scalar(2.0);
                let minor_term = last_coord.powf_scalar(2.0);
                let r2_sq = r2.clone().powf_scalar(2.0).unsqueeze();

                (major_term + minor_term - r2_sq).squeeze(1)
            }
        }
    }

    fn device(&self) -> &B::Device {
        match self {
            FieldND::Hyperplane { device, .. } => device,
            FieldND::Hypersphere { device, .. } => device,
            FieldND::Hypercone { device, .. } => device,
            FieldND::Hypercylinder { device, .. } => device,
            FieldND::Hypertorus { device, .. } => device,
        }
    }
}

impl<const N: usize, B: Backend> IntoIsosurface<N, B> for FieldND<N, B> {
    fn into_isosurface(self, constant: f32) -> Isosurface<N, B> {
        Isosurface {
            field: Box::new(self.clone()),
            ray_field: Some(Box::new(self.clone())),
            constant: Tensor::from_data([constant], ScalarField::device(&self)),
        }
    }
}

#[cfg(test)]
mod tests {

    use crate::csg::prelude::*;
    use crate::primitives::nvector::to_tensor;
    use crate::test_utils::{assert_tensor_almost_eq, assert_tensor_eq};
    use backend_macro::with_backend;
    use burn::prelude::*;
    use rstest::rstest;

    #[with_backend]
    #[rstest]
    // N-D
    #[case(
        "hyperplane",
        [[0.0, 0.0, 0.0], [1.0, 0.0, 0.0], [0.0, 1.0, 0.0]],
        [0.0, 1.0, 0.0]
    )]
    #[case(
        "hypersphere",  
        [[0.0, 0.0, 0.0], [1.0, 0.0, 0.0], [0.0, 1.0, 0.0]],
        [-1.0, 0.0, 0.0]
    )]
    #[case(
        "hypercone",
        [[0.0, 0.0, 0.0], [0.0, 0.0, 1.0], [0.5, 0.0, 0.866025]],
        [0.0, -0.25, 0.0]
    )]
    #[case(
        "hypercylinder",
        [[0.0, 0.0, 0.0], [1.0, 0.0, 0.0], [0.0, 1.0, 5.0]],
        [-1.0, 0.0, 0.0]
    )]
    fn test_scalar_field_nd<const N: usize>(
        #[case] field: &str,
        #[case] points: [[f32; N]; 3],
        #[case] expected: [f32; 3],
    ) {
        let field = match field {
            "hyperplane" => FieldND::<3, Backend>::hyperplane([1.0, 0.0, 0.0], device()),
            "hypersphere" => FieldND::<3, Backend>::hypersphere(1.0, device()),
            "hypercone" => {
                let axis = to_tensor([0.0, 0.0, 1.0], &device());
                let theta = to_tensor([std::f32::consts::FRAC_PI_6], &device());
                FieldND::<3, Backend>::Hypercone {
                    axis,
                    theta,
                    device: device(),
                }
            }
            "hypercylinder" => FieldND::<3, Backend>::hypercylinder(1.0, device()),
            _ => panic!("Invalid field"),
        };
        let points = Tensor::from_data(points, ScalarField::device(&field));
        let expected = Tensor::from_data(expected, ScalarField::device(&field));
        let values = field.evaluate(points);
        assert_tensor_almost_eq(values, expected, Some(EPSILON));
    }

    #[with_backend]
    #[rstest]
    #[case(
        "sphere",
        [[0.0, 0.0, 0.0], [1.0, 0.0, 0.0], [0.0, 1.0, 0.0]],
        [-1.0, 0.0, 0.0]
    )]
    #[case(
        "cone",
        [[0.0, 0.0, 0.0], [0.0, 0.0, 1.0], [0.5, 0.0, 0.866025]],
        [0.0, -0.25, 0.0]
    )]
    fn test_scalar_field_3d(
        #[case] field: &str,
        #[case] points: [[f32; 3]; 3],
        #[case] expected: [f32; 3],
    ) {
        let field = match field {
            "sphere" => Field3D::<Backend>::sphere(1.0, device()),
            "cone" => {
                Field3D::<Backend>::cone([0.0, 0.0, 1.0], std::f32::consts::FRAC_PI_6, device())
            }
            _ => panic!("Invalid field"),
        };
        let points = Tensor::from_data(points, ScalarField::device(&field));
        let expected = Tensor::from_data(expected, ScalarField::device(&field));
        let values = field.evaluate(points);
        println!("values: {}", values);
        println!("expected: {}", expected);
        assert_tensor_almost_eq(values, expected, Some(EPSILON));
    }

    #[with_backend]
    #[rstest]
    #[case(
        "line",
        [[0.0, 0.0], [1.0, 0.0], [0.0, 1.0]],
        [0.0, 1.0, 0.0]
    )]
    #[case(
        "circle",
        [[0.0, 0.0], [1.0, 0.0], [0.0, 1.0]],
        [-1.0, 0.0, 0.0]
    )]
    fn test_scalar_field_2d(
        #[case] field: &str,
        #[case] points: [[f32; 2]; 3],
        #[case] expected: [f32; 3],
    ) {
        let field = match field {
            "line" => Field2D::<Backend>::line([1.0, 0.0], device()),
            "circle" => Field2D::<Backend>::circle(1.0, device()),
            _ => panic!("Invalid field"),
        };
        let points = Tensor::from_data(points, &device());
        let expected = Tensor::from_data(expected, &device());
        let values = field.evaluate(points);
        assert_tensor_eq(values, expected);
    }

    #[with_backend]
    #[test]
    fn test_hypertorus_4d() {
        // Test 4D hypertorus with major radius 2.0, minor radius 0.5
        let hypertorus = FieldND::<4, Backend>::hypertorus(2.0, 0.5, device());

        let points = Tensor::from_data(
            [
                [2.0, 0.0, 0.0, 0.0], // On major circle in first 3 dims
                [2.0, 0.0, 0.0, 0.5], // On torus surface
                [0.0, 0.0, 0.0, 0.0], // At origin
                [1.0, 0.0, 0.0, 0.0], // Inside major radius
            ],
            &device(),
        );

        let values = hypertorus.evaluate(points);

        // Expected values based on hypertorus equation
        let expected = Tensor::from_data(
            [
                -0.25, // (2-2)² + 0² - 0.25 = -0.25
                0.0,   // (2-2)² + 0.25 - 0.25 = 0.0 (on surface)
                3.75,  // (0-2)² + 0² - 0.25 = 4 - 0.25 = 3.75
                0.75,  // (1-2)² + 0² - 0.25 = 1 - 0.25 = 0.75
            ],
            &device(),
        );

        assert_tensor_almost_eq(values, expected, Some(1e-5));
    }
}