wallswitch 0.60.8

//! 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 serde::{Deserialize, Serialize};
use std::{f32::consts::LOG2_E, io::Error, path::Path, thread};

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

/// 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: u32 = 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 {
    /// 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.
    ///
    /// Uses non-generic &Path slices to ensure the trait remains dyn compatible.
    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(())
    }
}

/// 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.
    ///
    /// * Characteristics: Generates glowing, wave-like atmospheric filaments with smooth edge transitions
    ///   mimicking planetary polar auroras.
    /// * Creator: Modeled using transcendental multi-frequency wave equations combined with coordinate-space
    ///   decay falloffs.
    /// * Generator function: Evaluated using composite sinusoidal waves mapping local coordinate inputs (x, y)
    ///   against continuous densities:
    ///   alpha = 0.25 * (sin(d_u * x) + cos(d_v * y) + sin(d_w * x + rho) + cos(sqrt(u^2 + v^2) * d_w4))
    ///   coupled with radial quadratic dampening to fade margins near boundaries:
    ///   Intensity = sin(alpha * pi)^2 * (1.0 - 0.45 * sqrt(delta_x^2 + delta_y^2))
    #[value(name = "aurora")]
    CosmicAurora,

    /// Procedural Starfield / Bokeh generator.
    ///
    /// * Characteristics: Projects soft, circular glowing stars and light orbs of varying intensities
    ///   and high-contrast neon colors.
    /// * Creator: Constructed via random coordinate sampling mapped against continuous, non-linear
    ///   light field structures.
    /// * Generator function: Calculated using spatial coordinates evaluated against an exponential Gaussian
    ///   falloff curve centered on randomized center points (x_s, y_s):
    ///   I(d) = I_0 * exp(-d^2 / (2 * sigma^2))
    ///   where d = sqrt((x - x_s)^2 + (y - y_s)^2) represents the Euclidean distance from the active pixel
    ///   to the stellar epicenter, sigma represents the star radius, and I_0 is the peak scalar intensity.
    #[value(name = "star")]
    Starfield,

    /// Fractal mode selector.
    ///
    /// * Characteristics: Randomly selects between the Julia Set, Newton-Raphson Basins, or Nova Julia
    ///   fractal overlays for the active generation cycle.
    #[value(name = "fractal")]
    Fractal,

    /// Fully randomized mode selector.
    ///
    /// * Characteristics: Dynamically and independently selects one of the active procedural overlays
    ///   for each physical monitor.
    #[value(name = "random")]
    Random,
}

impl ProceduralEffect {
    /// Returns the user-friendly display name of the active effect.
    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",
        }
    }

    /// Translates a logical placeholder effect (such as `Random` or `Fractal`) into a concrete, renderable variant.
    pub fn resolve(self) -> Self {
        match self {
            Self::Random => match get_random_integer(0, 5) {
                // Fractal
                0 => Self::JuliaSet,
                1 => Self::Mandelbrot,
                2 => Self::NewtonBasins,
                3 => Self::NovaJulia,
                // Others
                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,
        }
    }

    /// Dynamic factory method that constructs and returns the resolved image renderer.
    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
        )
    }
}

/// 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)
}

/// Applies power-law (Gamma) stretching to enhance the visual contrast of fractal filaments.
#[inline(always)]
pub fn stretch_potential(raw_t: f32) -> f32 {
    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) -> f32 {
    if i < max_iterations {
        let mag2 = z.norm_sq();
        let smooth_i = if mag2 > 4.0 {
            let log_zn = (mag2.ln() * 0.5) as f32; // ln(|z_n|)
            let nu = log_zn.ln() * LOG2_E;
            (i as f32 + 1.0 - nu).max(0.0)
        } else {
            i as f32
        };

        // Enforce a minimum escape iteration threshold to keep the outer space 100% transparent
        let min_render_iter = 32.0_f32;
        if smooth_i < min_render_iter {
            return 0.0;
        }

        // Map the active rendering interval linearly into [0.0, 1.0] before stretching
        let normalized = (smooth_i - min_render_iter) / (max_iterations as f32 - min_render_iter);
        stretch_potential(normalized)
    } else {
        1.0 // Set interior remains fully transparent
    }
}

