blue-noise 0.2.2

/**
 * Blue Noise Texture Generator using Void-and-Cluster Algorithm
 *
 * Implementation of Robert Ulichney's void-and-cluster method for generating
 * high-quality blue noise textures. Blue noise has evenly distributed energy
 * at high frequencies whilst minimising low-frequency content, avoiding the
 * clustering and void patterns found in white noise.
 *
 * ALGORITHM OVERVIEW
 * ==================
 *
 * Blue noise addresses issues with other dithering methods:
 * - White noise: Random distribution creates visible clusters and voids
 * - Bayer dithering: Regular patterns create repetitive artefacts
 * - Blue noise: Evenly distributed with minimal low-frequency patterns
 *
 * The void-and-cluster algorithm works by:
 * 1. Finding clusters (areas of high density) using Gaussian blur
 * 2. Finding voids (areas of low density) using the same blur
 * 3. Iteratively redistributing points to spread them evenly
 * 4. Assigning each pixel a rank based on its importance
 *
 * TORUS TOPOLOGY
 * ==============
 * All distance calculations wrap around at the edges (toroidal topology),
 * ensuring the resulting texture tiles seamlessly when repeated. This is
 * essential for dithering large images with small noise textures.
 *
 * FFT OPTIMISATION
 * ================
 * For power-of-two dimensions, Gaussian blur is performed in the frequency
 * domain using Fast Fourier Transform. Convolution becomes element-wise
 * multiplication in frequency space, providing ~50% performance improvement.
 *
 * GENERATION PHASES
 * =================
 *
 * Phase 0: Generate initial binary pattern
 *   - Place random points and redistribute until convergence
 *   - Convergence occurs when tightest cluster equals largest void
 *
 * Phase 1: Serialize initial points
 *   - Remove points from tightest clusters
 *   - Assign ranks from (initialPoints - 1) down to 0
 *
 * Phase 2: Fill to half capacity
 *   - Restore initial pattern and add points to largest voids
 *   - Assign ranks from initialPoints to area/2
 *
 * Phase 3: Fill to completion
 *   - Invert bitmap (0s become minority)
 *   - Remove minority points from tightest clusters
 *   - Assign ranks from area/2 to area-1
 *
 * Phase 4: Convert to threshold map
 *   - Map ranks [0, area-1] to threshold values [0, 255]
 *
 * REFERENCES
 * ==========
 * - Ulichney, R. (1993). "Void-and-cluster method for dither array generation"
 *   Proceedings of SPIE 1913, Human Vision, Visual Processing, and Digital
 *   Display IV. https://doi.org/10.1117/12.152707
 *
 * - Ulichney, R. (1988). "Dithering with blue noise"
 *   Proceedings of the IEEE, 76(1), 56-79.
 *
 * PERFORMANCE
 * ===========
 * 64�64 texture:  ~2-5 seconds
 * 128�128 texture: ~30-60 seconds
 * 256�256 texture: Several minutes
 *
 * For production use, pre-generate textures rather than generating at runtime.
 */
use image::{GrayImage, ImageBuffer, Luma};
use indicatif::{ProgressBar, ProgressStyle};
use rustfft::FftPlanner;
use rustfft::num_complex::Complex;
use std::path::Path;
use thiserror::Error;

/// Configuration for blue noise generation
#[derive(Debug, Clone)]
pub struct BlueNoiseConfig {
    /// Width of the texture in pixels
    pub width: usize,
    /// Height of the texture in pixels
    pub height: usize,
    /// Gaussian blur sigma value (typically 1.5-2.5)
    pub sigma: f32,
    /// Initial density of points (typically 0.1)
    pub initial_density: f32,
    /// Optional random seed for reproducible results
    pub seed: Option<u32>,
    /// Show progress indicators
    pub verbose: bool,
}

impl Default for BlueNoiseConfig {
    fn default() -> Self {
        Self {
            width: 64,
            height: 64,
            sigma: 1.9,
            initial_density: 0.1,
            seed: None,
            verbose: false,
        }
    }
}

/// Result of blue noise generation
#[derive(Debug, Clone)]
pub struct BlueNoiseResult {
    /// Grayscale threshold map data (0-255)
    pub data: Vec<u8>,
    /// Width of the generated texture
    pub width: usize,
    /// Height of the generated texture
    pub height: usize,
}

/// Error types for blue noise generation
#[derive(Error, Debug)]
pub enum GeneratorError {
    /// Width or height is zero or too large to process safely
    #[error("Width and height must be positive and small enough to process safely")]
    InvalidDimensions,

