wallswitch 0.60.12

//! Procedural wallpaper overlay common utilities and math structures.
//!
//! This module provides shared helper functions, coordinate system transformations,
//! blending equations, and escape-time evaluation loops used across all fractal
//! and wave-based procedural generation engines.

use crate::{
    AuroraGenerator, ColorRGB, Complex, JuliaGenerator, MandelbrotGenerator, Monitor, NeonColor,
    NewtonGenerator, NovaGenerator, StarfieldGenerator, WallSwitchError, WallSwitchResult,
    get_random_integer,
};
use clap::ValueEnum;
use image::RgbImage;
use rayon::prelude::*;
use serde::{Deserialize, Serialize};
use std::{
    f64::consts::{LOG2_E, PI},
    io::Error,
    path::Path,
};

/// The default minimum iteration limit for escape-time fractal calculation.
pub const MIN_ITERATIONS: u32 = 800;

/// The default maximum iteration limit for escape-time fractal calculation.
pub const MAX_ITERATIONS: u32 = 1200;

/// The number of angular steps used to evaluate structural rotations during optimization.
pub const ROTATION_STEPS: usize = 16;

/// Trait defining the behavior for any image processing effect.
///
/// This allows different procedural generators to be treated polymorphically,
/// keeping the rendering engine decoupled from specific math implementations.
pub trait ImageEffect: Sync + Send {
    /// Applies the procedural effect in-place to a mutable image buffer.
    fn apply(&self, rgb_img: &mut RgbImage);

    /// Returns a formatted string containing specific diagnostic details of the active effect.
    fn info(&self) -> String;

    /// Helper to process, render, and write output files directly to system storage.
    fn apply_effect(&self, input_path: &Path, output_path: &Path) -> WallSwitchResult<()> {
        let img = image::open(input_path)
            .map_err(|e| WallSwitchError::UnableToFind(format!("Failed to open image: {e}")))?;

        let mut rgb_img = img.to_rgb8();
        self.apply(&mut rgb_img);

        rgb_img
            .save(output_path)
            .map_err(|e| WallSwitchError::Io(Error::other(e)))?;

        Ok(())
    }
}

/// Unified generic configuration for procedural fractal viewports and rendering properties.
///
/// This structure centralizes shared parameters required to map and color escape-time
/// or convergence-based fractals, including iteration boundaries, visual color palettes,
/// and spatial transform variables (zoom scale and rotational phasor).
#[derive(Debug, Clone)]
pub struct FractalConfig {
    /// The maximum iteration limit for escape-time calculations.
    pub scan_iterations: u32,
    /// The base color palette selected for the neon glow.
    pub color_palette: NeonColor,
    /// The viewport zoom scale level.
    pub zoom: f64,
    /// The complex phasor representing the viewport rotation.
    pub rotation: Complex,
}

/// A polymorphic trait that defines the core algebraic structure for any procedural fractal.
/// Implementing this trait automatically provides an optimized ImageEffect implementation.
pub trait FractalDescriptor {
    /// Retrieves the shared viewport layout configuration.
    fn config(&self) -> &FractalConfig;
    /// Focus center point on the complex plane.
    fn center(&self) -> Complex;
    /// Toggles the origin mapping behavior used by Julia variants.
    fn is_julia(&self) -> bool;
    /// Translates a localized complex pixel coordinate into its blended color result.
    fn render_pixel(&self, z_init: Complex, scale: f64, max_radius: f64) -> (ColorRGB, f64, f64);
    /// Returns a comprehensive diagnostic string formatted for the fractal equation.
    fn info_text(&self) -> String;
}

// Blanket implementation enforcing DRY architectural principles across all generated fractals.
impl<T: FractalDescriptor + Sync + Send> ImageEffect for T {
    fn apply(&self, rgb_img: &mut RgbImage) {
        let cfg = self.config();
        render_fractal_parallel(
            rgb_img,
            cfg.zoom,
            cfg.rotation,
            self.center(),
            self.is_julia(),
            |z, scale, max_radius| self.render_pixel(z, scale, max_radius),
        );
    }
    fn info(&self) -> String {
        self.info_text()
    }
}