/// 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.norm();
        let dz_mag = dz.norm();
        if z_mag > 0.0 && dz_mag > 0.0 {
            // Standard DEM formula: d = 2 * |z| * ln(|z|) / |dz|
            return 2.0 * z_mag * z_mag.ln() / dz_mag;
        }
    }
    0.0
}

/// Linearizes sRGB values to perform mathematically correct alpha-blending,
/// preventing dark-boundary artifacts.
#[inline]
pub fn blend_channels_gamma(bg: u8, fg: f32, alpha: f32) -> u8 {
    let bg_f = bg as f32 / 255.0;
    let bg_linear = bg_f * bg_f; // Fast gamma = 2.0 approximation

    let fg_f = fg / 255.0;
    let fg_linear = fg_f * fg_f;

    let blended_linear = bg_linear * (1.0 - alpha) + fg_linear * alpha;

    (blended_linear.sqrt() * 255.0).clamp(0.0, 255.0) as u8 // Re-encode to sRGB
}

/// Standard smoothstep mathematical interpolation function.
#[inline]
pub fn smoothstep(edge0: f32, edge1: f32, x: f32) -> f32 {
    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]
pub fn blend_and_vignette(
    pixel: &mut ColorRGB,
    fractal_rgb: ColorRGB,
    alpha: f32,
    shadow_alpha: f32,
) {
    // 1. Apply drop shadow (ambient occlusion) using ColorRGB vector scaling
    if shadow_alpha > 0.005 {
        *pixel = pixel.scale(1.0 - shadow_alpha);
    }

    // 2. Perform smooth gamma-corrected linear interpolation (lerp)
    if alpha > 0.005 {
        let bg_linear = pixel.squared();
        let fg_linear = fractal_rgb.squared();

        // High-performance component-wise blending: bg_linear * (1 - alpha) + fg_linear * alpha
        let blended_linear = fg_linear.lerp(&bg_linear, alpha);

        *pixel = blended_linear.sqrt();
    }
}

/// Executes complex escape-time steps for both Julia and Mandelbrot fractal types.
///
/// Unifies the complex escape mathematics, returning the escape count,
/// final coordinate `z`, and boundary derivative `dz`.
#[inline(always)]
pub fn compute_escape_iterations(
    fractal_type: ProceduralEffect,
    init: Complex,
    c: Complex,
    scan_iterations: u32,
) -> (u32, Complex, Complex) {
    let (mut z, param) = if fractal_type == ProceduralEffect::JuliaSet {
        (init, c)
    } else {
        // Cardioid optimization boundary check
        let q = (init.re - 0.25) * (init.re - 0.25) + init.im * init.im;
        if q * (q + (init.re - 0.25)) < 0.25 * init.im * init.im {
            return (
                scan_iterations,
                Complex::new(0.0, 0.0),
                Complex::new(0.0, 0.0),
            );
        }
        // Period-2 bulb check
        if (init.re + 1.0) * (init.re + 1.0) + init.im * init.im < 0.0625 {
            return (
                scan_iterations,
                Complex::new(0.0, 0.0),
                Complex::new(0.0, 0.0),
            );
        }
        (Complex::new(0.0, 0.0), init)
    };

    let mut dz = if fractal_type == ProceduralEffect::JuliaSet {
        Complex::new(1.0, 0.0)
    } else {
        Complex::new(0.0, 0.0)
    };

    let mut i = 0;
    while i < scan_iterations {
        if z.norm_sq() > 4.0 {
            break;
        }

        let add_factor = if fractal_type == ProceduralEffect::JuliaSet {
            Complex::new(0.0, 0.0)
        } else {
            Complex::new(1.0, 0.0)
        };
        // Analytical boundary tracking: dz = 2 * z * dz + step_add
        dz = 2.0 * z * dz + add_factor;
        z = z * z + param;
        i += 1;
    }
    (i, z, dz)
}

