///|
/// A single interpolated segment emitted by a marching-squares pass.
pub(all) struct ContourSegment {
start : Vec2
end : Vec2
level : Double
cell_x : Int
cell_y : Int
} derive(Eq, Debug, ToJson)
///|
/// A stitched isoline with its source scalar level.
pub(all) struct ContourLine {
level : Double
points : Array[Vec2]
closed : Bool
length : Double
bounds : Rect
} derive(Eq, Debug, ToJson)
///|
/// A filled band between two neighboring contour levels.
pub(all) struct ContourBand {
low : Double
high : Double
polygons : Array[Array[Vec2]]
area : Double
} derive(Eq, Debug, ToJson)
///|
/// A compact field derivative used for shading and flow alignment.
pub(all) struct FieldGradient {
value : Double
dx : Double
dy : Double
magnitude : Double
direction : Double
} derive(Eq, Debug, ToJson)
///|
fn field_value(field : ScalarField, x : Int, y : Int) -> Double {
field.get(x, y).unwrap_or(0.0)
}
///|
fn interpolate_scalar(left : Double, right : Double, level : Double) -> Double {
let difference = right - left
if difference.abs() < 0.0000001 {
0.5
} else {
(level - left) / difference.max(-1000000.0).min(1000000.0)
}
}
///|
fn edge_point(
x : Int,
y : Int,
edge : Int,
top : Double,
right : Double,
bottom : Double,
left : Double,
level : Double,
) -> Vec2 {
let amount = match edge {
0 => interpolate_scalar(top, right, level)
1 => interpolate_scalar(right, bottom, level)
2 => interpolate_scalar(bottom, left, level)
_ => interpolate_scalar(left, top, level)
}
match edge {
0 => Vec2::new(x.to_double() + amount, y.to_double())
1 => Vec2::new(x.to_double() + 1.0, y.to_double() + amount)
2 => Vec2::new(x.to_double() + 1.0 - amount, y.to_double() + 1.0)
_ => Vec2::new(x.to_double(), y.to_double() + 1.0 - amount)
}
}
///|
fn marching_case(
top : Double,
right : Double,
bottom : Double,
left : Double,
level : Double,
) -> Int {
let mut result = 0
if top >= level {
result += 1
}
if right >= level {
result += 2
}
if bottom >= level {
result += 4
}
if left >= level {
result += 8
}
result
}
///|
fn add_contour_pair(
result : Array[ContourSegment],
x : Int,
y : Int,
level : Double,
first : Int,
second : Int,
top : Double,
right : Double,
bottom : Double,
left : Double,
) -> Unit {
result.push({
start: edge_point(x, y, first, top, right, bottom, left, level),
end: edge_point(x, y, second, top, right, bottom, left, level),
level,
cell_x: x,
cell_y: y,
})
}
///|
/// Emits all line segments for one scalar cell.
pub fn marching_cell_segments(
field : ScalarField,
x : Int,
y : Int,
level : Double,
) -> Array[ContourSegment] {
let result = []
let top = field_value(field, x, y)
let right = field_value(field, x + 1, y)
let bottom = field_value(field, x + 1, y + 1)
let left = field_value(field, x, y + 1)
match marching_case(top, right, bottom, left, level) {
0 | 15 => ()
1 => add_contour_pair(result, x, y, level, 3, 0, top, right, bottom, left)
2 => add_contour_pair(result, x, y, level, 0, 1, top, right, bottom, left)
3 => add_contour_pair(result, x, y, level, 3, 1, top, right, bottom, left)
4 => add_contour_pair(result, x, y, level, 1, 2, top, right, bottom, left)
5 => {
let center = (top + right + bottom + left) / 4.0
if center >= level {
add_contour_pair(result, x, y, level, 3, 2, top, right, bottom, left)
add_contour_pair(result, x, y, level, 0, 1, top, right, bottom, left)
} else {
add_contour_pair(result, x, y, level, 3, 0, top, right, bottom, left)
add_contour_pair(result, x, y, level, 2, 1, top, right, bottom, left)
}
}
6 => add_contour_pair(result, x, y, level, 0, 2, top, right, bottom, left)
7 => add_contour_pair(result, x, y, level, 3, 2, top, right, bottom, left)
8 => add_contour_pair(result, x, y, level, 2, 3, top, right, bottom, left)
9 => add_contour_pair(result, x, y, level, 0, 2, top, right, bottom, left)
10 => {
let center = (top + right + bottom + left) / 4.0
if center >= level {
add_contour_pair(result, x, y, level, 3, 0, top, right, bottom, left)
add_contour_pair(result, x, y, level, 1, 2, top, right, bottom, left)
} else {
add_contour_pair(result, x, y, level, 0, 1, top, right, bottom, left)
add_contour_pair(result, x, y, level, 2, 3, top, right, bottom, left)
}
}
11 => add_contour_pair(result, x, y, level, 1, 2, top, right, bottom, left)
12 => add_contour_pair(result, x, y, level, 3, 1, top, right, bottom, left)
13 => add_contour_pair(result, x, y, level, 0, 1, top, right, bottom, left)
_ => add_contour_pair(result, x, y, level, 3, 0, top, right, bottom, left)
}
result
}
///|
/// Extracts marching-squares segments across an entire scalar field.
pub fn marching_squares(
field : ScalarField,
level : Double,
) -> Array[ContourSegment] {
let result = []
if field.width > 1 && field.height > 1 {
for y in 0..<(field.height - 1) {
for x in 0..<(field.width - 1) {
for segment in marching_cell_segments(field, x, y, level) {
result.push(segment)
}
}
}
}
result
}
///|
fn close_enough(left : Vec2, right : Vec2, tolerance : Double) -> Bool {
left.distance(right) <= tolerance.max(0.000001)
}
///|
fn find_line_endpoint(
lines : Array[Array[Vec2]],
point : Vec2,
tolerance : Double,
) -> (Int, Bool)? {
for index, line in lines {
if line.length() > 0 {
if close_enough(line[0], point, tolerance) {
return Some((index, true))
}
if close_enough(line[line.length() - 1], point, tolerance) {
return Some((index, false))
}
}
}
None
}
///|
fn append_segment_to_lines(
lines : Array[Array[Vec2]],
segment : ContourSegment,
tolerance : Double,
) -> Unit {
let first = find_line_endpoint(lines, segment.start, tolerance)
let second = find_line_endpoint(lines, segment.end, tolerance)
match (first, second) {
(None, None) => lines.push([segment.start, segment.end])
(Some((index, at_start)), None) =>
if at_start {
lines[index].insert(0, segment.end)
} else {
lines[index].push(segment.end)
}
(None, Some((index, at_start))) =>
if at_start {
lines[index].insert(0, segment.start)
} else {
lines[index].push(segment.start)
}
(Some((first_index, first_start)), Some((second_index, second_start))) =>
if first_index == second_index {
if !close_enough(segment.start, segment.end, tolerance) {
if first_start {
lines[first_index].insert(0, segment.end)
} else {
lines[first_index].push(segment.end)
}
}
} else {
let first_line = lines[first_index].copy()
let second_line = lines[second_index].copy()
let merged : Array[Vec2] = []
if first_start {
for point in first_line.rev() {
merged.push(point)
}
} else {
for point in first_line {
merged.push(point)
}
}
if second_start {
for point in second_line {
merged.push(point)
}
} else {
for point in second_line.rev() {
merged.push(point)
}
}
lines[first_index] = merged
let _ = lines.remove(second_index)
}
}
}
///|
/// Stitches nearby marching-squares segments into ordered polylines.
pub fn stitch_contour_segments(
segments : Array[ContourSegment],
tolerance? : Double = 0.01,
) -> Array[Array[Vec2]] {
let lines : Array[Array[Vec2]] = []
for segment in segments {
append_segment_to_lines(lines, segment, tolerance)
}
lines
}
///|
/// Converts one polyline into a measured contour line.
pub fn contour_line(
level : Double,
points : Array[Vec2],
closed? : Bool = false,
) -> ContourLine {
let cleaned = deduplicate_polyline(points, 0.000001)
{
level,
points: cleaned,
closed,
length: polyline_length(cleaned, closed~),
bounds: bounds_of_points(cleaned),
}
}
///|
/// Extracts and stitches all isolines at one level.
pub fn isolines(
field : ScalarField,
level : Double,
tolerance? : Double = 0.01,
) -> Array[ContourLine] {
let result = []
for
points in stitch_contour_segments(
marching_squares(field, level),
tolerance~,
) {
let closed = points.length() > 2 &&
close_enough(points[0], points[points.length() - 1], tolerance)
let normalized = if closed {
let copy = []
for index in 0..<(points.length() - 1) {
copy.push(points[index])
}
copy
} else {
points
}
result.push(contour_line(level, normalized, closed~))
}
result
}
///|
/// Extracts isolines for every requested level in stable order.
pub fn field_isolines(
field : ScalarField,
levels : Array[Double],
tolerance? : Double = 0.01,
) -> Array[ContourLine] {
let result = []
for level in levels {
for line in isolines(field, level, tolerance~) {
result.push(line)
}
}
result
}
///|
/// Computes a central-difference gradient at a scalar field coordinate.
pub fn scalar_field_gradient(
field : ScalarField,
x : Int,
y : Int,
) -> FieldGradient {
let value = field_value(field, x, y)
let dx = (field_value(field, x + 1, y) - field_value(field, x - 1, y)) / 2.0
let dy = (field_value(field, x, y + 1) - field_value(field, x, y - 1)) / 2.0
{
value,
dx,
dy,
magnitude: (dx * dx + dy * dy).sqrt(),
direction: @math.atan2(dy, dx),
}
}
///|
/// Samples field gradients on a regular lattice.
pub fn scalar_field_gradients(
field : ScalarField,
stride : Int,
) -> Array[FieldGradient] {
let result = []
let step = stride.max(1)
for y in 0.. ScalarField {
let result = ScalarField::new(field.width, field.height)
for y in 0.. ScalarField {
let result = field.normalize()
result.map(fn(value) { low + (high - low) * value })
}
///|
/// Counts cells above a threshold, useful for occupancy diagnostics.
pub fn scalar_field_occupancy(field : ScalarField, threshold : Double) -> Int {
let mut count = 0
for y in 0..= threshold {
count += 1
}
}
}
count
}
///|
/// Converts contour lines into open or closed SVG paths.
pub fn contour_paths(lines : Array[ContourLine]) -> Array[Path] {
let result = []
for line in lines {
result.push(polyline(line.points, closed=line.closed))
}
result
}
///|
/// Creates an SVG document containing one style per contour level.
pub fn contour_document(
width : Double,
height : Double,
lines : Array[ContourLine],
palette : Palette,
stroke_width : Double,
) -> SvgDocument {
let document = SvgDocument::new(width, height)
for index, line in lines {
let color = palette
.at(index % palette.colors.length())
.unwrap_or(Color::rgb8(0, 0, 0))
document.add(
SvgElement::SvgPath(
polyline(line.points, closed=line.closed),
SvgStyle::new(fill="none", stroke=color.to_hex(), stroke_width~),
),
)
}
document
}
///|
/// Produces filled contour bands by rasterizing consecutive thresholds.
pub fn contour_bands(
field : ScalarField,
levels : Array[Double],
) -> Array[ContourBand] {
let sorted = levels.copy()
for i in 0..= 2 {
for index in 0..<(sorted.length() - 1) {
let low = sorted[index]
let high = sorted[index + 1]
let mask = field.map(fn(value) {
if value >= low && value < high {
1.0
} else {
0.0
}
})
let polygons = []
for component in raster_components_to_polygons(mask, 0.5) {
polygons.push(component)
}
result.push({
low,
high,
polygons,
area: polygons.fold(init=0.0, fn(total, polygon) {
total + polygon_area(polygon).abs()
}),
})
}
}
result
}
///|
/// Converts connected field cells into rectangular component polygons.
pub fn raster_components_to_polygons(
field : ScalarField,
threshold : Double,
) -> Array[Array[Vec2]] {
let result = []
for component in field_to_components(field, threshold) {
if component.length() > 0 {
let bounds = bounds_of_points(component)
result.push([
Vec2::new(bounds.x, bounds.y),
Vec2::new(bounds.right(), bounds.y),
Vec2::new(bounds.right(), bounds.bottom()),
Vec2::new(bounds.x, bounds.bottom()),
])
}
}
result
}
///|
/// Finds connected active cells without constructing a raster object.
pub fn field_to_components(
field : ScalarField,
threshold : Double,
) -> Array[Array[Vec2]] {
let visited = []
for _ in 0..<(field.width * field.height) {
visited.push(false)
}
let result = []
for y in 0..= threshold {
let queue : Array[(Int, Int)] = [(x, y)]
let component : Array[Vec2] = []
let mut head = 0
visited[index] = true
while head < queue.length() {
let current = queue[head]
head += 1
component.push(
Vec2::new(current.0.to_double(), current.1.to_double()),
)
for
neighbor in [
(current.0 + 1, current.1),
(current.0 - 1, current.1),
(current.0, current.1 + 1),
(current.0, current.1 - 1),
] {
if field.in_bounds(neighbor.0, neighbor.1) {
let ni = field.index(neighbor.0, neighbor.1)
if !visited[ni] &&
field_value(field, neighbor.0, neighbor.1) >= threshold {
visited[ni] = true
queue.push(neighbor)
}
}
}
}
result.push(component)
}
}
}
result
}
///|
/// Simplifies all contour lines by a distance tolerance.
pub fn simplify_contours(
lines : Array[ContourLine],
tolerance : Double,
) -> Array[ContourLine] {
let result = []
for line in lines {
let points = simplify_sampled_polyline(line.points, tolerance)
result.push(contour_line(line.level, points, closed=line.closed))
}
result
}
///|
/// Calculates a scalar summary for a contour collection.
pub fn contour_summary(lines : Array[ContourLine]) -> String {
let mut total_length = 0.0
let mut closed_count = 0
let mut point_count = 0
for line in lines {
total_length += line.length
point_count += line.points.length()
if line.closed {
closed_count += 1
}
}
"lines=\{lines.length()},closed=\{closed_count},points=\{point_count},length=\{total_length}"
}