/// Represents the supported procedural background overlay effects.
#[derive(Default, Debug, Clone, Copy, PartialEq, Eq, Serialize, Deserialize, ValueEnum)]
#[serde(rename_all = "lowercase")]
pub enum ProceduralEffect {
    /// No overlay effect is applied; displays the raw, unaltered wallpaper.
    #[value(name = "none")]
    #[default]
    None,

    /// Julia Set fractal overlay.
    ///
    /// * Characteristics: Rendered as thin, sharp, self-similar contour lines forming highly
    ///   symmetrical branching patterns. Depending on the selected complex constant, the lines trace
    ///   intricate shapes resembling swirling clouds, dendritic lace, spiral galaxy arms, leafy
    ///   filaments, or crystalline snowflakes.
    /// * Creator: Developed mathematically by the French mathematician Gaston Julia in 1918.
    /// * Generator function: Calculated by mapping the convergence boundary under the recursive function:
    ///   f(z) = z^2 + c
    ///   where c is a fixed complex constant perturbation and the initial coordinate z_0 varies across the viewport.
    #[value(name = "julia")]
    JuliaSet,

    /// Mandelbrot Set fractal overlay.
    ///
    /// * Characteristics: Rendered as thin, sharp, self-similar contour lines tracing the boundary of
    ///   the set. The lines expose highly detailed structural contours, including a main cardioid,
    ///   circular period bulbs, swirling spiral valleys, and repeating miniature copies of
    ///   the entire set connected by thin filaments.
    /// * Creator: First visualized and defined by the Polish-born French-American mathematician
    ///   Benoit Mandelbrot in 1980.
    /// * Generator function: Modeled using the quadratic recurrence equation starting from the origin:
    ///   z(n+1) = z(n)^2 + c
    ///   where z_0 = 0 and the complex parameter c varies across the viewport grid coordinates.
    #[value(name = "mandelbrot")]
    Mandelbrot,

    /// Newton-Raphson Basin of Attraction fractal overlay.
    ///
    /// * Characteristics: Symmetrical, kaleidoscope-like mandala structures representing root-finding
    ///   convergence fields across complex space boundaries. It maps the limits of convergence zones
    ///   where points migrate to specific roots of a polynomial equation.
    /// * Creator: Formulated based on Sir Isaac Newton's root-approximation methods (1690s) and Arthur
    ///   Cayley's subsequent complex-plane studies (1879).
    /// * Generator function: Computed using a relaxed Newton-Raphson recurrence formula:
    ///   z(n+1) = z(n) - lambda * f(z(n)) / f'(z(n))
    ///   on the polynomial f(z) = z^p - 1, where p is the integer polynomial power and lambda is a complex relaxation factor.
    #[value(name = "newton")]
    NewtonBasins,

    /// Nova Julia liquid fractal overlay.
    ///
    /// * Characteristics: Organic, flowing, fluid-like plumes resembling liquid mercury,
    ///   cosmic nebulae, or dynamic plasma current paths.
    /// * Creator: Developed by Paul Derbyshire in the late 1990s as a structural variation and
    ///   relaxation of the classic Newton-Raphson fractal.
    /// * Generator function: Evaluated using the relaxed Newton recurrence relation perturbed by a
    ///   dynamic additive complex value:
    ///   z(n+1) = z(n) - R * (z(n)^p - 1) / (p * z(n)^(p-1)) + c
    ///   where p is the polynomial exponent, R is a complex relaxation modifier, and c is a fixed perturbation coordinate.
    #[value(name = "nova")]
    NovaJulia,

    /// Procedural Cosmic Aurora wave generator.
    ///
    /// Generates glowing atmospheric wave filaments using multi-frequency wave mathematics.
    #[value(name = "aurora")]
    CosmicAurora,

    /// Procedural Starfield / Bokeh generator.
    ///
    /// Projects glowing, circular stars and light orbs with smooth Gaussian light falloffs.
    #[value(name = "star")]
    Starfield,