/// Computes a flat-center circular edge-fade window (vignette/containment)
/// with a smooth transition near the boundary radius.
#[inline]
pub fn calculate_circular_fade(z: Complex, max_radius: f32, flat_ratio: f32) -> f32 {
    let dist = z.norm() as f32;
    let r = dist / max_radius;

    if r < flat_ratio {
        1.0
    } else if r < 1.0 {
        // Smooth transition from 1.0 down to 0.0 using smoothstep
        let t = (r - flat_ratio) / (1.0 - flat_ratio);
        1.0 - t * t * (3.0 - 2.0 * t)
    } else {
        0.0
    }
}

/// Returns an iterator over unit complex phasors representing structural rotation angles.
#[inline]
pub fn get_rotation_phasors() -> impl Iterator<Item = Complex> {
    (0..ROTATION_STEPS).map(|step| {
        let angle_deg = (step * 360 / ROTATION_STEPS) as f64;
        let rad = angle_deg.to_radians();
        Complex::new(rad.cos(), rad.sin())
    })
}

/// Viewport configuration parameters representing current camera scale and rotation orientation.
pub struct ViewportSpecs {
    /// Focal complex center point.
    pub center: Complex,
    /// Zoom translation scaling index.
    pub zoom: f64,
    /// Cosine of rotation angle.
    pub cos_angle: f64,
    /// Sine of rotation angle.
    pub sin_angle: f64,
    /// Toggle determining whether mapping centers relative to Julia origins.
    pub is_julia: bool,
}

/// Viewport mapper that transforms physical coordinate grids into complex space.
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 {
    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;

        // dx is the horizontal scaling phasor: scale * e^(i * theta)
        let dx = scale * Complex::new(specs.cos_angle, specs.sin_angle);

        // dy is rotated counter-clockwise by exactly 90 degrees (multiplied by imaginary unit i)
        let dy = dx * Complex::new(0.0, 1.0);

        let v_center = if specs.is_julia {
            Complex::new(0.0, 0.0)
        } else {
            specs.center
        };

        // start = v_center - (cx_off * dx) - (cy_off * dy)
        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
    }
}

/// Escape results from relaxed Newton or Nova Julia iterations.
#[derive(Debug, Clone, Copy, PartialEq)]
pub struct RelaxedEscape {
    /// Number of iterations executed before escape or convergence.
    pub iterations: u32,
    /// Maximum iteration limit configured for the sweep.
    pub max_iterations: u32,
    /// Final successive term distance norm (|z_{n+1} - z_n|^2).
    pub diff_norm: f64,
    /// Final coordinate value reached.
    pub z_final: Complex,
}

impl RelaxedEscape {
    /// Creates a new EscapeState record representing loop execution metrics.
    pub fn new(iterations: u32, max_iterations: u32, diff_norm: f64, z_final: Complex) -> Self {
        Self {
            iterations,
            max_iterations,
            diff_norm,
            z_final,
        }
    }
}

