///|
/// Higher-level scalar-field tools for deterministic texture and contour work.
pub(all) struct FieldPlusStats {
  minimum : Double
  maximum : Double
  mean : Double
  variance : Double
  energy : Double
  positive_fraction : Double
  negative_fraction : Double
} derive(Eq, Debug, ToJson)

///|
pub(all) struct FieldVector {
  dx : Double
  dy : Double
} derive(Eq, Debug, ToJson)

///|
pub(all) enum FieldCombineMode {
  AddFields
  SubtractFields
  MultiplyFields
  MinimumField
  MaximumField
  DifferenceField
  AverageField
} derive(Eq, Debug, ToJson)

///|
pub(all) struct FieldSampleReport {
  width : Int
  height : Int
  sample_count : Int
  min : Double
  max : Double
  mean : Double
  rms : Double
  gradient_mean : Double
  gradient_max : Double
} derive(Eq, Debug, ToJson)

///|
fn fp_clamp(value : Double, low : Double, high : Double) -> Double {
  value.max(low).min(high)
}

///|
pub fn field_plus_stats(field : ScalarField) -> FieldPlusStats {
  let count = field.values.length()
  if count == 0 {
    {
      minimum: 0.0,
      maximum: 0.0,
      mean: 0.0,
      variance: 0.0,
      energy: 0.0,
      positive_fraction: 0.0,
      negative_fraction: 0.0,
    }
  } else {
    let mut minimum = field.values[0]
    let mut maximum = minimum
    let mut total = 0.0
    let mut energy = 0.0
    let mut positive = 0
    let mut negative = 0
    for value in field.values {
      minimum = minimum.min(value)
      maximum = maximum.max(value)
      total += value
      energy += value * value
      if value >= 0.0 {
        positive += 1
      } else {
        negative += 1
      }
    }
    let mean = total / count.to_double()
    let mut variance = 0.0
    for value in field.values {
      let delta = value - mean
      variance += delta * delta
    }
    {
      minimum,
      maximum,
      mean,
      variance: variance / count.to_double(),
      energy: energy / count.to_double(),
      positive_fraction: positive.to_double() / count.to_double(),
      negative_fraction: negative.to_double() / count.to_double(),
    }
  }
}

///|
pub fn field_plus_get(
  field : ScalarField,
  x : Int,
  y : Int,
  fallback : Double,
) -> Double {
  field.get(x, y).unwrap_or(fallback)
}

///|
pub fn field_plus_gradient(
  field : ScalarField,
  x : Int,
  y : Int,
) -> FieldVector {
  let left = field_plus_get(field, x - 1, y, field_plus_get(field, x, y, 0.0))
  let right = field_plus_get(field, x + 1, y, field_plus_get(field, x, y, 0.0))
  let top = field_plus_get(field, x, y - 1, field_plus_get(field, x, y, 0.0))
  let bottom = field_plus_get(field, x, y + 1, field_plus_get(field, x, y, 0.0))
  { dx: (right - left) / 2.0, dy: (bottom - top) / 2.0 }
}

///|
pub fn field_plus_gradient_magnitude(
  field : ScalarField,
  x : Int,
  y : Int,
) -> Double {
  let gradient = field_plus_gradient(field, x, y)
  (gradient.dx * gradient.dx + gradient.dy * gradient.dy).sqrt()
}

///|
pub fn field_plus_gradient_field(field : ScalarField) -> ScalarField {
  let result = ScalarField::new(field.width, field.height)
  for y in 0.. Double {
  let center = field_plus_get(field, x, y, 0.0)
  field_plus_get(field, x - 1, y, center) +
  field_plus_get(field, x + 1, y, center) +
  field_plus_get(field, x, y - 1, center) +
  field_plus_get(field, x, y + 1, center) -
  4.0 * center
}

///|
pub fn field_plus_laplacian_field(field : ScalarField) -> ScalarField {
  let result = ScalarField::new(field.width, field.height)
  for y in 0.. ScalarField {
  let width = left.width.min(right.width).max(0)
  let height = left.height.min(right.height).max(0)
  let result = ScalarField::new(width, height)
  for y in 0.. a + b
        SubtractFields => a - b
        MultiplyFields => a * b
        MinimumField => a.min(b)
        MaximumField => a.max(b)
        DifferenceField => (a - b).abs()
        AverageField => (a + b) / 2.0
      }
      let _ = result.set(x, y, value)
    }
  }
  result
}