    /// Fractal mode selector.
    ///
    /// Randomly selects between Julia, Mandelbrot, Newton-Raphson Basins, or Nova Julia fractal overlays.
    #[value(name = "fractal")]
    Fractal,

    /// Fully randomized mode selector.
    ///
    /// Automatically decides on a random overlay effect independently for each physical display.
    #[value(name = "random")]
    Random,
}

impl ProceduralEffect {
    pub fn get_name(self) -> &'static str {
        match self {
            Self::None => "None",
            Self::JuliaSet => "Julia Sets",
            Self::Mandelbrot => "Mandelbrot",
            Self::NewtonBasins => "Newton Basins",
            Self::NovaJulia => "Nova Julia",
            Self::CosmicAurora => "Cosmic Aurora",
            Self::Starfield => "Starfield",
            Self::Fractal => "Fractal",
            Self::Random => "Random",
        }
    }

    pub fn resolve(self) -> Self {
        match self {
            Self::Random => match get_random_integer(0, 5) {
                0 => Self::JuliaSet,
                1 => Self::Mandelbrot,
                2 => Self::NewtonBasins,
                3 => Self::NovaJulia,
                4 => Self::CosmicAurora,
                _ => Self::Starfield,
            },
            Self::Fractal => match get_random_integer(0, 3) {
                0 => Self::JuliaSet,
                1 => Self::Mandelbrot,
                2 => Self::NewtonBasins,
                _ => Self::NovaJulia,
            },
            concrete => concrete,
        }
    }

    pub fn get_renderer(self, monitor: &Monitor) -> Option<Box<dyn ImageEffect>> {
        match self {
            Self::JuliaSet => Some(Box::new(JuliaGenerator::random(monitor))),
            Self::Mandelbrot => Some(Box::new(MandelbrotGenerator::random(monitor))),
            Self::NewtonBasins => Some(Box::new(NewtonGenerator::random(monitor))),
            Self::NovaJulia => Some(Box::new(NovaGenerator::random(monitor))),
            Self::Starfield => Some(Box::new(StarfieldGenerator::random(monitor))),
            Self::CosmicAurora => Some(Box::new(AuroraGenerator::random(monitor))),
            _ => None,
        }
    }
}

/// Represents a mathematical coordinate preset for procedural fractal effects.
#[derive(Debug, Clone, Copy, PartialEq, Serialize, Deserialize)]
pub struct FractalPreset {
    /// The target focal center coordinate.
    pub center: Complex,
    /// Friendly descriptive name of the structural pattern.
    pub fractal_name: &'static str,
    /// Associated category of procedural effect.
    pub effect_name: ProceduralEffect,
}

impl std::fmt::Display for FractalPreset {
    fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
        write!(
            f,
            "{} ({:+.5} {:+.5}i) under {:?}",
            self.fractal_name, self.center.re, self.center.im, self.effect_name
        )
    }
}

/// Configuration options for optimizing relaxed fractal viewports (e.g. Newton and Nova).
#[derive(Debug, Clone, Copy)]
pub struct RelaxedViewportConfig {
    /// The horizontal physical display resolution.
    pub width: u32,
    /// The vertical physical display resolution.
    pub height: u32,
    /// Maximum search boundary in complex coordinates.
    pub search_limit: f64,
    /// Number of steps in the coordinate search grid.
    pub steps: usize,
    /// Clamping zoom boundaries.
    pub zoom_range: [f64; 2],
    /// Margin multipliers applied to randomized fitting.
    pub rand_range: [f64; 2],
    /// Fallback boundaries used when search sweeps return empty pools.
    pub fallback_range: [f64; 2],
}

/// Partitions an RGB image buffer into mutable row segments for thread-safe parallel processing.
pub fn partition_rows(rgb_img: &mut RgbImage) -> (Vec<(usize, &mut [u8])>, usize) {
    let (width, _) = rgb_img.dimensions();
    let width_usize = width as usize;
    let row_stride = width_usize * 3;
    let pixels_buffer = rgb_img.as_mut();

    let rows: Vec<(usize, &mut [u8])> = pixels_buffer
        .chunks_exact_mut(row_stride)
        .enumerate()
        .collect();

    (rows, width_usize)
}

