pristine/src/main.rs

// Copyright 2022 The Fuchsia Authors. All rights reserved.
// Use of this source code is governed by a BSD-style license that can be
// found in the LICENSE file.

use std::{fs::File, io::Write, path::PathBuf};

use arrayvec::ArrayVec;
use clap::{CommandFactory, ErrorKind, Parser};
use rayon::prelude::*;
use roots::Roots;

const LINES_PER_PIXEL: usize = 58;
const INCREMENT: f64 = 1.0 / (LINES_PER_PIXEL as f64);

#[derive(Parser)]
#[clap(about = "high-precision Bézier renderer")]
struct Config {
    /// Width of the image
    #[clap(short, long, default_value = "1000")]
    width: usize,
    /// Height of the image
    #[clap(short, long, default_value = "1000")]
    height: usize,
    /// .pgm output file
    #[clap(short, long, parse(from_os_str), default_value = "ref.pgm")]
    file: PathBuf,
    points: Vec<f64>,
}

fn roots_at(y: f64, points: &[f64]) -> ArrayVec<f64, 3> {
    let roots = match *points {
        [a, b] => roots::find_roots_linear(a, b - y),
        [a, b, c] => roots::find_roots_quadratic(a, b, c - y),
        [a, b, c, d] => roots::find_roots_cubic(a, b, c, d - y),
        _ => unreachable!(),
    };

    match roots {
        Roots::No(_) => ArrayVec::new(),
        Roots::One(roots) => roots.as_slice().try_into().unwrap(),
        roots::Roots::Two(roots) => roots.as_slice().try_into().unwrap(),
        roots::Roots::Three(roots) => roots.as_slice().try_into().unwrap(),
        _ => unreachable!(),
    }
}

fn lerp(t: f64, a: f64, b: f64) -> f64 {
    t.mul_add(b, (-t).mul_add(a, a))
}

fn eval(t: f64, points: &[f64]) -> f64 {
    match points {
        [a, b] => lerp(t, *a, *b),
        points => lerp(t, eval(t, &points[..points.len() - 1]), eval(t, &points[1..])),
    }
}

fn intersections_at(y: f64, points: &[f64]) -> ArrayVec<f64, 3> {
    let roots = match *points {
        [_, y0, _, y1] => roots_at(y, &[y1 - y0, y0]),
        [_, y0, _, y1, _, y2] => roots_at(y, &[y0 - 2.0 * y1 + y2, 2.0 * (y1 - y0), y0]),
        [_, y0, _, y1, _, y2, _, y3] => roots_at(
            y,
            &[y3 + 3.0 * (y1 - y2) - y0, 3.0 * (y0 - 2.0 * y1 + y2), 3.0 * (y1 - y0), y0],
        ),
        _ => unreachable!(),
    };

    roots
        .into_iter()
        .filter_map(|t| {
            (0.0..=1.0).contains(&t).then(|| match *points {
                [x0, _, x1, _] => eval(t, &[x0, x1]),
                [x0, _, x1, _, x2, _] => eval(t, &[x0, x1, x2]),
                [x0, _, x1, _, x2, _, x3, _] => eval(t, &[x0, x1, x2, x3]),
                _ => unreachable!(),
            })
        })
        .collect()
}

fn linear_to_srgb(l: f32) -> u8 {
    // From Hacker's Delight, p. 378-380. 2 ^ 23 as f32.
    const C23: u32 = 0x4B00_0000;

    let l = if l <= 0.0031308 { 12.92 * l } else { 1.055 * l.powf(1.0 / 2.4) - 0.055 };

    let max = u8::MAX as f32;

    let scaled = (l * max).clamp(0.0, max);
    let val = scaled + f32::from_bits(C23);

    val.to_bits() as u8
}

fn render_line(line: &mut [u8], y: usize, points: &[f64]) {
    let mut coverage = vec![0.0; line.len()];
    let intersections: Vec<_> = (0..LINES_PER_PIXEL)
        .into_iter()
        .map(|i| {
            let y = y as f64 + i as f64 * INCREMENT;

            let mut intersections = intersections_at(y, points);

            if intersections.len() % 2 != 0 {
                if !intersections_at(y - INCREMENT, points).is_empty()
                    && !intersections_at(y + INCREMENT, points).is_empty()
                {
                    intersections.push(line.len() as f64);
                } else {
                    intersections.clear();
                }
            }

            intersections
        })
        .flatten()
        .collect();

    for x in intersections.chunks(2) {
        for c in &mut coverage[x[0].ceil() as usize..x[1].floor() as usize] {
            *c += INCREMENT as f32;
        }

        if x[0].floor() != x[0].ceil() {
            coverage[x[0].floor() as usize] += ((1.0 - x[0].fract()) * INCREMENT) as f32;
        }

        if x[1].floor() != x[1].ceil() {
            coverage[x[1].floor() as usize] += (x[1].fract() * INCREMENT) as f32;
        }
    }

    for (x, c) in line.iter_mut().zip(coverage.into_iter()) {
        *x = linear_to_srgb(1.0 - c)
    }
}

fn main() {
    let config = Config::parse();

    if config.points.len() < 4 {
        Config::command().error(ErrorKind::TooFewValues, "[POINTS] need to be at least 2").exit();
    }

    if config.points.len() > 8 {
        Config::command().error(ErrorKind::TooFewValues, "[POINTS] need to be at most 8").exit();
    }

    if config.points.len() % 2 != 0 {
        Config::command()
            .error(ErrorKind::TooFewValues, "[POINTS] need to come in x y pairs")
            .exit();
    }

    let mut buffer = vec![255u8; config.width * config.height];

    buffer.par_chunks_mut(config.width).enumerate().for_each(|(y, line)| {
        render_line(line, y, &config.points);
    });

    let mut output = File::options().write(true).create(true).open(config.file).unwrap();
    output.write_all(format!("P5\n{} {}\n255\n", config.width, config.height).as_bytes()).unwrap();
    output.write_all(&buffer).unwrap();
}