    /// Sigma parameter is zero or negative
    #[error("Sigma must be positive")]
    InvalidSigma,

    /// Initial density is out of valid range (0, 1)
    #[error("Initial density must be between 0 and 1")]
    InvalidDensity,

    /// Failed to save generated image
    #[error("Failed to save image: {0}")]
    ImageSaveError(#[from] image::ImageError),

    /// Algorithm failed to converge
    #[error("Generation failed to converge")]
    ConvergenceError,
}

/// Result type for generator operations
pub type Result<T> = std::result::Result<T, GeneratorError>;

/**
 * Mulberry32 seeded random number generator
 * Fast, high-quality PRNG for reproducible results
 */
struct SeededRandom {
    seed: u32,
}

impl SeededRandom {
    fn new(seed: Option<u32>) -> Self {
        Self {
            seed: seed.unwrap_or_else(|| {
                std::time::SystemTime::now()
                    .duration_since(std::time::UNIX_EPOCH)
                    .unwrap()
                    .as_millis() as u32
            }),
        }
    }

    fn next(&mut self) -> f32 {
        self.seed = self.seed.wrapping_add(0x6D2B79F5);
        let mut t = self.seed ^ (self.seed >> 15);
        t = t.wrapping_mul(1 | self.seed);
        t ^= t.wrapping_add(t.wrapping_mul(t ^ (t >> 7)).wrapping_mul(61 | t));
        let bits = t ^ (t >> 14);
        Self::normalize(bits)
    }

    fn normalize(bits: u32) -> f32 {
        // Map a u32 into [0, 1) using only the top 24 bits, which fit exactly
        // in an f32 mantissa. Dividing a full u32 by 2^32 as f32 can round up to
        // exactly 1.0 for the top 128 values, which then makes
        // `(next() * area) as usize == area` and panics with an out-of-bounds
        // index. 2^24 is exactly representable, so the result is always < 1.0.
        ((bits >> 8) as f32) / 16_777_216.0
    }

    fn index_below(&mut self, upper_bound: usize) -> usize {
        debug_assert!(upper_bound > 0);
        ((self.next() * upper_bound as f32) as usize).min(upper_bound - 1)
    }
}

/**
 * Main class for generating blue noise textures
 */
pub struct BlueNoiseGenerator {
    // Constants
    max_iterations_multiplier: usize,
    threshold_map_levels: usize,

    // Configuration
    width: usize,
    height: usize,
    area: usize,
    sigma: f32,
    initial_density: f32,
    verbose: bool,
    random: SeededRandom,

    // Working arrays
    bitmap: Vec<u8>,
    rank: Vec<i32>,
    energy: Vec<f32>,

    // Cached values for performance
    ones_count: usize,

    // FFT optimization
    use_fft: bool,
    gaussian_kernel: Vec<f32>,
    gaussian_kernel_freq: Option<Vec<Complex<f32>>>,

    // Progress bar
    progress: Option<ProgressBar>,
}

impl BlueNoiseGenerator {
    const MAX_ITERATIONS_MULTIPLIER: usize = 10;
    const THRESHOLD_MAP_LEVELS: usize = 256;

    /// Check if a number is a power of two
    fn is_power_of_two(n: usize) -> bool {
        n > 0 && (n & (n - 1)) == 0
    }

    /// Create a new generator with the given configuration
    pub fn new(config: BlueNoiseConfig) -> Result<Self> {
        // Validation
        if config.width == 0 || config.height == 0 {
            return Err(GeneratorError::InvalidDimensions);
        }
        if config.sigma <= 0.0 {
            return Err(GeneratorError::InvalidSigma);
        }
        if config.initial_density <= 0.0 || config.initial_density >= 1.0 {
            return Err(GeneratorError::InvalidDensity);
        }

        if config.width > u32::MAX as usize || config.height > u32::MAX as usize {
            return Err(GeneratorError::InvalidDimensions);
        }

        let area = config
            .width
            .checked_mul(config.height)
            .ok_or(GeneratorError::InvalidDimensions)?;

        if area > i32::MAX as usize || area.checked_mul(Self::MAX_ITERATIONS_MULTIPLIER).is_none() {
            return Err(GeneratorError::InvalidDimensions);
        }

        let use_fft = Self::is_power_of_two(config.width) && Self::is_power_of_two(config.height);

        let progress = if config.verbose {
            Some(ProgressBar::new(100))
        } else {
            None
        };

        if let Some(pb) = &progress {
            pb.set_style(
                ProgressStyle::default_bar()
                    .template("[{elapsed_precise}] {bar:40.cyan/blue} {pos:>3}% {msg}")
                    .unwrap()
                    .progress_chars("##-"),
            );
        }

        let mut generator = Self {
            max_iterations_multiplier: Self::MAX_ITERATIONS_MULTIPLIER,
            threshold_map_levels: Self::THRESHOLD_MAP_LEVELS,
            width: config.width,
            height: config.height,
            area,
            sigma: config.sigma,
            initial_density: config.initial_density,
            verbose: config.verbose,
            random: SeededRandom::new(config.seed),
            bitmap: vec![0; area],
            rank: vec![0; area],
            energy: vec![0.0; area],
            ones_count: 0,
            use_fft,
            gaussian_kernel: Vec::new(),
            gaussian_kernel_freq: None,
            progress,
        };

        generator.gaussian_kernel = generator.create_gaussian_kernel();

        if use_fft {
            generator.gaussian_kernel_freq =
                Some(generator.fft_2d_forward(&generator.gaussian_kernel));
        }

        Ok(generator)
    }

