numerical_analysis 0.1.2

fn multilinear_interpolation_weights<const N: usize>(
    point: &[f32; N],
    starts: &[f32; N],
    ends: &[f32; N],
) -> [f32; N]
{
    let mut weights = [0.0; N];
    for i in 0..N
    {
        weights[i] = (point[i] - starts[i]) / (ends[i] - starts[i]);
    }

    weights
}

fn grid_coordinates<const N: usize>(
    point: &[f32; N],
    corner: &[f32; N],
    extents: &[f32; N],
    step_size: f32,
) -> [usize; N]
{
    let mut coords = [0; N];
    for i in 0..N
    {
        coords[i] = ((point[i] - corner[i]) / (step_size)) as usize;
        coords[i] = coords[i].max(0).min((extents[i] / step_size) as usize - 1);
    }

    coords
}

pub fn sample_multilinear<const N: usize>(
    point: &[f32; N],
    corner: &[f32; N],
    extents: &[f32; N],
    step_size: f32,
    data: &[f32],
) -> f32
{
    let mut resolutions = [0; N];
    for i in 0..N
    {
        resolutions[i] = (extents[i] / step_size) as usize;
    }

    assert_eq!(resolutions.iter().copied().product::<usize>(), data.len());

    sample_multilinear_with_resolutions(point, corner, extents, &resolutions, data)
}

pub fn sample_multilinear_with_resolutions<const N: usize>(
    point: &[f32; N],
    corner: &[f32; N],
    extents: &[f32; N],
    resolutions: &[usize; N],
    data: &[f32],
) -> f32
{
    let mut max_extent = extents[0];
    for f in extents.iter().skip(1)
    {
        max_extent = max_extent.max(*f);
    }
    let largest_index = extents.iter().position(|f| *f == max_extent).unwrap();
    let step_size = extents[largest_index] / resolutions[largest_index] as f32;
    let coords = grid_coordinates(point, corner, extents, step_size);

    let mut starts = [0.0; N];
    for i in 0..N
    {
        starts[i] = coords[i] as f32 * step_size + corner[i];
    }

    let mut ends = [0.0; N];
    for i in 0..N
    {
        ends[i] = (coords[i] + 1) as f32 * step_size + corner[i];
    }

    let weights = multilinear_interpolation_weights(point, &starts, &ends);

    let corner_total = 1 << N;
    let mut result = 0.0;
    // Iterate over all 2^n corners of the hypercube.
    for corner in 0..corner_total
    {
        let mut coord = coords;
        // Evaluate the polynomial for this corner.
        let mut polynomial = 1.0;
        for bit_index in 0..N
        {
            // Depending on the corner the current factor of the polynomial is either
            // t or (1 - t). Use the bit pattern of the corner index to identify which is
            // it.
            let bit = ((corner >> (N - 1 - bit_index)) & 1) as i32;

            coord[bit_index] += bit as usize;
            coord[bit_index] = coord[bit_index].clamp(0, resolutions[bit_index] - 1);

            polynomial *= if bit == 1
            {
                weights[bit_index]
            }
            else
            {
                1.0 - weights[bit_index]
            };
        }

        // Index into the N dimensional array.
        let mut index = 0;
        for i in (0..N).rev()
        {
            index *= resolutions[i];
            index += coord[i];
        }

        // Multilinear interpolation of the data.
        result += data[index] * polynomial;
    }

    result
}

// +| Tests |+ =================================================================
#[cfg(test)]
mod tests
{
    use image::{ImageBuffer, Rgba, RgbaImage};
    use misc_iterators::IterDimensions;
    use noise::perlin;

    use super::sample_multilinear;

    #[test]
    fn test_2d_sampling()
    {
        let resolution = 64;
        let function_sampling = |x: f32, y: f32| {
            let value = perlin(&[x, y], 0);
            value
        };

        let mut img: RgbaImage = ImageBuffer::new(resolution as u32, resolution as u32);
        let mut image_data = Vec::with_capacity(resolution * resolution);
        for [i, j] in IterDimensions::new([resolution, resolution])
        {
            let value = function_sampling(
                i as f32 * 4.0 / resolution as f32,
                j as f32 * 4.0 / resolution as f32,
            );
            let pixel = (value) * 255.0;
            image_data.push(pixel);

            let pixel = pixel as u8;
            img.put_pixel(i as u32, j as u32, Rgba::<u8>([pixel, pixel, pixel, 255]));
        }

        img.save("input_test.png").unwrap();

        let out_res = 512;
        let mut img: RgbaImage = ImageBuffer::new(out_res as u32, out_res as u32);
        for [i, j] in IterDimensions::new([out_res, out_res])
        {
            let x = i as f32 / out_res as f32;
            let y = j as f32 / out_res as f32;

            let value = sample_multilinear(
                &[x, y],
                &[0.0, 0.0],
                &[1.0, 1.0],
                1.0 / resolution as f32,
                &image_data,
            );

            let pixel = value as u8;
            img.put_pixel(i as u32, j as u32, Rgba::<u8>([pixel, pixel, pixel, 255]));
        }

        img.save("output_test.png").unwrap();
    }
}