/// 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, f32, f32) {
    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) as f32;

    let thickness = 2.5_f32;
    let max_radius = thickness * 2.0;
    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;

    let core = if dist_pixels < 1.2 {
        1.0 - (dist_pixels / 1.2)
    } else {
        0.0
    };

    let ripple_freq = 12.0_f32;
    let ripple_wave = (t * std::f32::consts::PI * ripple_freq).sin().abs();
    let nested_detail = (1.0 - smoothstep(0.0, 0.4, 1.0 - ripple_wave)) * (1.0 - norm_dist);

    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;
    let shadow_profile = (1.0 - norm_shadow * norm_shadow).powi(2) * 0.35;

    let angle = z.arg() as f32;
    let light = 0.65_f32 + 0.35_f32 * (angle * 4.0).cos().abs();

    let t_cycled = (t * 2.0) % 1.0;

    // Smooth dual-tone gradient interpolation using NeonColor's new lerp API
    let secondary = color_palette.rotated();
    let core_color = if t_cycled < 0.5 {
        let factor = t_cycled * 2.0;
        color_palette.color_rgb.lerp(&secondary, 1.0 - factor)
    } else {
        let factor = (t_cycled - 0.5) * 2.0;
        secondary.lerp(&color_palette.color_rgb, 1.0 - factor)
    };

    // 1. Calculate complementary border and normalize it to maximize vibrancy
    let border_color = core_color.complementary().saturate_components();

    // 2. Blend core and border colors linearly
    let color_blend = norm_dist.powi(2);
    let blended = border_color.lerp(&core_color, color_blend);

    // 3. Apply local illumination shading and global brightness boosts
    let brightness_boost = 1.20_f32;
    let rgb = blended.scale(light * brightness_boost).clamp_bounds();

    let iteration_fade = if i < 16 { (i as f32 - 3.0) / 13.0 } else { 1.0 };

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

/// Unifies relaxed Newton-Raphson/Nova Julia rendering math to ensure highly consistent visualization.
#[inline(always)]
pub fn color_relaxed_newton_fractal(
    escape: &RelaxedEscape,
    color_palette: NeonColor,
    edge_fade: f32,
    ln_epsilon: f32,
    is_nova: bool,
) -> (ColorRGB, f32, f32) {
    if escape.iterations >= escape.max_iterations {
        return (ColorRGB::default(), 0.0, 0.0);
    }

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

    let ripple_frequency = 0.50_f32;
    let raw_wave = (smooth_i * ripple_frequency * std::f32::consts::PI)
        .sin()
        .abs();

    let norm_dist = if is_nova {
        raw_wave.powf(2.5)
    } else {
        raw_wave
    };

    let core = if is_nova {
        if norm_dist > 0.92 {
            (norm_dist - 0.92) / 0.08
        } else {
            0.0
        }
    } else {
        if norm_dist > 0.95 {
            (norm_dist - 0.95) / 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
    };

    let angle = escape.z_final.arg() as f32;
    let light = if is_nova {
        0.75_f32 + 0.25_f32 * (angle * 4.0).cos().abs()
    } else {
        0.70_f32 + 0.30_f32 * (angle * 3.0).cos().abs()
    };

    let t_cycled = (smooth_i * 0.08) % 1.0;

    // Smooth dual-tone gradient interpolation using NeonColor's new lerp API
    let secondary = color_palette.rotated();
    let core_color = if is_nova {
        let t_cos = (t_cycled * std::f32::consts::PI).cos() * 0.5 + 0.5;
        color_palette.color_rgb.lerp(&secondary, t_cos)
    } else {
        color_palette.color_rgb.lerp(&secondary, 1.0 - t_cycled)
    };

    // 1. Calculate complementary border and normalize it to maximize vibrancy
    let border_color = core_color.complementary().saturate_components();

    // 2. Blend core and border colors linearly
    let blended = if is_nova {
        let color_blend = norm_dist.powf(3.0);
        border_color.lerp(&core_color, color_blend)
    } else {
        border_color.lerp(&core_color, norm_dist)
    };

    // 3. Apply local illumination shading and global brightness boosts
    let brightness_boost = if is_nova { 1.45_f32 } else { 1.25_f32 };
    let rgb = blended.scale(light * brightness_boost).clamp_bounds();

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

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

/// 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,
    cos_angle: f64,
    sin_angle: f64,
    center: Complex,
    is_julia: bool,
    pixel_fn: F,
) where
    F: Fn(Complex, f64) -> (ColorRGB, f32, f32) + Send + Sync,
{
    let (width, height) = rgb_img.dimensions();
    let w_f = width as f64;
    let h_f = height as f64;

    let specs = ViewportSpecs {
        center,
        zoom,
        cos_angle,
        sin_angle,
        is_julia,
    };
    let viewport = Viewport::new(w_f, h_f, &specs);

    let (mut rows, _) = partition_rows(rgb_img);

    let cores = thread::available_parallelism()
        .map(|n| n.get())
        .unwrap_or(4);
    let chunk_size = (rows.len() / cores).max(1);

    let min_dim = w_f.min(h_f);
    let scale = zoom / min_dim;

    thread::scope(|scope| {
        let viewport_ref = &viewport;
        let pixel_fn_ref = &pixel_fn;
        for chunk in rows.chunks_mut(chunk_size) {
            scope.spawn(move || {
                for (y, row_data) in chunk.iter_mut() {
                    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_ref.map(x_f, y_f);

                        let (fractal_rgb, alpha, s_alpha) = pixel_fn_ref(z_init, scale);

                        // Load background pixel, mapping [0..255] u8 range to [0.0..1.0] linear scale
                        let mut pixel_color = ColorRGB::from_slice(pixel_slice);

                        // Execute procedural rendering and mixing operations directly on the ColorRGB representation
                        blend_and_vignette(&mut pixel_color, fractal_rgb, alpha, s_alpha);

                        // Persist processed color channels back to the image buffer
                        pixel_color.write_to_slice(pixel_slice);
                    }
                }
            });
        }
    });
}

