use image::{GrayImage, Luma};
use rayon::prelude::*;

const TAU: f32 = 6.283_185_5;
const EROSION_LACUNARITY: f32 = 2.0;
const EROSION_GAIN: f32 = 0.5;
const EROSION_CELL_SCALE: f32 = 0.7;
const EROSION_NORMALIZATION: f32 = 0.5;
const EROSION_ROUNDING: [f32; 4] = [0.1, 0.0, 0.1, 2.0];
const EROSION_ONSET: [f32; 4] = [1.25, 1.25, 2.8, 1.5];
const EROSION_ASSUMED_SLOPE: [f32; 2] = [0.7, 1.0];

#[derive(Clone, Debug)]
pub struct GenerationParams {
    pub image_width: usize,
    pub image_height: usize,
    pub generation_scale: f32,
    pub erosion_scale: f32,
    pub random_seed: u64,
    pub erosion_strength: f32,
    pub erosion_gully_weight: f32,
    pub erosion_detail: f32,
    pub erosion_octaves: usize,
}

impl Default for GenerationParams {
    fn default() -> Self {
        Self {
            image_width: 512,
            image_height: 512,
            generation_scale: 1.0,
            erosion_scale: 0.15,
            random_seed: 0,
            erosion_strength: 0.05,
            erosion_gully_weight: 0.5,
            erosion_detail: 1.5,
            erosion_octaves: 5,
        }
    }
}

#[derive(Clone)]
pub struct GeneratedImages {
    pub width: usize,
    pub height: usize,
    pub heightmap_image: GrayImage,
    pub bumpmap_image: GrayImage,
}

#[derive(Copy, Clone, Debug, Default)]
struct Vec2 {
    x: f32,
    y: f32,
}

impl Vec2 {
    fn new(x: f32, y: f32) -> Self {
        Self { x, y }
    }

    fn dot(self, other: Self) -> f32 {
        self.x * other.x + self.y * other.y
    }

    fn length(self) -> f32 {
        self.dot(self).sqrt()
    }

    fn normalize_safe(self) -> Self {
        let length = self.length();
        if length > 1e-10 {
            Self::new(self.x / length, self.y / length)
        } else {
            self
        }
    }
}

impl std::ops::Add for Vec2 {
    type Output = Self;

    fn add(self, rhs: Self) -> Self::Output {
        Self::new(self.x + rhs.x, self.y + rhs.y)
    }
}

impl std::ops::AddAssign for Vec2 {
    fn add_assign(&mut self, rhs: Self) {
        self.x += rhs.x;
        self.y += rhs.y;
    }
}

impl std::ops::Sub for Vec2 {
    type Output = Self;

    fn sub(self, rhs: Self) -> Self::Output {
        Self::new(self.x - rhs.x, self.y - rhs.y)
    }
}

impl std::ops::Mul<f32> for Vec2 {
    type Output = Self;

    fn mul(self, rhs: f32) -> Self::Output {
        Self::new(self.x * rhs, self.y * rhs)
    }
}

fn mix(a: f32, b: f32, t: f32) -> f32 {
    a + (b - a) * t
}

fn mix2(a: Vec2, b: Vec2, t: f32) -> Vec2 {
    Vec2::new(mix(a.x, b.x, t), mix(a.y, b.y, t))
}

fn clamp(x: f32, lo: f32, hi: f32) -> f32 {
    x.clamp(lo, hi)
}

fn clamp01(x: f32) -> f32 {
    clamp(x, 0.0, 1.0)
}

fn hash22(p: Vec2, random_seed: u64) -> Vec2 {
    let scrambled_x = splitmix64(random_seed ^ 0x9E37_79B9_7F4A_7C15);
    let scrambled_y = splitmix64(random_seed ^ 0xBF58_476D_1CE4_E5B9);
    let seed_x = u64_to_unit_f64(scrambled_x) * 1024.0;
    let seed_y = u64_to_unit_f64(scrambled_y) * 1024.0;
    let kx = 0.318_309_9_f64;
    let ky = 0.367_879_4_f64;

    let px = ((p.x as f64) + seed_x) * kx + ky;
    let py = ((p.y as f64) + seed_y) * ky + kx;
    let val = px * py * (px + py);
    let f_val = val - val.floor();

    let rx = -1.0 + 2.0 * ((16.0 * kx * f_val) - (16.0 * kx * f_val).floor());
    let ry = -1.0 + 2.0 * ((16.0 * ky * f_val) - (16.0 * ky * f_val).floor());
    Vec2::new(rx as f32, ry as f32)
}

