leibniz 0.2.0

//! Grayscale rasterization.

mod sampling;

use ::burn::tensor::{Tensor, backend::Backend};

use crate::base::Config;

use super::{
    filter,
    geometry::{Contour, contains},
};

const EPSILON: f32 = 1e-6;

/// Grayscale raster with shape `[height, width]`.
pub type Raster<B> = Tensor<B, 2>;

/// Render contours.
///
/// # Panics
///
/// Panics if there are no contours, if the target dimensions are zero, or if
/// any contour is empty or does not have segment shape `[segments, 2, 2]`.
pub fn render<B: Backend>(
    contours: &[Contour<B>],
    height: usize,
    width: usize,
    config: Config,
) -> Raster<B> {
    let x_sample_count = config.x_sample_count();
    let y_sample_count = config.y_sample_count();
    let radius = config.radius();

    assert!(
        !contours.is_empty() && height > 0 && width > 0,
        "contours and target dimensions must be non-empty"
    );

    for contour in contours {
        let [segment_count, point_count, coordinate_count] = contour.dims();

        assert!(
            segment_count > 0 && point_count == 2 && coordinate_count == 2,
            "contours must have shape [segments, 2, 2]"
        );
    }

    // Evaluate coverage on the fine sub-sample grid, then reconstruct each pixel
    // as the box-filtered average, mirroring DiffVG's render_kernel: the pixel
    // value is the filter-weighted sub-sample sum divided by the same
    // weight_image used by the boundary gather. The box filter is separable, so
    // the reconstruction is a pair of per-axis selection matmuls.
    let device = contours[0].device();
    let (x, y) =
        sampling::fine_positions::<B>(height, width, x_sample_count, y_sample_count, &device);
    let coverage = contains(contours.iter().cloned(), x, y)
        .float()
        .reshape([height * y_sample_count, width * x_sample_count]);
    let (columns, column_counts) = filter::axis::<B>(width, x_sample_count, radius, &device);
    let (rows, row_counts) = filter::axis::<B>(height, y_sample_count, radius, &device);
    let weight = row_counts.unsqueeze_dim::<2>(1) * column_counts.unsqueeze_dim::<2>(0);
    let coverage = rows.matmul(coverage).matmul(columns.transpose());

    coverage / weight.clamp_min(EPSILON)
}

#[cfg(test)]
mod tests {
    use ::burn::tensor::{Tensor, TensorData};

    use super::render;
    use crate::{
        base::{Config, geometry::Command},
        burn::{geometry::Contour, tests::Backend},
    };

    #[test]
    #[should_panic]
    fn rejects_empty_contours() {
        let _ = render::<Backend>(&[], 1, 1, Config::default());
    }

    #[test]
    #[should_panic]
    fn rejects_zero_config_x_sample_count() {
        let _ = render::<Backend>(
            &[square([0.0, 0.0], [1.0, 1.0])],
            1,
            1,
            Config::new(0, 2, 0.5),
        );
    }

    #[test]
    #[should_panic]
    fn rejects_zero_config_y_sample_count() {
        let _ = render::<Backend>(
            &[square([0.0, 0.0], [1.0, 1.0])],
            1,
            1,
            Config::new(2, 0, 0.5),
        );
    }

    #[test]
    #[should_panic]
    fn rejects_zero_height() {
        let _ = render::<Backend>(&[square([0.0, 0.0], [1.0, 1.0])], 0, 1, Config::default());
    }

    #[test]
    #[should_panic]
    fn rejects_zero_width() {
        let _ = render::<Backend>(&[square([0.0, 0.0], [1.0, 1.0])], 1, 0, Config::default());
    }

    #[test]
    fn renders_full_coverage() {
        let raster = render::<Backend>(&[square([0.0, 0.0], [2.0, 2.0])], 2, 2, Config::default());

        assert_raster(raster, [1.0, 1.0, 1.0, 1.0]);
    }

    #[test]
    fn renders_partial_coverage() {
        let raster = render::<Backend>(&[square([0.0, 0.0], [0.5, 1.0])], 1, 1, Config::default());

        assert_raster(raster, [0.5]);
    }

    #[test]
    fn renders_ring_with_nonzero_hole() {
        let raster = render::<Backend>(
            &[
                square([0.0, 0.0], [3.0, 3.0]),
                clockwise_square([1.0, 1.0], [2.0, 2.0]),
            ],
            3,
            3,
            Config::default(),
        );

        assert_raster(raster, [1.0, 1.0, 1.0, 1.0, 0.0, 1.0, 1.0, 1.0, 1.0]);
    }

    #[test]
    fn renders_with_antialiasing() {
        let contours = [square([0.0, 0.0], [0.4, 1.0])];
        let coarse = render::<Backend>(&contours, 1, 1, Config::new(1, 1, 0.5));
        let finer = render::<Backend>(&contours, 1, 1, Config::new(2, 2, 0.5));

        assert_raster(coarse, [0.0]);
        assert_raster(finer, [0.5]);
    }

    fn assert_close(actual: f32, expected: f32) {
        assert!((actual - expected).abs() < 1e-6);
    }

    fn assert_raster<const N: usize>(raster: Tensor<Backend, 2>, expected: [f32; N]) {
        let actual = raster.into_data().to_vec::<f32>().unwrap();

        assert_eq!(actual.len(), expected.len());
        for (actual, expected) in actual.into_iter().zip(expected) {
            assert_close(actual, expected);
        }
    }

    fn clockwise_square(min: [f32; 2], max: [f32; 2]) -> Contour<Backend> {
        Contour::new(
            vec![Command::Linear; 4],
            Tensor::<Backend, 3>::from_data(
                TensorData::from(square_segments(
                    [min[0], min[1]],
                    [min[0], max[1]],
                    [max[0], max[1]],
                    [max[0], min[1]],
                )),
                &Default::default(),
            ),
        )
    }

    fn interpolate(start: [f32; 2], end: [f32; 2], t: f32) -> [f32; 2] {
        [
            start[0] + (end[0] - start[0]) * t,
            start[1] + (end[1] - start[1]) * t,
        ]
    }

    fn square(min: [f32; 2], max: [f32; 2]) -> Contour<Backend> {
        Contour::new(
            vec![Command::Linear; 4],
            Tensor::<Backend, 3>::from_data(
                TensorData::from(square_segments(
                    [min[0], min[1]],
                    [max[0], min[1]],
                    [max[0], max[1]],
                    [min[0], max[1]],
                )),
                &Default::default(),
            ),
        )
    }

    fn square_segments(a: [f32; 2], b: [f32; 2], c: [f32; 2], d: [f32; 2]) -> [[[f32; 2]; 2]; 4] {
        [
            [a, interpolate(a, b, 0.5)],
            [b, interpolate(b, c, 0.5)],
            [c, interpolate(c, d, 0.5)],
            [d, interpolate(d, a, 0.5)],
        ]
    }
}