/// Executes row-by-row processing in parallel using Rayon's work-stealing thread pool.
/// This abstracts the threading scope across effects, satisfying DRY principles.
pub fn process_rows_parallel_scoped<F>(rgb_img: &mut RgbImage, row_processor: F)
where
    F: Fn(u32, &mut [u8]) + Send + Sync,
{
    let (rows, _) = partition_rows(rgb_img);
    rows.into_par_iter().for_each(|(y, row_data)| {
        row_processor(y as u32, row_data);
    });
}

/// Applies power-law (Gamma) stretching to enhance the visual contrast of fractal filaments.
#[inline(always)]
pub fn stretch_potential(raw_t: f64) -> f64 {
    raw_t.clamp(0.0, 1.0).powf(0.35)
}

/// Calculates the continuous potential (smooth coloring) value for quadratic escape-time fractals.
///
/// Filters out low escape iterations to guarantee complete transparency in the far exterior.
#[inline]
pub fn calculate_smooth_potential(i: u32, max_iterations: u32, z: Complex) -> f64 {
    if i < max_iterations {
        let mag2 = z.abs_sq();
        let smooth_i = if mag2 > 4.0 {
            let log_zn = mag2.ln() * 0.5;
            let nu = log_zn.ln() * LOG2_E;
            (i as f64 + 1.0 - nu).max(0.0)
        } else {
            i as f64
        };

        let min_render_iter = 32.0_f64;
        if smooth_i < min_render_iter {
            return 0.0;
        }

        let normalized = (smooth_i - min_render_iter) / (max_iterations as f64 - min_render_iter);
        stretch_potential(normalized)
    } else {
        1.0
    }
}

/// Calculates the analytical distance estimator (DEM) to the boundary of the fractal set.
/// Returns a coordinate-independent distance value.
#[inline]
pub fn calculate_distance_estimator(i: u32, max_iterations: u32, z: Complex, dz: Complex) -> f64 {
    if i < max_iterations {
        let z_mag = z.abs();
        let dz_mag = dz.abs();
        if z_mag > 0.0 && dz_mag > 0.0 {
            return 2.0 * z_mag * z_mag.ln() / dz_mag;
        }
    }
    0.0
}

/// Standard smoothstep mathematical interpolation function.
#[inline(always)]
pub fn smoothstep(edge0: f64, edge1: f64, x: f64) -> f64 {
    let t = ((x - edge0) / (edge1 - edge0)).clamp(0.0, 1.0);
    t * t * (3.0 - 2.0 * t)
}

/// Blends the computed fractal color and vignette shadow onto a mutable [`ColorRGB`] pixel.
///
/// Implements component-wise graphics mathematics to cleanly blend color spaces.
#[inline(always)]
pub fn blend_and_vignette(
    pixel: &mut ColorRGB,
    fractal_rgb: ColorRGB,
    alpha: f64,
    shadow_alpha: f64,
) {
    if shadow_alpha > 0.005 {
        *pixel = pixel.scale(1.0 - shadow_alpha);
    }
    if alpha > 0.005 {
        *pixel = pixel.blend(fractal_rgb, alpha);
    }
}

/// Pure evaluation loop specifically optimized for Julia Sets, guaranteeing branchless execution in the hot loop.
#[inline(always)]
pub fn julia_escape(z_init: Complex, c: Complex, max_iter: u32) -> (u32, Complex, Complex) {
    let mut z = z_init;
    let mut dz = Complex::one();
    let mut i = 0;
    while i < max_iter {
        if z.abs_sq() > 4.0 {
            break;
        }
        dz = dz * z * 2.0;
        z = z.square() + c;
        i += 1;
    }
    (i, z, dz)
}