    /**
     * Create Gaussian kernel in frequency domain for FFT convolution
     *
     * Pre-computes the Gaussian kernel and transforms it to frequency space.
     * The kernel uses toroidal distance to ensure seamless tiling. By keeping
     * the kernel in frequency space, we only need to compute FFT once during
     * initialization rather than for every blur operation.
     */
    fn create_gaussian_kernel(&self) -> Vec<f32> {
        let mut kernel = vec![0.0f32; self.area];
        let divisor = 2.0 * self.sigma * self.sigma;

        for y in 0..self.height {
            for x in 0..self.width {
                // Compute distance from center with wrapping
                let dx = x.min(self.width - x) as f32;
                let dy = y.min(self.height - y) as f32;
                let dist_sq = dx * dx + dy * dy;

                kernel[y * self.width + x] = (-dist_sq / divisor).exp();
            }
        }

        // Normalize kernel
        let sum: f32 = kernel.iter().sum();
        for val in kernel.iter_mut() {
            *val /= sum;
        }

        kernel
    }

    /**
     * Perform 2D FFT on real-valued data
     */
    fn fft_2d_forward(&self, data: &[f32]) -> Vec<Complex<f32>> {
        let mut complex_data: Vec<Complex<f32>> =
            data.iter().map(|&x| Complex::new(x, 0.0)).collect();

        // FFT on rows
        let mut planner = FftPlanner::new();
        let fft = planner.plan_fft_forward(self.width);

        for y in 0..self.height {
            let start = y * self.width;
            let end = start + self.width;
            fft.process(&mut complex_data[start..end]);
        }

        // FFT on columns
        let fft = planner.plan_fft_forward(self.height);
        let mut column = vec![Complex::new(0.0, 0.0); self.height];

        for x in 0..self.width {
            for y in 0..self.height {
                column[y] = complex_data[y * self.width + x];
            }
            fft.process(&mut column);
            for y in 0..self.height {
                complex_data[y * self.width + x] = column[y];
            }
        }

        complex_data
    }

    /**
     * Perform 2D inverse FFT
     */
    fn fft_2d_inverse(&self, complex_data: &[Complex<f32>]) -> Vec<f32> {
        let mut data = complex_data.to_vec();

        // Inverse FFT on columns
        let mut planner = FftPlanner::new();
        let ifft = planner.plan_fft_inverse(self.height);
        let mut column = vec![Complex::new(0.0, 0.0); self.height];

        for x in 0..self.width {
            for y in 0..self.height {
                column[y] = data[y * self.width + x];
            }
            ifft.process(&mut column);
            for y in 0..self.height {
                data[y * self.width + x] = column[y];
            }
        }

        // Inverse FFT on rows
        let ifft = planner.plan_fft_inverse(self.width);

        for y in 0..self.height {
            let start = y * self.width;
            let end = start + self.width;
            ifft.process(&mut data[start..end]);
        }

        // Extract real parts and normalize
        data.iter().map(|c| c.re / (self.area as f32)).collect()
    }

    /**
     * Apply Gaussian blur using FFT (frequency domain convolution)
     *
     * FFT optimization: Convolution in the spatial domain (O(n� � k�) where k is
     * kernel size) becomes element-wise multiplication in the frequency domain
     * (O(n� log n) for FFT). This provides ~50% performance improvement for
     * power-of-two dimensions.
     *
     * The convolution theorem states: convolution(A, B) = IFFT(FFT(A) � FFT(B))
     */
    fn gaussian_blur_fft(&self, data: &[u8]) -> Vec<f32> {
        // Convert to float and transform to frequency domain
        let float_data: Vec<f32> = data.iter().map(|&x| x as f32).collect();
        let data_freq = self.fft_2d_forward(&float_data);

        // Element-wise multiplication in frequency domain (convolution)
        let kernel_freq = self.gaussian_kernel_freq.as_ref().unwrap();
        let result_freq: Vec<Complex<f32>> = data_freq
            .iter()
            .zip(kernel_freq.iter())
            .map(|(d, k)| d * k)
            .collect();

        // Transform back to spatial domain
        self.fft_2d_inverse(&result_freq)
    }

