polydf 0.1.1

// Copyright (c) 2025 William Bradley
use std::error::Error;

use image::{ImageBuffer, Luma};
use indicatif::{ProgressBar, ProgressStyle};
use nalgebra::{Vector2, vector};
use polydf::{NearestOptions, is_definitely_far, nearest_t_with_derivative};

// Curve: C(t) = t [cos(theta), sin(theta)], theta(t) = sin(0.04 * sqrt(t) + 5)
fn curve(t: f32) -> Vector2<f32> {
    let s = 0.04 * t.sqrt() + 5.0;
    let theta = s.sin();
    let (st, ct) = theta.sin_cos();
    vector![t * ct, t * st]
}

// Derivative: C'(t) = [cos(theta), sin(theta)] + (0.02*sqrt(t)*cos(s)) * [-sin(theta), cos(theta)]
fn curve_d(t: f32) -> Vector2<f32> {
    let sqrt_t = t.sqrt();
    let s = 0.04 * sqrt_t + 5.0;
    let theta = s.sin();
    let (st, ct) = theta.sin_cos();
    let k = 0.02 * sqrt_t * s.cos();
    vector![ct - k * st, st + k * ct]
}

fn main() -> Result<(), Box<dyn Error>> {
    // Image parameters (overridable via env for quick tests)
    let width: u32 = std::env::var("RENDER_WIDTH")
        .ok()
        .and_then(|s| s.parse().ok())
        .unwrap_or(1024);
    let height: u32 = std::env::var("RENDER_HEIGHT")
        .ok()
        .and_then(|s| s.parse().ok())
        .unwrap_or(1024);

    // World domain
    let x_min = 0.0f32;
    let x_max = 500_000.0f32; // 5e5
    let y_min = -250_000.0f32; // -2.5e5
    let y_max = 250_000.0f32; // 2.5e5

    let threshold: f32 = std::env::var("RENDER_THRESHOLD")
        .ok()
        .and_then(|s| s.parse().ok())
        .unwrap_or(50_000.0);

    // Parameter search domain
    let t_range = 0.0f32..=500_000.0f32;

    // Conservative upper bound on |C'(t)|; derivative is still used for accuracy/speed
    let speed_upper_bound = Some(16.0f32);

    // Use a modest number of samples; increase for higher accuracy
    let samples: usize = std::env::var("RENDER_SAMPLES")
        .ok()
        .and_then(|s| s.parse().ok())
        .unwrap_or(16);
    let opts = NearestOptions {
        samples,
        ..Default::default()
    };

    let mut img: ImageBuffer<Luma<u8>, Vec<u8>> = ImageBuffer::new(width, height);

    // Precompute scales to map pixels to world coords
    let sx = if width > 1 {
        (x_max - x_min) / (width - 1) as f32
    } else {
        0.0
    };
    let sy = if height > 1 {
        (y_max - y_min) / (height - 1) as f32
    } else {
        0.0
    };

    let pb = ProgressBar::new(height as u64);
    pb.set_style(
        ProgressStyle::with_template(
            "{spinner:.green} [{elapsed_precise}] [{bar:40.cyan/blue}] {pos}/{len} rows ({eta})",
        )
        .unwrap()
        .progress_chars("#>-")
        .tick_chars("⠁⠂⠄⡀⢀⠠⠐⠈ "),
    );

    for y in 0..height {
        let wy = y_min + (y as f32) * sy;
        for x in 0..width {
            let wx = x_min + (x as f32) * sx;
            let p = vector![wx, wy];

            // Early-out using the speed bound; if inconclusive, compute with derivative
            let lum = if is_definitely_far(
                p,
                curve,
                t_range.clone(),
                threshold,
                speed_upper_bound,
                opts.samples,
            ) {
                255u8
            } else {
                let res = nearest_t_with_derivative(p, curve, curve_d, t_range.clone(), opts);
                let d = res.distance;
                ((d / threshold).clamp(0.0, 1.0) * 255.0).round() as u8
            };

            img.put_pixel(x, height - 1 - y, Luma([lum]));
        }
        pb.inc(1);
    }
    pb.finish_with_message("render complete");

    let out_path =
        std::env::var("RENDER_OUT").unwrap_or_else(|_| "curve_derivative.png".to_string());
    img.save(&out_path)?;
    eprintln!("Wrote {}", out_path);
    Ok(())
}