/// Pure evaluation loop specifically optimized for the Mandelbrot Set, incorporating cardioid and period-bulb shortcuts.
#[inline(always)]
pub fn mandelbrot_escape(c: Complex, max_iter: u32) -> (u32, Complex, Complex) {
    let q = (c - Complex::new(0.25, 0.0)).abs_sq();
    if q * (q + (c.re - 0.25)) < 0.25 * c.im * c.im {
        return (max_iter, Complex::zero(), Complex::zero());
    }
    if (c + Complex::one()).abs_sq() < 0.0625 {
        return (max_iter, Complex::zero(), Complex::zero());
    }

    let mut z = Complex::zero();
    let mut dz = Complex::zero();
    let mut i = 0;
    while i < max_iter {
        if z.abs_sq() > 4.0 {
            break;
        }
        dz = dz * z * 2.0 + Complex::one();
        z = z.square() + c;
        i += 1;
    }
    (i, z, dz)
}

/// Returns an iterator over unit complex phasors representing structural rotation angles.
#[inline(always)]
pub fn get_rotation_phasors(rotations: usize) -> impl Iterator<Item = Complex> {
    (0..rotations).map(move |step| {
        let angle = (step as f64) * 2.0 * PI / (rotations as f64);
        Complex::cis(angle)
    })
}

pub struct ViewportSpecs {
    /// Focal complex center point.
    pub center: Complex,
    /// Zoom translation scaling index.
    pub zoom: f64,
    /// Rotational coordinate transformation phasor.
    pub rotation: Complex,
    /// Toggle determining whether mapping centers relative to Julia origins.
    pub is_julia: bool,
}

pub struct Viewport {
    /// The complex coordinate representing the starting point of the viewport.
    pub start: Complex,
    /// The complex step offset increment per pixel along the screen's X-axis.
    pub dx: Complex,
    /// The complex step offset increment per pixel along the screen's Y-axis.
    pub dy: Complex,
}

impl Viewport {
    /// Creates a new coordinate viewport mapper.
    pub fn new(width: f64, height: f64, specs: &ViewportSpecs) -> Self {
        let min_dim = width.min(height);
        let scale = specs.zoom / min_dim;

        let cx_off = width / 2.0;
        let cy_off = height / 2.0;

        let dx = specs.rotation * scale;
        let dy = dx * Complex::i();

        let v_center = if specs.is_julia {
            Complex::zero()
        } else {
            specs.center
        };
        let start = v_center - dx * cx_off - dy * cy_off;

        Self { start, dx, dy }
    }

    /// Maps physical coordinate coordinates (x, y) into complex space.
    #[inline(always)]
    pub fn map(&self, x: f64, y: f64) -> Complex {
        self.start + self.dx * x + self.dy * y
    }
}

#[derive(Debug, Clone, Copy, PartialEq)]
pub struct RelaxedEscape {
    pub iterations: u32,
    pub max_iterations: u32,
    pub diff_norm: f64,
    pub z_final: Complex,
}