///|
pub fn field_plus_scale(
  field : ScalarField,
  factor : Double,
  bias : Double,
) -> ScalarField {
  field.map(fn(value) { value * factor + bias })
}

///|
pub fn field_plus_remap(
  field : ScalarField,
  low : Double,
  high : Double,
) -> ScalarField {
  let stats = field_plus_stats(field)
  if stats.maximum == stats.minimum {
    field.map(fn(_value) { low })
  } else {
    field.map(fn(value) {
      low +
      (value - stats.minimum) / (stats.maximum - stats.minimum) * (high - low)
    })
  }
}

///|
pub fn field_plus_clamp(
  field : ScalarField,
  low : Double,
  high : Double,
) -> ScalarField {
  field.map(fn(value) { fp_clamp(value, low, high) })
}

///|
pub fn field_plus_threshold(
  field : ScalarField,
  threshold : Double,
  low : Double,
  high : Double,
) -> ScalarField {
  field.map(fn(value) { if value >= threshold { high } else { low } })
}

///|
pub fn field_plus_smoothstep(
  field : ScalarField,
  edge0 : Double,
  edge1 : Double,
) -> ScalarField {
  field.map(fn(value) {
    let t = fp_clamp(
      (value - edge0) / (edge1 - edge0).abs().max(0.000001),
      0.0,
      1.0,
    )
    t * t * (3.0 - 2.0 * t)
  })
}

///|
pub fn field_plus_pow(field : ScalarField, exponent : Double) -> ScalarField {
  field.map(fn(value) {
    let sign = if value < 0.0 { -1.0 } else { 1.0 }
    sign * @math.pow(value.abs(), exponent)
  })
}

///|
pub fn field_plus_abs(field : ScalarField) -> ScalarField {
  field.map(fn(value) { value.abs() })
}

///|
pub fn field_plus_invert(field : ScalarField) -> ScalarField {
  field.map(fn(value) { 1.0 - value })
}

///|
pub fn field_plus_rotate(
  field : ScalarField,
  quarter_turns : Int,
) -> ScalarField {
  let turns = (quarter_turns % 4 + 4) % 4
  if turns == 0 {
    field
  } else {
    let width = if turns % 2 == 0 { field.width } else { field.height }
    let height = if turns % 2 == 0 { field.height } else { field.width }
    let result = ScalarField::new(width, height)
    for y in 0.. (field.height - 1 - y, x)
          2 => (field.width - 1 - x, field.height - 1 - y)
          _ => (y, field.width - 1 - x)
        }
        let _ = result.set(
          destination.0,
          destination.1,
          field_plus_get(field, x, y, 0.0),
        )
      }
    }
    result
  }
}

