///|
/// Return the dimensions of a scalar field.
pub fn Field2D::dimensions(self : Field2D) -> Size {
  self.size
}

///|
/// Sum a half-open integer region after clipping to the field.
pub fn Field2D::region_sum(
  self : Field2D,
  x0~ : Int,
  y0~ : Int,
  x1~ : Int,
  y1~ : Int,
) -> Double {
  let left = x0.clamp(min=0, max=self.size.width)
  let top = y0.clamp(min=0, max=self.size.height)
  let right = x1.clamp(min=left, max=self.size.width)
  let bottom = y1.clamp(min=top, max=self.size.height)
  let mut total = 0.0
  for y in top.. Double {
  let left = x0.clamp(min=0, max=self.size.width)
  let top = y0.clamp(min=0, max=self.size.height)
  let right = x1.clamp(min=left, max=self.size.width)
  let bottom = y1.clamp(min=top, max=self.size.height)
  let count = (right - left) * (bottom - top)
  safe_divide(
    self.region_sum(x0=left, y0=top, x1=right, y1=bottom),
    count.to_double(),
    fallback=0.0,
  )
}

///|
/// Count values in a half-open region that pass a predicate.
pub fn Field2D::region_count(
  self : Field2D,
  x0~ : Int,
  y0~ : Int,
  x1~ : Int,
  y1~ : Int,
  level~ : Double,
) -> Int {
  let left = x0.clamp(min=0, max=self.size.width)
  let top = y0.clamp(min=0, max=self.size.height)
  let right = x1.clamp(min=left, max=self.size.width)
  let bottom = y1.clamp(min=top, max=self.size.height)
  let mut count = 0
  for y in top..= level {
        count += 1
      }
    }
  }
  count
}

///|
/// Extract a clipped scalar subfield.
pub fn Field2D::crop(
  self : Field2D,
  x0~ : Int,
  y0~ : Int,
  width~ : Int,
  height~ : Int,
) -> Field2D {
  let result = Field2D::new(
    size=Size::new(width=width.max(0), height=height.max(0)),
  )
  for y in 0.. Field2D {
  let result = Field2D::new(
    size=Size::new(width=self.size.height, height=self.size.width),
  )
  for y in 0.. Field2D {
  let result = Field2D::new(
    size=Size::new(width=self.size.height, height=self.size.width),
  )
  for y in 0.. Field2D {
  let result = Field2D::new(size=self.size)
  for y in 0.. Field2D {
  let result = Field2D::new(size=self.size)
  for y in 0.. Field2D {
  let result = Field2D::new(size~)
  for y in 0.. Field2D {
  let b = border.max(0)
  let result = Field2D::new(
    size=Size::new(
      width=self.size.width + b * 2,
      height=self.size.height + b * 2,
    ),
    initial=value,
  )
  for y in 0.. Field2D {
  let side = kernel.length().to_double().sqrt().round().to_int()
  if side <= 0 || side * side != kernel.length() || side % 2 == 0 {
    self.map((_, _, value) => value)
  } else {
    let result = Field2D::new(size=self.size)
    let radius = side / 2
    for y in 0.. Field2D {
  let side = radius.max(0) * 2 + 1
  self.convolve(
    kernel=Array::make(side * side, 1.0 / (side * side).to_double()),
  )
}

///|
/// Compute a radial average profile around a point.
pub fn Field2D::radial_average(
  self : Field2D,
  center~ : Point,
  bins~ : Int,
) -> Array[Double] {
  let n = bins.max(0)
  let sums = Array::make(n, 0.0)
  let counts = Array::make(n, 0)
  let radius = center.distance(
    Point::new(x=self.size.width.to_double(), y=self.size.height.to_double()),
  )
  for y in 0..= 0 && index < n {
        sums[index] += self.get(x, y)
        counts[index] += 1
      }
    }
  }
  let result = Array::make(n, 0.0)
  for i in 0.. Double {
  let mut total = 0.0
  for y in 0.. Double {
  let mut total = 0.0
  for y in 0.. DomainMask {
  let result = DomainMask::new(size=self.size)
  for y in 0..= threshold {
        result.set(x~, y~, solid=true)
      }
    }
  }
  result
}

///|
/// Return a pointwise square of a scalar field.
pub fn Field2D::square(self : Field2D) -> Field2D {
  self.map((_, _, value) => value * value)
}

///|
/// Add a constant to every field value.
pub fn Field2D::add_constant(self : Field2D, amount~ : Double) -> Field2D {
  self.map((_, _, value) => value + amount)
}

///|
/// Scale every field value.
pub fn Field2D::scale_values(self : Field2D, factor~ : Double) -> Field2D {
  self.map((_, _, value) => value * factor)
}

///|
/// Compute a scalar field residual against a constant target.
pub fn Field2D::target_residual(self : Field2D, target~ : Double) -> Double {
  let target_field = Field2D::new(size=self.size, initial=target)
  field_residual(self, target_field)
}

///|
/// Return the mean value of the outermost one-cell ring.
pub fn Field2D::boundary_mean(self : Field2D) -> Double {
  let values = Array::new()
  for y in 0..