    /**
     * Apply Gaussian blur using spatial domain convolution (fallback)
     *
     * Used when dimensions are not powers of two. Applies the Gaussian kernel
     * directly in the spatial domain by computing weighted sums of neighbouring
     * pixels. Coordinates wrap around (torus topology) to ensure seamless tiling.
     */
    fn gaussian_blur_spatial(&self, data: &[u8]) -> Vec<f32> {
        let mut blurred = vec![0.0f32; self.area];

        for (source_idx, &source) in data.iter().enumerate().take(self.area) {
            let value = source as f32;
            if value == 0.0 {
                continue;
            }

            for (target_idx, target) in blurred.iter_mut().enumerate() {
                let kernel_idx = self.kernel_index_for_delta(source_idx, target_idx);
                *target += value * self.gaussian_kernel[kernel_idx];
            }
        }

        blurred
    }

    /**
     * Apply Gaussian blur (chooses FFT or spatial based on size)
     */
    fn gaussian_blur(&self, data: &[u8]) -> Vec<f32> {
        if self.use_fft {
            self.gaussian_blur_fft(data)
        } else {
            self.gaussian_blur_spatial(data)
        }
    }

    /**
     * Find the tightest cluster: pixel with highest energy among all 1s
     *
     * A "cluster" is a region where pixels are densely packed. The Gaussian
     * blur creates an energy field where areas with many nearby 1s have high
     * energy. The tightest cluster is the 1-pixel with the most neighbours.
     */
    fn find_tightest_cluster(&self) -> Option<usize> {
        let mut max_energy = f32::NEG_INFINITY;
        let mut max_idx = None;

        for i in 0..self.area {
            if self.bitmap[i] == 1 && self.energy[i] > max_energy {
                max_energy = self.energy[i];
                max_idx = Some(i);
            }
        }

        max_idx
    }

    /**
     * Find the largest void: pixel with lowest energy among all 0s
     *
     * A "void" is a region where pixels are sparse. In the energy field,
     * areas with few nearby 1s have low energy. The largest void is the
     * 0-pixel with the fewest neighbours.
     */
    fn find_largest_void(&self) -> Option<usize> {
        let mut min_energy = f32::INFINITY;
        let mut min_idx = None;

        for i in 0..self.area {
            if self.bitmap[i] == 0 && self.energy[i] < min_energy {
                min_energy = self.energy[i];
                min_idx = Some(i);
            }
        }

        min_idx
    }

    /**
     * Count number of 1s in the bitmap (returns cached value)
     */
    fn count_ones(&self) -> usize {
        self.ones_count
    }

    /**
     * Set a bit in the bitmap and update the cache
     */
    fn set_bit(&mut self, idx: usize, value: u8) {
        let old_value = self.bitmap[idx];
        if old_value == value {
            return;
        }

        self.bitmap[idx] = value;
        self.ones_count = self.ones_count + value as usize - old_value as usize;
    }

    fn set_bit_and_update_energy(&mut self, idx: usize, value: u8) {
        let old_value = self.bitmap[idx];
        if old_value == value {
            return;
        }

        self.set_bit(idx, value);
        self.apply_energy_delta(idx, value as f32 - old_value as f32);
    }

    fn kernel_index_for_delta(&self, source_idx: usize, target_idx: usize) -> usize {
        let source_x = source_idx % self.width;
        let source_y = source_idx / self.width;
        let target_x = target_idx % self.width;
        let target_y = target_idx / self.width;

        let dx = if target_x >= source_x {
            target_x - source_x
        } else {
            self.width - source_x + target_x
        };
        let dy = if target_y >= source_y {
            target_y - source_y
        } else {
            self.height - source_y + target_y
        };

        dy * self.width + dx
    }

    fn apply_energy_delta(&mut self, idx: usize, delta: f32) {
        for target_idx in 0..self.area {
            let kernel_idx = self.kernel_index_for_delta(idx, target_idx);
            self.energy[target_idx] += delta * self.gaussian_kernel[kernel_idx];
        }
    }