///|
pub fn field_plus_flip_x(field : ScalarField) -> ScalarField {
  let result = ScalarField::new(field.width, field.height)
  for y in 0.. ScalarField {
  let result = ScalarField::new(field.width, field.height)
  for y in 0.. ScalarField {
  let result = ScalarField::new(width, height)
  if width > 0 && height > 0 && field.width > 0 && field.height > 0 {
    for y in 0.. ScalarField {
  let result = ScalarField::new(width, height)
  if width > 0 && height > 0 && field.width > 0 && field.height > 0 {
    for y in 0.. Array[Int] {
  let result = []
  let count = bins.max(0)
  for _ in 0.. 0 && high > low {
    for value in field.values {
      let index = ((value - low) / (high - low) * count.to_double())
        .to_int()
        .min(count - 1)
        .max(0)
      result[index] += 1
    }
  }
  result
}

///|
pub fn field_plus_quantile(field : ScalarField, quantile : Double) -> Double {
  if field.values.length() == 0 {
    0.0
  } else {
    let values = field.values.copy()
    for i in 1.. 0 && values[position - 1] > value {
        values[position] = values[position - 1]
        position -= 1
      }
      values[position] = value
    }
    values[(fp_clamp(quantile, 0.0, 1.0) * (values.length() - 1).to_double()).to_int()]
  }
}

///|
pub fn field_plus_median(field : ScalarField) -> Double {
  field_plus_quantile(field, 0.5)
}

///|
pub fn field_plus_sample_bilinear(
  field : ScalarField,
  x : Double,
  y : Double,
) -> Double {
  if field.width <= 0 || field.height <= 0 {
    0.0
  } else {
    let gx = fp_clamp(x, 0.0, (field.width - 1).to_double())
    let gy = fp_clamp(y, 0.0, (field.height - 1).to_double())
    let x0 = gx.to_int()
    let y0 = gy.to_int()
    let x1 = (x0 + 1).min(field.width - 1)
    let y1 = (y0 + 1).min(field.height - 1)
    let tx = gx - x0.to_double()
    let ty = gy - y0.to_double()
    let top = field_plus_get(field, x0, y0, 0.0) * (1.0 - tx) +
      field_plus_get(field, x1, y0, 0.0) * tx
    let bottom = field_plus_get(field, x0, y1, 0.0) * (1.0 - tx) +
      field_plus_get(field, x1, y1, 0.0) * tx
    top * (1.0 - ty) + bottom * ty
  }
}

///|
pub fn field_plus_sample_report(field : ScalarField) -> FieldSampleReport {
  let stats = field_plus_stats(field)
  let gradient = field_plus_gradient_field(field)
  let gradient_stats = field_plus_stats(gradient)
  {
    width: field.width,
    height: field.height,
    sample_count: field.values.length(),
    min: stats.minimum,
    max: stats.maximum,
    mean: stats.mean,
    rms: stats.energy.sqrt(),
    gradient_mean: gradient_stats.mean,
    gradient_max: gradient_stats.maximum,
  }
}

///|
pub fn field_plus_report_csv(report : FieldSampleReport) -> String {
  "width,height,sample_count,min,max,mean,rms,gradient_mean,gradient_max\n" +
  report.width.to_string() +
  "," +
  report.height.to_string() +
  "," +
  report.sample_count.to_string() +
  "," +
  report.min.to_string() +
  "," +
  report.max.to_string() +
  "," +
  report.mean.to_string() +
  "," +
  report.rms.to_string() +
  "," +
  report.gradient_mean.to_string() +
  "," +
  report.gradient_max.to_string() +
  "\n"
}

///|
pub fn field_plus_to_raster(
  field : ScalarField,
  low : Color,
  high : Color,
) -> Raster {
  let result = Raster::new(field.width, field.height)
  let normalized = field.normalize()
  for y in 0.. Array[ContourLine] {
  let result = []
  for level in levels {
    for line in isolines(field, level, tolerance=1.0) {
      result.push(line)
    }
  }
  result
}

///|
pub fn[N : NoiseSampler] field_plus_noise_stack(
  width : Int,
  height : Int,
  sampler : N,
  base_scale : Double,
  octaves : Int,
  persistence : Double,
  lacunarity : Double,
) -> ScalarField {
  let result = ScalarField::new(width, height)
  for y in 0.. ScalarField {
  let result = ScalarField::new(field.width, field.height)
  for y in 0.. ScalarField {
  let result = ScalarField::new(field.width, field.height)
  match direction {
    Row =>
      for y in 0..
      for x in 0.. ScalarField {
  let result = ScalarField::new(field.width, field.height, fill=1000000.0)
  let seeds = field.to_points(threshold)
  for y in 0.. ScalarField {
  let result = ScalarField::new(width, height)
  for y in 0..