impl RelaxedEscape {
    /// Computes a highly consistent dual-tone visual rendering of relaxed Newton or Nova Julia boundaries.
    ///
    /// Maps converged coordinate sectors, iteration speeds, and convergence errors to a visually
    /// pleasing linear color gradient blended over the active desktop background.
    ///
    /// Uses the phase angle (argument) of the final complex coordinate to calculate lighting sectors.
    #[inline(always)]
    pub fn color(
        &self,
        color_palette: NeonColor,
        edge_fade: f64,
        ln_epsilon: f64,
        is_nova: bool,
    ) -> (ColorRGB, f64, f64) {
        if self.iterations >= self.max_iterations {
            return (ColorRGB::default(), 0.0, 0.0);
        }

        let smooth_i = self.iterations as f64 + (self.diff_norm.ln() / ln_epsilon).clamp(0.0, 1.0);

        let ripple_frequency = 0.50_f64;
        let raw_wave = (smooth_i * ripple_frequency * std::f64::consts::PI)
            .sin()
            .abs();
        let norm_dist = if is_nova {
            raw_wave.powf(2.5)
        } else {
            raw_wave
        };

        let core = if norm_dist > (if is_nova { 0.92 } else { 0.95 }) {
            (norm_dist - if is_nova { 0.92 } else { 0.95 }) / if is_nova { 0.08 } else { 0.05 }
        } else {
            0.0
        };

        let glow = if is_nova {
            norm_dist.powi(6) * 0.52
        } else {
            norm_dist.powi(5) * 0.40
        };
        let profile = if is_nova {
            core * 0.78 + glow * 0.22
        } else {
            core * 0.70 + glow * 0.30
        };
        let shadow_profile = if is_nova {
            (1.0 - norm_dist).powi(3) * 0.48
        } else {
            (1.0 - norm_dist).powi(2) * 0.35
        };

        // Leverage polar coordinates algebra to extract phase argument cleanly
        let angle = self.z_final.arg();
        let light = if is_nova {
            0.75 + 0.25 * (angle * 4.0).cos().abs()
        } else {
            0.70 + 0.30 * (angle * 3.0).cos().abs()
        };

        let t_cycled = (smooth_i * 0.08) % 1.0;
        let secondary = color_palette.rotated();
        let core_color = if is_nova {
            let t_cos = (t_cycled * PI).cos() * 0.5 + 0.5;
            secondary.lerp(color_palette.color_rgb, t_cos)
        } else {
            secondary.lerp(color_palette.color_rgb, t_cycled)
        };

        let border_color = core_color.complementary().saturate_components();
        let blended = if is_nova {
            core_color.lerp(border_color, norm_dist.powf(3.0))
        } else {
            core_color.lerp(border_color, norm_dist)
        };

        let brightness_boost = if is_nova { 1.45 } else { 1.25 };
        let rgb = blended.scale(light * brightness_boost).clamp_bounds();

        let limit_fade_iter = if is_nova { 6 } else { 8 };
        let iteration_fade = if self.iterations < limit_fade_iter {
            self.iterations as f64 / limit_fade_iter as f64
        } else {
            1.0
        };

        (
            rgb,
            profile * 0.95 * iteration_fade * edge_fade,
            shadow_profile * iteration_fade * edge_fade,
        )
    }
}

/// Unifies distance estimator coloring logic shared by Julia and Mandelbrot generators.
#[inline(always)]
pub fn color_distance_estimator(
    i: u32,
    scan_iterations: u32,
    z: Complex,
    dz: Complex,
    scale: f64,
    color_palette: NeonColor,
) -> (ColorRGB, f64, f64) {
    let t = calculate_smooth_potential(i, scan_iterations, z);
    if t <= 0.005 || i >= scan_iterations {
        return (ColorRGB::default(), 0.0, 0.0);
    }

    let dist_complex = calculate_distance_estimator(i, scan_iterations, z, dz);
    let dist_pixels = dist_complex / scale;

    let max_radius = 5.0_f64;
    let shadow_radius = max_radius * 1.5;

    if dist_pixels >= shadow_radius {
        return (ColorRGB::default(), 0.0, 0.0);
    }

    let norm_dist = (dist_pixels / max_radius).clamp(0.0, 1.0);
    let core = if dist_pixels < 1.2 {
        1.0 - (dist_pixels / 1.2)
    } else {
        0.0
    };

    let ripple_freq = 12.0_f64;
    let ripple_wave = (t * std::f64::consts::PI * ripple_freq).sin().abs();
    let nested_detail =
        (1.0 - smoothstep(0.0, 0.4, 1.0 - ripple_wave)) * (1.0 - norm_dist).max(0.0);
    let glow = if dist_pixels < max_radius {
        (1.0 - norm_dist * norm_dist).powi(6) * 0.45
    } else {
        0.0
    };

    let profile = core * 0.65 + nested_detail * 0.20 + glow * 0.15;
    let norm_shadow = (dist_pixels / shadow_radius).clamp(0.0, 1.0);
    let shadow_profile = (1.0 - norm_shadow * norm_shadow).powi(2) * 0.35;

    // Leverage polar coordinates algebra to extract phase argument cleanly
    let angle = z.arg();
    let light = 0.65 + 0.35 * (angle * 4.0).cos().abs();
    let t_cycled = (t * 2.0) % 1.0;

    let secondary = color_palette.rotated();
    let core_color = if t_cycled < 0.5 {
        secondary.lerp(color_palette.color_rgb, t_cycled * 2.0)
    } else {
        color_palette
            .color_rgb
            .lerp(secondary, (t_cycled - 0.5) * 2.0)
    };

    let border_color = core_color.complementary().saturate_components();
    let blended = core_color.lerp(border_color, norm_dist.powi(2));
    let rgb = blended.scale(light * 1.20).clamp_bounds();
    let iteration_fade = if i < 16 { (i as f64 - 3.0) / 13.0 } else { 1.0 };

    (
        rgb,
        profile * 0.95 * iteration_fade,
        shadow_profile * iteration_fade,
    )
}