    /**
     * Recalculate ones count from scratch (used for initialization)
     */
    fn recalculate_ones_count(&mut self) {
        self.ones_count = self.bitmap.iter().map(|&x| x as usize).sum();
    }

    /**
     * Recalculate energy field by applying Gaussian blur to bitmap
     *
     * The energy field is the key to finding clusters and voids. By applying
     * a Gaussian blur to the binary pattern, we create a smooth field where:
     * - High values indicate clusters (many nearby 1s)
     * - Low values indicate voids (few nearby 1s)
     *
     * The blur uses torus wrapping so edges connect seamlessly.
     */
    fn recalculate_energy(&mut self) {
        self.energy = self.gaussian_blur(&self.bitmap);
    }

    /**
     * Phase 0: Generate Initial Binary Pattern
     *
     * Creates a well-distributed set of initial points by:
     * 1. Randomly placing points (10% density by default)
     * 2. Iteratively swapping clustered points with void positions
     * 3. Stopping when convergence is reached (cluster = void position)
     *
     * This phase establishes the foundation for even distribution.
     */
    fn phase0_generate_initial_pattern(&mut self) -> Result<()> {
        if self.verbose
            && let Some(pb) = &self.progress
        {
            pb.set_message("Phase 0: Generating initial pattern");
            pb.set_position(0);
        }

        let target_points = (self.area as f32 * self.initial_density) as usize;

        // Randomly place initial points
        while self.count_ones() < target_points {
            let idx = self.random.index_below(self.area);
            if self.bitmap[idx] == 0 {
                self.set_bit(idx, 1);
            }
        }

        self.recalculate_energy();

        // Redistribute points until convergence
        let max_iterations = self
            .area
            .checked_mul(self.max_iterations_multiplier)
            .ok_or(GeneratorError::InvalidDimensions)?;
        let mut iterations = 0;

        while iterations < max_iterations {
            iterations += 1;

            // Find tightest cluster and remove it
            let cluster_idx = self
                .find_tightest_cluster()
                .ok_or(GeneratorError::ConvergenceError)?;
            self.set_bit_and_update_energy(cluster_idx, 0);

            // Find largest void with updated energy
            let void_idx = self
                .find_largest_void()
                .ok_or(GeneratorError::ConvergenceError)?;

            // Check for convergence
            if void_idx == cluster_idx {
                self.set_bit_and_update_energy(cluster_idx, 1);
                break;
            }

            // Place point in void
            self.set_bit_and_update_energy(void_idx, 1);
        }

        if self.verbose
            && let Some(pb) = &self.progress
        {
            pb.set_position(20);
        }

        Ok(())
    }

    /**
     * Phase 1: Serialize Initial Points
     *
     * Assigns ranks to the initial minority pattern by removing points from
     * tightest clusters first. These points are ranked from (initialPoints - 1)
     * down to 0, establishing which pixels are most important for creating
     * the blue noise distribution.
     */
    fn phase1_serialize_initial_points(&mut self) -> Result<()> {
        if self.verbose
            && let Some(pb) = &self.progress
        {
            pb.set_message("Phase 1: Serializing initial points");
            pb.set_position(20);
        }

        let mut rank_counter = self.count_ones() as i32 - 1;

        while self.count_ones() > 0 {
            let cluster_idx = self
                .find_tightest_cluster()
                .ok_or(GeneratorError::ConvergenceError)?;
            self.rank[cluster_idx] = rank_counter;
            rank_counter -= 1;

            self.set_bit_and_update_energy(cluster_idx, 0);
        }

        if self.verbose
            && let Some(pb) = &self.progress
        {
            pb.set_position(40);
        }

        Ok(())
    }

    /**
     * Phase 2: Fill to Half Capacity
     *
     * Restores the initial pattern and continues adding points to the largest
     * voids until the bitmap is 50% full. Ranks continue from initialPoints
     * to area/2. This builds up a minority pattern (less than half full).
     */
    fn phase2_fill_to_half(&mut self, prototype: &[u8], initial_points: usize) -> Result<()> {
        if self.verbose
            && let Some(pb) = &self.progress
        {
            pb.set_message("Phase 2: Filling to half capacity");
            pb.set_position(40);
        }

        self.bitmap.copy_from_slice(prototype);
        self.recalculate_ones_count();
        self.recalculate_energy();

        let mut rank_counter = initial_points as i32;
        let half_area = self.area / 2;

        while self.count_ones() < half_area {
            let void_idx = self
                .find_largest_void()
                .ok_or(GeneratorError::ConvergenceError)?;
            self.rank[void_idx] = rank_counter;
            rank_counter += 1;

            self.set_bit_and_update_energy(void_idx, 1);
        }

        if self.verbose
            && let Some(pb) = &self.progress
        {
            pb.set_position(60);
        }

        Ok(())
    }