/// Helper to find the optimal viewport rotation and zoom for a set of active complex points.
///
/// Returns `(best_zoom, best_cos, best_sin)` representing the optimal framing parameters.
pub fn find_optimal_framing(
    active_points: &[Complex],
    width: u32,
    height: u32,
    default_cos: f64,
    default_sin: f64,
) -> (f64, f64, f64) {
    if active_points.is_empty() {
        return (f64::MAX, default_cos, default_sin);
    }

    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_cos = default_cos;
    let mut best_sin = default_sin;

    for phasor in get_rotation_phasors() {
        let mut max_cx_abs = 0.0_f64;
        let mut max_cy_abs = 0.0_f64;

        let inverse_phasor = phasor.conj();

        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 zoom_x = 2.0 * max_cx_abs * min_dim / w_f;
        let zoom_y = 2.0 * max_cy_abs * min_dim / h_f;
        let required_zoom = zoom_x.max(zoom_y);

        if required_zoom < best_zoom {
            best_zoom = required_zoom;
            best_cos = phasor.re;
            best_sin = phasor.im;
        }
    }

    (best_zoom, best_cos, best_sin)
}

#[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,
            cos_angle: 1.0,
            sin_angle: 0.0,
            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().collect();
        assert_eq!(phasors.len() as u32, ROTATION_STEPS);
        for phasor in phasors {
            let magnitude = phasor.norm();
            assert!((magnitude - 1.0).abs() < 1e-9);
        }
    }

    #[test]
    fn test_render_fractal_parallel() {
        let mut img = RgbImage::new(10, 10);
        let center = Complex::new(0.0, 0.0);
        render_fractal_parallel(&mut img, 2.0, 1.0, 0.0, center, true, |_z_init, _scale| {
            (ColorRGB::new(1.0, 0.0, 0.0), 1.0, 0.0)
        });
        let pixel = img.get_pixel(5, 5);
        assert!(pixel[0] > 0);
    }

    #[test]
    fn test_find_optimal_framing_empty() {
        let (zoom, cos, sin) = find_optimal_framing(&[], 100, 100, 1.0, 0.0);
        assert_eq!(zoom, f64::MAX);
        assert_eq!(cos, 1.0);
        assert_eq!(sin, 0.0);
    }

    #[test]
    fn test_find_optimal_framing_with_points() {
        let points = vec![Complex::new(1.0, 1.0), Complex::new(-1.0, -1.0)];
        let (zoom, cos, sin) = find_optimal_framing(&points, 100, 100, 1.0, 0.0);
        assert!(zoom > 0.0);
        assert!(zoom < f64::MAX);
        assert!((cos * cos + sin * sin - 1.0).abs() < 1e-9);
    }
}