/// Helper to scan a complex coordinate grid and optimize the viewport framing for any fractal type.
pub fn optimize_fractal_viewport<F>(
    width: u32,
    height: u32,
    search_limit: f64,
    steps: usize,
    rotation: Complex,
    mut escape_check: F,
) -> (f64, Complex)
where
    F: FnMut(Complex) -> bool,
{
    let inv_steps_minus_1 = 1.0 / (steps - 1) as f64;
    let range = 2.0 * search_limit;
    let mut active_points = Vec::with_capacity(steps * steps);

    for step_y in 0..steps {
        let ry = -search_limit + (step_y as f64 * inv_steps_minus_1) * range;
        for step_x in 0..steps {
            let rx = -search_limit + (step_x as f64 * inv_steps_minus_1) * range;
            let z = Complex::new(rx, ry);
            if escape_check(z) {
                active_points.push(z);
            }
        }
    }

    if !active_points.is_empty() {
        find_optimal_framing(&active_points, width, height, rotation)
    } else {
        (f64::MAX, rotation)
    }
}

/// Unifies fitting and zoom scaling with a randomized factor for relaxed convergence fractals (Newton/Nova).
pub fn optimize_relaxed_viewport<F>(
    config: RelaxedViewportConfig,
    rotation: Complex,
    mut escape_check: F,
) -> (f64, Complex)
where
    F: FnMut(Complex) -> bool,
{
    let (best_zoom, best_rotation) = optimize_fractal_viewport(
        config.width,
        config.height,
        config.search_limit,
        config.steps,
        rotation,
        &mut escape_check,
    );

    if best_zoom < f64::MAX {
        let rand_factor = get_random_integer::<_, f64>(
            (config.rand_range[0] * 100.0) as u64,
            (config.rand_range[1] * 100.0) as u64,
        ) / 100.0;
        let zoom = (best_zoom * rand_factor).clamp(config.zoom_range[0], config.zoom_range[1]);
        (zoom, best_rotation)
    } else {
        let flat_rand = get_random_integer::<_, f64>(
            (config.fallback_range[0] * 100.0) as u64,
            (config.fallback_range[1] * 100.0) as u64,
        ) / 100.0;
        (
            flat_rand.clamp(config.zoom_range[0], config.zoom_range[1]),
            rotation,
        )
    }
}

/// Helper function that performs parallel rendering over a viewport mapping physical pixels to complex coordinates.
pub fn render_fractal_parallel<F>(
    rgb_img: &mut RgbImage,
    zoom: f64,
    rotation: Complex,
    center: Complex,
    is_julia: bool,
    pixel_fn: F,
) where
    F: Fn(Complex, f64, f64) -> (ColorRGB, f64, f64) + Send + Sync,
{
    let (width, height) = rgb_img.dimensions();
    let w_f = width as f64;
    let h_f = height as f64;
    let aspect_ratio = w_f.max(h_f) / w_f.min(h_f);

    let max_radius = 0.98 * 0.5 * zoom * aspect_ratio;

    let specs = ViewportSpecs {
        center,
        zoom,
        rotation,
        is_julia,
    };
    let viewport = Viewport::new(w_f, h_f, &specs);
    let min_dim = w_f.min(h_f);
    let scale = zoom / min_dim;

    process_rows_parallel_scoped(rgb_img, |y, row_data| {
        let y_f = y as f64;
        for (x, pixel_slice) in row_data.chunks_exact_mut(3).enumerate() {
            let x_f = x as f64;
            let z_init = viewport.map(x_f, y_f);

            let (fractal_rgb, alpha, s_alpha) = pixel_fn(z_init, scale, max_radius);

            let mut pixel_color = ColorRGB::from_slice(pixel_slice);
            blend_and_vignette(&mut pixel_color, fractal_rgb, alpha, s_alpha);
            pixel_color.write_to_slice(pixel_slice);
        }
    });
}