    /**
     * Phase 3: Fill to Completion
     *
     * Inverts the bitmap so 0s become the minority (more than half full).
     * Then removes the remaining minority pixels from their tightest clusters,
     * ranking them from area/2 to area-1. This clever inversion allows the
     * algorithm to work symmetrically for both minority and majority patterns.
     */
    fn phase3_fill_to_completion(&mut self, mut rank_counter: i32) -> Result<()> {
        if self.verbose
            && let Some(pb) = &self.progress
        {
            pb.set_message("Phase 3: Filling to completion");
            pb.set_position(60);
        }

        // Invert bitmap
        for i in 0..self.area {
            self.bitmap[i] = 1 - self.bitmap[i];
        }
        self.recalculate_ones_count();
        self.recalculate_energy();

        while rank_counter < self.area as i32 {
            let cluster_idx = self
                .find_tightest_cluster()
                .ok_or(GeneratorError::ConvergenceError)?;
            self.rank[cluster_idx] = rank_counter;
            rank_counter += 1;

            self.set_bit_and_update_energy(cluster_idx, 0);
        }

        if self.verbose
            && let Some(pb) = &self.progress
        {
            pb.set_position(80);
        }

        Ok(())
    }

    /**
     * Phase 4: Convert Ranks to Threshold Map
     *
     * Maps the rank values (0 to area-1) to grayscale threshold values (0-255).
     * Each pixel's rank determines its threshold value for dithering. Lower
     * ranked pixels will be turned "on" first when dithering bright images.
     */
    fn phase4_convert_to_threshold_map(&self) -> Vec<u8> {
        if self.verbose
            && let Some(pb) = &self.progress
        {
            pb.set_message("Phase 4: Converting to threshold map");
            pb.set_position(80);
        }

        let output: Vec<u8> = self
            .rank
            .iter()
            .map(|&r| ((r as usize * self.threshold_map_levels) / self.area) as u8)
            .collect();

        if self.verbose
            && let Some(pb) = &self.progress
        {
            pb.set_position(100);
            pb.finish_with_message("Blue noise generation complete");
        }

        output
    }

    /**
     * Generate the blue noise texture
     */
    pub fn generate(mut self) -> Result<BlueNoiseResult> {
        let start_time = std::time::Instant::now();

        if self.verbose {
            println!(
                "Generating {}�{} blue noise texture...",
                self.width, self.height
            );
            println!(
                "Using {} Gaussian blur",
                if self.use_fft {
                    "FFT-optimized"
                } else {
                    "spatial"
                }
            );
        }

        self.phase0_generate_initial_pattern()?;

        let prototype = self.bitmap.clone();
        let initial_points = self.count_ones();

        if self.verbose {
            println!("Initial pattern: {} points", initial_points);
        }

        self.phase1_serialize_initial_points()?;
        self.phase2_fill_to_half(&prototype, initial_points)?;

        let half_area = self.area / 2;
        self.phase3_fill_to_completion(half_area as i32)?;

        let data = self.phase4_convert_to_threshold_map();

        if self.verbose {
            let elapsed = start_time.elapsed();
            println!(
                "Blue noise generation complete in {:.2}s",
                elapsed.as_secs_f32()
            );
        }

        Ok(BlueNoiseResult {
            data,
            width: self.width,
            height: self.height,
        })
    }
}

/**
 * Save blue noise texture to PNG file
 */
pub fn save_blue_noise_to_png<P: AsRef<Path>>(result: &BlueNoiseResult, filename: P) -> Result<()> {
    let img: GrayImage = ImageBuffer::from_fn(result.width as u32, result.height as u32, |x, y| {
        let idx = y as usize * result.width + x as usize;
        Luma([result.data[idx]])
    });

    img.save(&filename)?;
    println!(
        "Saved blue noise texture to {}",
        filename.as_ref().display()
    );

    Ok(())
}

#[cfg(test)]
mod tests {
    use super::*;

    #[test]
    fn test_seeded_random_deterministic() {
        let mut rng1 = SeededRandom::new(Some(42));
        let mut rng2 = SeededRandom::new(Some(42));

        // Same seed should produce same sequence
        for _ in 0..100 {
            assert_eq!(rng1.next(), rng2.next());
        }
    }