fn splitmix64(mut x: u64) -> u64 {
    x = x.wrapping_add(0x9E37_79B9_7F4A_7C15);
    x = (x ^ (x >> 30)).wrapping_mul(0xBF58_476D_1CE4_E5B9);
    x = (x ^ (x >> 27)).wrapping_mul(0x94D0_49BB_1331_11EB);
    x ^ (x >> 31)
}

fn u64_to_unit_f64(value: u64) -> f64 {
    const SCALE: f64 = 1.0 / (u64::MAX as f64);
    (value as f64) * SCALE
}

fn noised(p: Vec2, random_seed: u64) -> (f32, Vec2) {
    let i = Vec2::new(p.x.floor(), p.y.floor());
    let f = Vec2::new(p.x - i.x, p.y - i.y);

    let u = Vec2::new(
        f.x * f.x * f.x * (f.x * (f.x * 6.0 - 15.0) + 10.0),
        f.y * f.y * f.y * (f.y * (f.y * 6.0 - 15.0) + 10.0),
    );
    let du = Vec2::new(
        30.0 * f.x * f.x * (f.x * (f.x - 2.0) + 1.0),
        30.0 * f.y * f.y * (f.y * (f.y - 2.0) + 1.0),
    );

    let ga = hash22(i + Vec2::new(0.0, 0.0), random_seed);
    let gb = hash22(i + Vec2::new(1.0, 0.0), random_seed);
    let gc = hash22(i + Vec2::new(0.0, 1.0), random_seed);
    let gd = hash22(i + Vec2::new(1.0, 1.0), random_seed);

    let va = ga.dot(f - Vec2::new(0.0, 0.0));
    let vb = gb.dot(f - Vec2::new(1.0, 0.0));
    let vc = gc.dot(f - Vec2::new(0.0, 1.0));
    let vd = gd.dot(f - Vec2::new(1.0, 1.0));

    let val = va + u.x * (vb - va) + u.y * (vc - va) + u.x * u.y * (va - vb - vc + vd);

    let deriv_x = ga.x
        + u.x * (gb.x - ga.x)
        + u.y * (gc.x - ga.x)
        + u.x * u.y * (ga.x - gb.x - gc.x + gd.x)
        + du.x * (u.y * (va - vb - vc + vd) + (vb - va));
    let deriv_y = ga.y
        + u.x * (gb.y - ga.y)
        + u.y * (gc.y - ga.y)
        + u.x * u.y * (ga.y - gb.y - gc.y + gd.y)
        + du.y * (u.x * (va - vb - vc + vd) + (vc - va));

    (val, Vec2::new(deriv_x, deriv_y))
}

fn fractal_noise(p: Vec2, freq: f32, octaves: usize, lacunarity: f32, gain: f32, random_seed: u64) -> (f32, Vec2) {
    let mut val = 0.0;
    let mut deriv = Vec2::default();
    let mut nf = freq;
    let mut na = 1.0;
    for _ in 0..octaves {
        let (v, d) = noised(p * nf, random_seed);
        val += v * na;
        deriv += d * (na * nf);
        na *= gain;
        nf *= lacunarity;
    }
    (val, deriv)
}