/// Helper to find the optimal viewport rotation and zoom for a set of active complex points.
pub fn find_optimal_framing(
    active_points: &[Complex],
    width: u32,
    height: u32,
    default_rotation: Complex,
) -> (f64, Complex) {
    if active_points.is_empty() {
        return (f64::MAX, default_rotation);
    }

    let w_f = width as f64;
    let h_f = height as f64;
    let min_dim = w_f.min(h_f);

    let mut best_zoom = f64::MAX;
    let mut best_rotation = default_rotation;

    for phasor in get_rotation_phasors(ROTATION_STEPS) {
        let inverse_phasor = phasor.conj();
        let mut max_cx_abs = 0.0_f64;
        let mut max_cy_abs = 0.0_f64;

        for &point in active_points {
            let rotated = point * inverse_phasor;
            max_cx_abs = max_cx_abs.max(rotated.re.abs());
            max_cy_abs = max_cy_abs.max(rotated.im.abs());
        }

        let required_zoom =
            (2.0 * max_cx_abs * min_dim / w_f).max(2.0 * max_cy_abs * min_dim / h_f);

        if required_zoom < best_zoom {
            best_zoom = required_zoom;
            best_rotation = phasor;
        }
    }
    (best_zoom, best_rotation)
}

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

    #[test]
    fn test_procedural_effect_resolution() {
        let effect_rand = ProceduralEffect::Random;
        let resolved_rand = effect_rand.resolve();
        assert_ne!(resolved_rand, ProceduralEffect::Random);

        let effect_fractal = ProceduralEffect::Fractal;
        let resolved_fractal = effect_fractal.resolve();
        assert_ne!(resolved_fractal, ProceduralEffect::Fractal);

        let resolved_julia = ProceduralEffect::JuliaSet.resolve();
        assert_eq!(resolved_julia, ProceduralEffect::JuliaSet);
    }

    #[test]
    fn test_smooth_potential_clamping() {
        let z = Complex::new(5.0, 5.0);
        let t = calculate_smooth_potential(50, 100, z);
        assert!((0.0..=1.0).contains(&t));
    }

    #[test]
    fn test_viewport_complex_operations() {
        let specs = ViewportSpecs {
            center: Complex::new(0.5, -0.5),
            zoom: 2.0,
            rotation: Complex::one(),
            is_julia: false,
        };
        let viewport = Viewport::new(100.0, 100.0, &specs);
        let mapped = viewport.map(50.0, 50.0);
        assert!((mapped.re - 0.5).abs() < 1e-9);
    }

    #[test]
    fn test_complex_phasors() {
        let phasors: Vec<Complex> = get_rotation_phasors(ROTATION_STEPS).collect();
        assert_eq!(phasors.len(), ROTATION_STEPS);
        for phasor in phasors {
            let magnitude = phasor.abs();
            assert!((magnitude - 1.0).abs() < 1e-9);
        }
    }

    #[test]
    fn test_optimize_relaxed_viewport() {
        let config = RelaxedViewportConfig {
            width: 100,
            height: 100,
            search_limit: 1.5,
            steps: 10,
            zoom_range: [1.0, 3.0],
            rand_range: [0.9, 1.1],
            fallback_range: [1.2, 2.0],
        };
        let (zoom, rotation) = optimize_relaxed_viewport(config, Complex::one(), |z| z.abs() < 1.0);
        assert!((1.0..=3.0).contains(&zoom));
        assert!((rotation.abs() - 1.0).abs() < 1e-9);
    }
}