    #[test]
    fn test_seeded_random_range() {
        let mut rng = SeededRandom::new(Some(12345));

        // All values should be in [0, 1), and crucially strictly < 1.0 so that
        // `(next() * area) as usize` can never equal `area` (see issue #3).
        for _ in 0..5_000_000 {
            let val = rng.next();
            assert!((0.0..1.0).contains(&val), "out of range: {val}");
        }

        // The bug was that `next()` cannot return exactly 1.0. `next()` draws
        // from internal state, so we can't feed it the boundary `u32` values
        // directly. Instead, verify the conversion formula `next()` uses is
        // safe for the worst-case inputs: the top u32 values are exactly the
        // ones that previously rounded up to 1.0f32 when divided by 2^32.
        for bits in [
            u32::MAX,
            u32::MAX - 1,
            u32::MAX - 127,
            4_294_967_168, // u32::MAX - 127, the first value that used to overflow to 1.0
            0,
            1,
        ] {
            let val = SeededRandom::normalize(bits);
            assert!((0.0..1.0).contains(&val), "boundary {bits} -> {val}");
        }
        // The single largest possible output must still be < 1.0.
        assert!(SeededRandom::normalize(u32::MAX) < 1.0);
    }

    #[test]
    fn test_seeded_random_index_below_never_returns_upper_bound() {
        let mut rng = SeededRandom::new(Some(12345));

        assert_eq!(rng.index_below(1), 0);

        for _ in 0..1000 {
            assert!(rng.index_below(262_144) < 262_144);
        }
    }

    #[test]
    fn test_is_power_of_two() {
        assert!(BlueNoiseGenerator::is_power_of_two(1));
        assert!(BlueNoiseGenerator::is_power_of_two(2));
        assert!(BlueNoiseGenerator::is_power_of_two(4));
        assert!(BlueNoiseGenerator::is_power_of_two(8));
        assert!(BlueNoiseGenerator::is_power_of_two(16));
        assert!(BlueNoiseGenerator::is_power_of_two(64));
        assert!(BlueNoiseGenerator::is_power_of_two(128));
        assert!(BlueNoiseGenerator::is_power_of_two(256));

        assert!(!BlueNoiseGenerator::is_power_of_two(0));
        assert!(!BlueNoiseGenerator::is_power_of_two(3));
        assert!(!BlueNoiseGenerator::is_power_of_two(5));
        assert!(!BlueNoiseGenerator::is_power_of_two(100));
        assert!(!BlueNoiseGenerator::is_power_of_two(255));
    }

    #[test]
    fn test_config_validation() {
        // Valid config should work
        let config = BlueNoiseConfig {
            width: 64,
            height: 64,
            sigma: 1.9,
            initial_density: 0.1,
            seed: Some(42),
            verbose: false,
        };
        assert!(BlueNoiseGenerator::new(config).is_ok());

        // Zero width should fail
        let config = BlueNoiseConfig {
            width: 0,
            height: 64,
            ..Default::default()
        };
        assert!(BlueNoiseGenerator::new(config).is_err());

        // Zero height should fail
        let config = BlueNoiseConfig {
            width: 64,
            height: 0,
            ..Default::default()
        };
        assert!(BlueNoiseGenerator::new(config).is_err());

        // Negative sigma should fail
        let config = BlueNoiseConfig {
            width: 64,
            height: 64,
            sigma: -1.0,
            ..Default::default()
        };
        assert!(BlueNoiseGenerator::new(config).is_err());

        // Invalid density (too low) should fail
        let config = BlueNoiseConfig {
            width: 64,
            height: 64,
            initial_density: 0.0,
            ..Default::default()
        };
        assert!(BlueNoiseGenerator::new(config).is_err());

        // Invalid density (too high) should fail
        let config = BlueNoiseConfig {
            width: 64,
            height: 64,
            initial_density: 1.0,
            ..Default::default()
        };
        assert!(BlueNoiseGenerator::new(config).is_err());

        // Overflowing area should fail before allocation
        let config = BlueNoiseConfig {
            width: usize::MAX,
            height: 2,
            ..Default::default()
        };
        assert!(BlueNoiseGenerator::new(config).is_err());

        // Rank counters use i32 internally, so larger areas are rejected
        let config = BlueNoiseConfig {
            width: i32::MAX as usize + 1,
            height: 1,
            ..Default::default()
        };
        assert!(BlueNoiseGenerator::new(config).is_err());
    }

    #[test]
    fn test_generate_small_texture() {
        let config = BlueNoiseConfig {
            width: 16,
            height: 16,
            sigma: 1.5,
            seed: Some(42),
            verbose: false,
            ..Default::default()
        };

        let generator = BlueNoiseGenerator::new(config).unwrap();
        let result = generator.generate().unwrap();

        assert_eq!(result.width, 16);
        assert_eq!(result.height, 16);
        assert_eq!(result.data.len(), 256);

        assert!(result.data.iter().any(|&val| val > 0));
        assert!(result.data.iter().any(|&val| val < 255));
    }