fn phacelle_noise(p: Vec2, norm_dir: Vec2, freq: f32, offset: f32, normalization: f32, random_seed: u64) -> (Vec2, Vec2) {
    let side_dir = Vec2::new(-norm_dir.y, norm_dir.x) * (freq * TAU);
    let offset_tau = offset * TAU;

    let p_int = Vec2::new(p.x.floor(), p.y.floor());
    let p_frac = p - p_int;

    let mut phase_dir = Vec2::default();
    let mut weight_sum = 0.0;

    for i in -1..3 {
        for j in -1..3 {
            let grid_offset = Vec2::new(i as f32, j as f32);
            let grid_point = p_int + grid_offset;
            let random_offset = hash22(grid_point, random_seed) * 0.5;
            let vector_from_cell_point = p_frac - grid_offset - random_offset;
            let sqr_dist = vector_from_cell_point.dot(vector_from_cell_point);

            let mut weight = (-sqr_dist * 2.0).exp();
            weight = weight.max(0.0) - 0.01111;
            weight = weight.max(0.0);
            weight_sum += weight;

            let wave_input = vector_from_cell_point.dot(side_dir) + offset_tau;
            phase_dir.x += wave_input.cos() * weight;
            phase_dir.y += wave_input.sin() * weight;
        }
    }

    if weight_sum < 1e-6 {
        return (Vec2::default(), side_dir);
    }

    let interpolated = phase_dir * (1.0 / weight_sum);
    let magnitude = interpolated.dot(interpolated).sqrt().max(1.0 - normalization);
    (interpolated * (1.0 / magnitude), side_dir)
}

fn pow_inv(t: f32, power: f32) -> f32 {
    1.0 - (1.0 - clamp01(t)).powf(power)
}

fn ease_out(t: f32) -> f32 {
    let v = 1.0 - clamp01(t);
    1.0 - v * v
}

fn smooth_start(t: f32, smoothing: f32) -> f32 {
    if t >= smoothing {
        t - 0.5 * smoothing
    } else {
        0.5 * t * t / smoothing
    }
}

pub fn generate(params: &GenerationParams) -> GeneratedImages {
    let width = params.image_width;
    let height = params.image_height;
    let scale = params.erosion_scale;
    let strength_base = params.erosion_strength * scale;

    let mut heightmap = vec![0.0_f32; width * height];
    heightmap
        .par_chunks_mut(width)
        .enumerate()
        .for_each(|(y, row)| {
            for (x, value) in row.iter_mut().enumerate() {
                let p = Vec2::new(
                    x as f32 / width as f32,
                    y as f32 / height as f32,
                ) * params.generation_scale;

                let height_freq = 3.0;
                let height_amp = 0.125;
                let (val, deriv) = fractal_noise(p, height_freq, 3, 2.0, 0.1, params.random_seed);

                let n_height = val * height_amp + 0.5;
                let curr_slope = deriv * height_amp;

                let mut fade_target = clamp((n_height - 0.5) / (height_amp * 0.6), -1.0, 1.0);
                let mut freq = 1.0 / (scale * EROSION_CELL_SCALE);
                let slope_length = curr_slope.length().max(1e-10);

                let rounding_for_input = mix(
                    EROSION_ROUNDING[1],
                    EROSION_ROUNDING[0],
                    clamp01(fade_target + 0.5),
                ) * EROSION_ROUNDING[2];
                let mut combi_mask = ease_out(smooth_start(
                    slope_length * EROSION_ONSET[0],
                    rounding_for_input * EROSION_ONSET[0],
                ));

                let mut gully_slope = mix2(
                    curr_slope,
                    curr_slope * (EROSION_ASSUMED_SLOPE[0] / slope_length),
                    EROSION_ASSUMED_SLOPE[1],
                );

                let mut strength = strength_base;
                let mut rounding_mult = 1.0;
                let mut total_h_delta = 0.0;
                let mut total_strength = 0.0;

                for _ in 0..params.erosion_octaves {
                    let (phacelle_vec, side_dir) = phacelle_noise(
                        p * freq,
                        gully_slope.normalize_safe(),
                        EROSION_CELL_SCALE,
                        0.25,
                        EROSION_NORMALIZATION,
                        params.random_seed,
                    );

                    let p_deriv_dir = side_dir * -freq;
                    let sloping = phacelle_vec.y.abs();
                    gully_slope += p_deriv_dir
                        * phacelle_vec.y.signum()
                        * strength
                        * params.erosion_gully_weight;

                    let gullies_h = phacelle_vec.x;
                    let faded_gullies_h = mix(
                        fade_target,
                        gullies_h * params.erosion_gully_weight,
                        combi_mask,
                    );
                    total_h_delta += faded_gullies_h * strength;
                    total_strength += strength;
                    fade_target = faded_gullies_h;

                    let rounding_for_octave = mix(
                        EROSION_ROUNDING[1],
                        EROSION_ROUNDING[0],
                        clamp01(phacelle_vec.x + 0.5),
                    ) * rounding_mult;
                    let new_mask = ease_out(smooth_start(
                        sloping * EROSION_ONSET[1],
                        rounding_for_octave * EROSION_ONSET[1],
                    ));
                    combi_mask = pow_inv(combi_mask, params.erosion_detail) * new_mask;
                    strength *= EROSION_GAIN;
                    freq *= EROSION_LACUNARITY;
                    rounding_mult *= EROSION_ROUNDING[3];
                }

                let terrain_height_offset = [-0.65, 0.0];
                let final_fade_target = clamp((n_height - 0.5) / (height_amp * 0.6), -1.0, 1.0);
                let offset = mix(
                    terrain_height_offset[0],
                    -final_fade_target,
                    terrain_height_offset[1],
                ) * total_strength;

                *value = n_height + total_h_delta + offset;
            }
        });

    let heightmap_normalized = normalize_values(&heightmap);
    let bumpmap = apply_shading(&heightmap, width, height, 0.5);
    let heightmap_image = grayscale_image_from_values(&heightmap_normalized, width, height);
    let bumpmap_image = grayscale_image_from_values(&bumpmap, width, height);

    GeneratedImages {
        width,
        height,
        heightmap_image,
        bumpmap_image,
    }
}