    #[test]
    fn test_incremental_energy_matches_full_recalculation() {
        let config = BlueNoiseConfig {
            width: 8,
            height: 8,
            sigma: 1.5,
            seed: Some(42),
            verbose: false,
            ..Default::default()
        };

        let mut generator = BlueNoiseGenerator::new(config).unwrap();
        generator.set_bit(0, 1);
        generator.recalculate_energy();

        generator.set_bit_and_update_energy(17, 1);
        let incremental = generator.energy.clone();

        generator.recalculate_energy();

        for (incremental, recalculated) in incremental.iter().zip(generator.energy.iter()) {
            assert!((incremental - recalculated).abs() < 0.000_01);
        }
    }

    #[test]
    fn test_generate_reproducible() {
        let config1 = BlueNoiseConfig {
            width: 32,
            height: 32,
            seed: Some(12345),
            verbose: false,
            ..Default::default()
        };

        let config2 = BlueNoiseConfig {
            width: 32,
            height: 32,
            seed: Some(12345),
            verbose: false,
            ..Default::default()
        };

        let gen1 = BlueNoiseGenerator::new(config1).unwrap();
        let gen2 = BlueNoiseGenerator::new(config2).unwrap();

        let result1 = gen1.generate().unwrap();
        let result2 = gen2.generate().unwrap();

        // Same seed should produce identical results
        assert_eq!(result1.data, result2.data);
    }

    #[test]
    fn test_generate_different_seeds() {
        let config1 = BlueNoiseConfig {
            width: 32,
            height: 32,
            seed: Some(111),
            verbose: false,
            ..Default::default()
        };

        let config2 = BlueNoiseConfig {
            width: 32,
            height: 32,
            seed: Some(222),
            verbose: false,
            ..Default::default()
        };

        let gen1 = BlueNoiseGenerator::new(config1).unwrap();
        let gen2 = BlueNoiseGenerator::new(config2).unwrap();

        let result1 = gen1.generate().unwrap();
        let result2 = gen2.generate().unwrap();

        // Different seeds should produce different results
        assert_ne!(result1.data, result2.data);
    }

    #[test]
    fn test_generate_power_of_two_uses_fft() {
        let config = BlueNoiseConfig {
            width: 64,
            height: 64,
            seed: Some(42),
            verbose: false,
            ..Default::default()
        };

        let generator = BlueNoiseGenerator::new(config).unwrap();
        assert!(generator.use_fft);
        assert!(generator.gaussian_kernel_freq.is_some());
    }

    #[test]
    fn test_generate_non_power_of_two_no_fft() {
        let config = BlueNoiseConfig {
            width: 50,
            height: 50,
            seed: Some(42),
            verbose: false,
            ..Default::default()
        };

        let generator = BlueNoiseGenerator::new(config).unwrap();
        assert!(!generator.use_fft);
        assert!(generator.gaussian_kernel_freq.is_none());
    }

    #[test]
    fn test_threshold_map_distribution() {
        let config = BlueNoiseConfig {
            width: 32,
            height: 32,
            seed: Some(99),
            verbose: false,
            ..Default::default()
        };

        let generator = BlueNoiseGenerator::new(config).unwrap();
        let result = generator.generate().unwrap();

        // Check that we have good distribution across the range
        let mut histogram = vec![0usize; 256];
        for &val in &result.data {
            histogram[val as usize] += 1;
        }

        // Should have values across the spectrum
        let non_empty_bins = histogram.iter().filter(|&&count| count > 0).count();
        assert!(non_empty_bins > 200, "Expected diverse distribution");
    }

    #[test]
    fn test_rectangular_texture() {
        let config = BlueNoiseConfig {
            width: 32,
            height: 16,
            seed: Some(42),
            verbose: false,
            ..Default::default()
        };

        let generator = BlueNoiseGenerator::new(config).unwrap();
        let result = generator.generate().unwrap();

        assert_eq!(result.width, 32);
        assert_eq!(result.height, 16);
        assert_eq!(result.data.len(), 512);
    }

    #[test]
    fn test_default_config() {
        let config = BlueNoiseConfig::default();
        assert_eq!(config.width, 64);
        assert_eq!(config.height, 64);
        assert_eq!(config.sigma, 1.9);
        assert_eq!(config.initial_density, 0.1);
        assert_eq!(config.seed, None);
        assert!(!config.verbose);
    }
}