fn normalize_values(values: &[f32]) -> Vec<f32> {
    let mut min_val = f32::INFINITY;
    let mut max_val = f32::NEG_INFINITY;
    for value in values {
        min_val = min_val.min(*value);
        max_val = max_val.max(*value);
    }

    if (max_val - min_val).abs() < 1e-10 {
        return vec![0.0; values.len()];
    }

    values
        .iter()
        .map(|value| (value - min_val) / (max_val - min_val))
        .collect()
}

fn grayscale_image_from_values(values: &[f32], width: usize, height: usize) -> GrayImage {
    let mut image = GrayImage::new(width as u32, height as u32);
    for y in 0..height {
        for x in 0..width {
            let value = values[y * width + x].clamp(0.0, 1.0);
            image.put_pixel(x as u32, y as u32, Luma([(value * 255.0).round() as u8]));
        }
    }
    image
}

fn apply_shading(heightmap: &[f32], width: usize, height: usize, strength: f32) -> Vec<f32> {
    let mut output = vec![0.0; width * height];

    let lx = -1.0_f32;
    let ly = -1.0_f32;
    let lz = 1.0_f32;
    let l_len = (lx * lx + ly * ly + lz * lz).sqrt();
    let (lx, ly, lz) = (lx / l_len, ly / l_len, lz / l_len);

    for y in 0..height {
        for x in 0..width {
            let left = sample_height(heightmap, width, height, x.saturating_sub(1), y);
            let right = sample_height(heightmap, width, height, (x + 1).min(width - 1), y);
            let up = sample_height(heightmap, width, height, x, y.saturating_sub(1));
            let down = sample_height(heightmap, width, height, x, (y + 1).min(height - 1));

            let mut dx = (right - left) * 0.5 * strength;
            let mut dy = (down - up) * 0.5 * strength;
            dx *= width as f32;
            dy *= height as f32;

            let n_len = (dx * dx + dy * dy + 1.0).sqrt();
            let nx = -dx / n_len;
            let ny = -dy / n_len;
            let nz = 1.0 / n_len;
            let dot = nx * lx + ny * ly + nz * lz;
            output[y * width + x] = dot.clamp(0.0, 1.0);
        }
    }

    output
}

fn sample_height(heightmap: &[f32], width: usize, _height: usize, x: usize, y: usize) -> f32 {
    heightmap[y * width + x]
}