Marching squares
Generate contours for a two-dimensional scalar field.
- Task
See: Marching squares
11l
-V
DIR_E = ( 1, 0)
DIR_NE = ( 1, 1)
DIR_N = ( 0, 1)
DIR_NW = (-1, 1)
DIR_W = (-1, 0)
DIR_SW = (-1, -1)
DIR_S = ( 0, -1)
DIR_SE = ( 1, -1)
F path_str(origin_x, origin_y, directions)
-V :dirn = [:DIR_E = ‘E’,
:DIR_NE = ‘NE’,
:DIR_N = ‘N’,
:DIR_NW = ‘NW’,
:DIR_W = ‘W’,
:DIR_SW = ‘SW’,
:DIR_S = ‘S’,
:DIR_SE = ‘SE’]
R ‘X: ’origin_x‘, Y: ’origin_y‘, Path: ’directions.map(d -> @:dirn[d])
T MarchingSquares
. Int width, height
. [Int] data
F (width, height, data)
.width = width
.height = height
.data = copy(data)
. F is_set(x, y)
I x <= 0 | x > .width | y <= 0 | y > .height
R 0B
R .data[(y - 1) * .width + (x - 1)] != 0
. F value(x, y)
V sum = 0
I .is_set(x, y) {sum [|]= 1}
I .is_set(x + 1, y) {sum [|]= 2}
I .is_set(x, y + 1) {sum [|]= 4}
I .is_set(x + 1, y + 1) {sum [|]= 8}
R sum
. F identify_perimeter_(=initial_x, =initial_y)
I initial_x < 0 {initial_x = 0}
I initial_x > .width {initial_x = .width}
I initial_y < 0 {initial_y = 0}
I initial_y > .height {initial_y = .height}
V initial_value = .value(initial_x, initial_y)
I initial_value C (0, 15)
X.throw RuntimeError(‘Supplied initial coordinates (#., #.) do not lie on a perimeter.’.format(initial_x, initial_y))
[IVec2] directions
V x = initial_x
V y = initial_y
V previous = (0, 0)
L
IVec2 direction
S .value(x, y)
1, 5, 13
direction = :DIR_N
2, 3, 7
direction = :DIR_E
4, 12, 14
direction = :DIR_W
8, 10, 11
direction = :DIR_S
6
direction = I previous == :DIR_N {:DIR_W} E :DIR_E
9
direction = I previous == :DIR_E {:DIR_N} E :DIR_S
E
X.throw RuntimeError(‘Illegal state.’)
directions.append(direction)
x += direction.x
y -= direction.y
previous = direction
I x == initial_x & y == initial_y
L.break
R path_str(initial_x, -initial_y, directions)
F identify_perimeter()
V size = .width * .height
L(i) 0 .< size
I .data[i] != 0
R .identify_perimeter_(i % .width, i I/ .width)
V example = [
0, 0, 0, 0, 0,
0, 0, 0, 0, 0,
0, 0, 1, 1, 0,
0, 0, 1, 1, 0,
0, 0, 0, 1, 0,
0, 0, 0, 0, 0
]
V ms = MarchingSquares(5, 6, example)
print(ms.identify_perimeter())
- Output:
X: 2, Y: -2, Path: [S, S, E, S, E, N, N, N, W, W]
FreeBASIC
Type IVec2
x As Integer
y As Integer
End Type
Dim Shared DIR_E As IVec2 = (Type<IVec2>(1, 0))
Dim Shared DIR_NE As IVec2 = (Type<IVec2>(1, 1))
Dim Shared DIR_N As IVec2 = (Type<IVec2>(0, 1))
Dim Shared DIR_NW As IVec2 = (Type<IVec2>(-1, 1))
Dim Shared DIR_W As IVec2 = (Type<IVec2>(-1, 0))
Dim Shared DIR_SW As IVec2 = (Type<IVec2>(-1, -1))
Dim Shared DIR_S As IVec2 = (Type<IVec2>(0, -1))
Dim Shared DIR_SE As IVec2 = (Type<IVec2>(1, -1))
Function path_str(origin_x As Integer, origin_y As Integer, directions() As IVec2) As String
Dim As String result = "X: " & origin_x & ", Y: " & origin_y & ", Path: "
For i As Integer = 1 To Ubound(directions)
result &= "[" & directions(i).x & "," & directions(i).y & "] "
Next
Return result
End Function
Type MarchingSquares
ancho As Integer
alto As Integer
dato(Any) As Integer
Declare Constructor (ancho As Integer, alto As Integer, dato() As Integer)
Declare Function is_set(x As Integer, y As Integer) As Boolean
Declare Function value(x As Integer, y As Integer) As Integer
Declare Function identify_perimeter_(initial_x As Integer, initial_y As Integer) As String
Declare Function identify_perimeter() As String
End Type
Constructor MarchingSquares (ancho As Integer, alto As Integer, dato() As Integer)
This.ancho = ancho
This.alto = alto
Redim This.dato(Ubound(dato))
For i As Integer = 0 To Ubound(dato)
This.dato(i) = dato(i)
Next
End Constructor
Function MarchingSquares.is_set(x As Integer, y As Integer) As Boolean
If x <= 0 Or x > ancho Or y <= 0 Or y > alto Then Return False
Return dato((y - 1) * ancho + (x - 1)) <> 0
End Function
Function MarchingSquares.value(x As Integer, y As Integer) As Integer
Dim As Integer sum = 0
If is_set(x , y ) Then sum Or= 1
If is_set(x + 1, y ) Then sum Or= 2
If is_set(x , y + 1) Then sum Or= 4
If is_set(x + 1, y + 1) Then sum Or= 8
Return sum
End Function
Function MarchingSquares.identify_perimeter_(initial_x As Integer, initial_y As Integer) As String
If initial_x < 0 Then initial_x = 0
If initial_x > ancho Then initial_x = ancho
If initial_y < 0 Then initial_y = 0
If initial_y > alto Then initial_y = alto
Dim As Integer initial_value = value(initial_x, initial_y)
If initial_value = 0 Or initial_value = 15 Then
Print "Supplied initial coordinates (" & initial_x & ", " & initial_y & ") do not lie on a perimeter."
Return ""
End If
Dim As IVec2 directions()
Redim directions(0)
Dim As Integer x = initial_x
Dim As Integer y = initial_y
Dim As IVec2 previous = (Type<IVec2>(0, 0))
Do
Dim As IVec2 direction
Select Case value(x, y)
Case 1, 5, 13
direction = DIR_N
Case 2, 3, 7
direction = DIR_E
Case 4, 12, 14
direction = DIR_W
Case 8, 10, 11
direction = DIR_S
Case 6
direction = Iif(previous.x = DIR_N.x And previous.y = DIR_N.y, DIR_W, DIR_E)
Case 9
direction = Iif(previous.x = DIR_E.x And previous.y = DIR_E.y, DIR_N, DIR_S)
Case Else
Print "Illegal state."
Return ""
End Select
Redim Preserve directions(Ubound(directions) + 1)
directions(Ubound(directions)) = direction
x += direction.x
y -= direction.y
previous = direction
If x = initial_x And y = initial_y Then Exit Do
Loop
Return path_str(initial_x, -initial_y, directions())
End Function
Function MarchingSquares.identify_perimeter() As String
Dim As Integer size = ancho * alto
For i As Integer = 0 To size - 1
If dato(i) <> 0 Then
Return identify_perimeter_(i Mod ancho, i \ ancho)
End If
Next
Return ""
End Function
Dim As Integer example(30) = { _
0, 0, 0, 0, 0, _
0, 0, 0, 0, 0, _
0, 0, 1, 1, 0, _
0, 0, 1, 1, 0, _
0, 0, 0, 1, 0, _
0, 0, 0, 0, 0 }
Dim As MarchingSquares ms = MarchingSquares(5, 6, example())
Print ms.identify_perimeter()
Sleep
- Output:
X: 2, Y: -2, Path: [0,-1] [0,-1] [1,0] [0,-1] [1,0] [0,1] [0,1] [0,1] [-1,0] [-1,0]
J
This is a partial implementation, see the talk page for some discussion of the untouched issues.
step1=: {{2 2 #.@(0 1 3 2&{)@,;._3 ,&0^:2@|:@|.^:4 y}}
step2=: {{(($y)#:i.$y) {{<x step2a y{::LUT}}"1 0 y}}
step2a=: {{ if. #y do. x+"1 y else. y end. }}
LUT=: <@".;._2 {{)n
EMPTY NB. 0
0 0,:1 1 NB. 1
0 1,:1 0 NB. 2
0 0,:0 1 NB. 3
0 0,:1 1 NB. 4
0 1,:1 0 NB. 5
0 0,:1 0 NB. 6
EMPTY NB. 7 0 1,:1 0
1 0,:0 1 NB. 8
1 1,:0 1 NB. 9
0 0,:1 1 NB. 10
EMPTY NB. 11 0 0,:1 1
1 0,:1 1 NB. 12
EMPTY NB. 13 0 1,:1 0
EMPTY NB. 14 0 0,:1 1
EMPTY NB. 15
}}
unwind=: {{
near=. 6 7 8 5 2 3 1 0 {,(+/~ *&({:$y))i:1
r=., c=. EMPTY
TODO=. I.(<EMPTY)~:Y=.,y
j=. _
while.#TODO=. TODO-.j do.
adj=. (j+near) ([-.-.) TODO
if. #adj do.
j=. {.adj
else.
if. #c do. c=.EMPTY [r=. r,<~.c end.
j=. {.TODO
end.
c=. c, j{::Y
end.
r,<~.c
}}
Task example:
img=: 4~:i.3 2
img
1 1
1 1
0 1
unwind step2 step1 img
┌───┐
│1 2│
│2 1│
│3 1│
│4 2│
│5 3│
│4 4│
│3 4│
│2 4│
│1 3│
└───┘
Here, img
is a bitmap. We pad the bitmap by surrounding it with zeros during processing. The box at the end contains a contour corresponding to the bitmap. Here, the first column represents row number (row 0 at the top) and the second column represents column number (column 0 at the left). Each contour represents a closed loop (so the final coordinate pair would connect with the first coordinate pair).
While the above implementation is incomplete, it seems to adequately handle an oval of cassini with focal points at X=1, -1 and Y=0:
a=: 1
X=: |:Y=:201#"0]0.02*i:100
Z0=: (*:(*:X)+(*:Y)) + (_2*(*:a)*X -&*: Y) + *:a
Z=: (Z0>:0.8)*Z0<:1.2
C=: unwind step2 step1 Z
Here, Z is a bitmap, and C is a list of three contours (one with 336 distinct coordinate pairs, and two with 134 distinct coordinate pairs) which surround that bitmap. These can be inspected (after require'viewmat'
) with viewmat Z
and viewmat 1 (<"1;C)} 200 200$0
, and look plausible.
(Presumably the above implementation would fail if a threshold had been picked such that the bitmap exhibited a saddlepoint at the origin.)
Julia
Uses the marching squares algorithm: see github.com/JuliaGeometry/Contour.jl/blob/master/src/Contour.jl See the discussion page for the Oval of Cassini example
using Contour
import GLMakie as GM # GLMakie also defines Contour so import, not using
const example = Float64.([
0 0 0 0 0;
0 0 0 0 0;
0 0 1 1 0;
0 0 1 1 0;
0 0 0 1 0;
0 0 0 0 0;
])
const cl = first(levels(contours(1:6, 1:5, example)))
const xs, ys = coordinates(first(lines(cl)))
# Showing the points of the contour as origin (0, 0) at bottom left
const points = [(Int(round(ys[i])) - 1, 6 - Int(round(xs[i]))) for i in eachindex(xs)]
@show points
# oval of Cassini formula in z below, see formula at en.wikipedia.org/wiki/Cassini_oval#Equations
const xarray, yarray, a = -2.0:0.02:2.0, -2.0:0.02:2.0, 1.0
const z = [isapprox((x^2 + y^2)^2 - 2 * a^2 * (x^2 - y^2) + a^2, 1.0, atol=0.2) ? 1.0 : 0.0
for x in xarray, y in yarray]
# The first (and pehaps only significant) level is the 0 <-> 1 transition border
# There are 3 separate contours at that level, on outside and 2 holes
const figeight = levels(contours(1:size(z, 1), 1:size(z, 2), z))
const ovalxs, ovalys = coordinates(lines(figeight[1])[1])
const ovalxs2, ovalys2 = coordinates(lines(figeight[1])[2])
const ovalxs3, ovalys3 = coordinates(lines(figeight[1])[3])
const oplot = GM.plot(z)
GM.lines!(ovalxs, ovalys, color = :red, linewidth = 4)
GM.lines!(ovalxs2, ovalys2, color = :white, linewidth = 4)
GM.lines!(ovalxs3, ovalys3, color = :lightgreen, linewidth = 4)
GM.display(oplot)
- Output:
points = [(3, 4), (4, 3), (4, 2), (4, 1), (3, 0), (2, 1), (2, 1), (1, 2), (1, 3), (2, 4), (3, 4)]
Lua
Based on the Phix and Wren solutions.
-- positive directions: right, down, clockwise
local Directions = { -- clockwise from North
N = {x= 0, y=-1},
E = {x= 1, y= 0},
S = {x= 0, y= 1},
W = {x=-1, y= 0},
}
local dxdybList = {
{0,0,1}, -- same position
{1,0,2}, -- right
{0,1,4}, -- down
{1,1,8}, -- right and down
}
local cases = {
"W", "N", "W", "S",
"S", nil, "S", "E",
nil, "N", "W", "E",
"E", "N",
}
local function identifyPerimeter(iLayer, startingX, startingY, data)
local resultDirections = {}
local resultPositions = {}
local currentX, currentY = startingX, startingY
local direction, prevDirection
while true do
local mask = 0
for _, d in ipairs (dxdybList) do
local dx, dy, b = d[1], d[2], d[3]
local mx, my = currentX+dx, currentY+dy
if mx>1 and my>1
and data[my-1] and data[my-1][mx-1]
and data[my-1][mx-1] == iLayer then
mask = mask + b
end
end
direction = cases[mask]
if not direction then
if mask == 6 then
direction = (prevDirection == "E") and "N" or "S"
elseif mask == 9 then
direction = (prevDirection == "S") and "E" or "W"
else
error ('no mask: '.. mask..' by x:'..currentX..' y:'..currentY, 1)
end
end
table.insert (resultDirections, direction)
table.insert (resultPositions, currentX)
table.insert (resultPositions, currentY)
local vDir = Directions[direction]
currentX, currentY = currentX+vDir.x, currentY+vDir.y
prevDirection = direction
if startingX == currentX and startingY == currentY then
return resultDirections, resultPositions
end
end
end
local function findFirstOnLayer (iLayer, data)
for y = 1, #data do -- from 1 to hight
for x = 1, #data[1] do -- from 1 to width
local value = data[y][x]
if value == iLayer then
return x, y -- only one contour
end
end
end
end
local function msMain (iLayer, data)
local rootX, rootY = findFirstOnLayer (iLayer, data)
print ("root: x="..rootX..' y='..rootY)
local directions, positions = identifyPerimeter(iLayer, rootX, rootY, data)
print ('directions amount: '..#directions)
print ('directions: '.. table.concat (directions, ','))
local strPositions = ""
for i = 1, #positions-1, 2 do
strPositions = strPositions..positions[i]..','..positions[i+1]..', '
end
print ('positions: {' .. strPositions..'}')
end
local example = {
{0, 0, 0, 0, 0, 0},
{1, 0, 0, 0, 0, 1},
{0, 1, 1, 0, 1, 0},
{0, 1, 1, 1, 1, 0},
{0, 1, 0, 1, 1, 0},
{1, 0, 0, 0, 0, 1},
}
msMain (1, example)
- Output:
root: x=1 y=2 directions amount: 34 directions: E,S,E,E,S,E,N,E,N,E,S,W,S,S,S,E,S,W,N,W,W,N,W,S,W,S,W,N,E,N,N,N,W,N positions: {1,2, 2,2, 2,3, 3,3, 4,3, 4,4, 5,4, 5,3, 6,3, 6,2, 7,2, 7,3, 6,3, 6,4, 6,5, 6,6, 7,6, 7,7, 6,7, 6,6, 5,6, 4,6, 4,5, 3,5, 3,6, 2,6, 2,7, 1,7, 1,6, 2,6, 2,5, 2,4, 2,3, 1,3, }
Perl
use v5.36;
no warnings 'experimental::for_list';
use List::Util 'any';
use enum <E N W S>;
sub X ($a,$b) { my @c; for my $aa (0..$a) { for my $bb (0..$b) { push @c, $aa, $bb } } @c }
sub identify_perimeter(@data) {
for my ($x,$y) (X $#{$data[0]}, $#data) {
next unless $data[$y][$x] and $data[$y][$x] != 0;
my ($path,$cx,$cy,$d,$p) = ('', $x, $y);
do {
my $mask;
for my($dx,$dy,$b) (0,0,1, 1,0,2, 0,1,4, 1,1,8) {
my ($mx, $my) = ($cx+$dx, $cy+$dy);
$mask += $b if $mx>1 and $my>1 and $data[$my-1][$mx-1] != 0
}
$d = N if any { $mask == $_ } (1, 5,13);
$d = E if any { $mask == $_ } (2, 3, 7);
$d = W if any { $mask == $_ } (4,12,14);
$d = S if any { $mask == $_ } (8,10,11);
$d = $p == N ? W : E if $mask == 6;
$d = $p == E ? N : S if $mask == 9;
$path .= $p = (<E N W S>)[$d];
$cx += (1, 0,-1,0)[$d];
$cy += (0,-1, 0,1)[$d];
} until $cx == $x and $cy == $y;
return $x, -$y, $path
}
exit 'That did not work out...';
}
my @M = ([0, 0, 0, 0, 0],
[0, 0, 0, 0, 0],
[0, 0, 1, 1, 0],
[0, 0, 1, 1, 0],
[0, 0, 0, 1, 0],
[0, 0, 0, 0, 0]);
printf "X: %d, Y: %d, Path: %s\n", identify_perimeter(@M);
- Output:
X: 2, Y: -2, Path: SSESENNNWW
Phix
Based on the same code as the Wren example.
with javascript_semantics enum E, N, W, S constant dx = {1,0,-1,0}, dy = {0,-1,0,1} function identifyPerimeter(sequence data) integer height = length(data), width = length(data[1]) for x=1 to width do for y=1 to height do if data[y][x]!=0 then string directions = "" integer cx = x, cy = y, direction, previous = null; while true do integer mask = 0 for dxyb in {{0,0,1},{1,0,2},{0,1,4},{1,1,8}} do integer {dx,dy,b} = dxyb, mx = cx+dx, my = cy+dy if mx>1 and my>1 and data[my-1,mx-1]!=0 then mask += b end if end for switch mask case 1,5,13 : direction = N case 2,3,7 : direction = E case 4,12,14: direction = W case 8,10,11: direction = S case 6: direction = iff(previous == N ? W : E) case 9: direction = iff(previous == E ? N : S) end switch directions &= "ENWS"[direction] cx += dx[direction] cy += dy[direction] previous = direction if cx=x and cy=y then exit end if end while -- return 0-based indexes to match other entries return {x-1, -(y-1), directions} end if end for end for return {-1,-1,"Not found!"} end function constant example = {{0, 0, 0, 0, 0}, {0, 0, 0, 0, 0}, {0, 0, 1, 1, 0}, {0, 0, 1, 1, 0}, {0, 0, 0, 1, 0}, {0, 0, 0, 0, 0}} printf(1,"X: %d, Y: %d, Path: %s\n",identifyPerimeter(example))
- Output:
X: 2, Y: -2, Path: SSESENNNWW
Python
from numpy import array, round
from skimage import measure
example = array([
[0, 0, 0, 0, 0],
[0, 0, 0, 0, 0],
[0, 0, 1, 1, 0],
[0, 0, 1, 1, 0],
[0, 0, 0, 1, 0],
[0, 0, 0, 0, 0]
])
# Find contours at a constant value of 0.1 and extract the first one found
contours = round(measure.find_contours(example, 0.1))[0]
print('[', ', '.join([str((p[1], 5 - p[0])) for p in contours]), ']')
- Output:
[ (3.0, 0.0), (2.0, 1.0), (2.0, 1.0), (1.0, 2.0), (1.0, 3.0), (2.0, 4.0), (3.0, 4.0), (4.0, 3.0), (4.0, 2.0), (4.0, 1.0), (3.0, 0.0) ]
Raku
# 20220708 Raku programming solution
enum <E N W S>;
my (@dx,@dy) := (1,0,-1,0), (0,-1,0,1);
my \example = ((0, 0, 0, 0, 0),
(0, 0, 0, 0, 0),
(0, 0, 1, 1, 0),
(0, 0, 1, 1, 0),
(0, 0, 0, 1, 0),
(0, 0, 0, 0, 0));
printf("X: %d, Y: %d, Path: %s\n", identifyPerimeter(example));
sub identifyPerimeter(\data) {
my (\height,\width) = { .elems, .first.elems }(data);
for ^width X ^height -> (\x,\y) {
if data[y;x] {
my ($cx,$cy,$directions,$previous) = x, y;
repeat {
my $mask;
for (0,0,1),(1,0,2),(0,1,4),(1,1,8) -> (\dx,\dy,\b) {
my ($mx,$my) = $cx+dx,$cy+dy;
$mask += b if all $mx>1, $my>1, data[$my-1;$mx-1]
}
given do given $mask {
when * ∈ ( 1, 5, 13 ) { N }
when * ∈ ( 2, 3, 7 ) { E }
when * ∈ ( 4, 12, 14 ) { W }
when * ∈ ( 8, 10, 11 ) { S }
when * ∈ ( 6 ) { $previous == N ?? W !! E }
when * ∈ ( 9 ) { $previous == E ?? N !! S }
} {
$directions ~= $previous = $_ ;
($cx,$cy) <<+=<< (@dx[.value], @dy[.value])
}
} until $cx==x and $cy==y ;
return x, -y, $directions
}
}
return -1, -1, "Not found!"
}
Output is the same as the Phix entry.
Wren
This is a translation of the public domain Java code, written by Tom Gibara, which is linked to from the Wikipedia article. It also uses his example to test the code.
import "./seq" for Lst, FrozenList
/* A direction in the plane. */
class Direction {
// statics
static E { new_( 1, 0) }
static NE { new_( 1, 1) }
static N { new_( 0, 1) }
static NW { new_(-1, 1) }
static W { new_(-1, 0) }
static SW { new_(-1, -1) }
static S { new_( 0, -1) }
static SE { new_( 1, -1) }
// private constructor
construct new_(x, y) {
_planeX = x
_planeY = y
_screenX = x
_screenY = -y
_length = (x != 0 && y != 0) ? 2.sqrt/2 : 1
}
// property getters
planeX { _planeX } // horizontal distance moved in this direction within the plane
planeY { _planeY } // vertical distance moved in this direction within the plane
screenX { _screenX } // horizontal distance moved in this direction in screen coordinates
screenY { _screenY } // vertical distance moved in this direction in screen coordinates
length { _length } // euclidean length of this direction's vectors
// equality override
==(that) {
if (Object.same(this, that)) return true
return _planeX == that.planeX && _planeY == that.planeY &&
_screenX == that.screenX && _screenY == that.screenY &&
_length == that.length
}
// string representation
toString {
if (this == Direction.E) return "E"
if (this == Direction.NE) return "NE"
if (this == Direction.N) return "N"
if (this == Direction.NW) return "NW"
if (this == Direction.W) return "W"
if (this == Direction.SW) return "SW"
if (this == Direction.S) return "S"
if (this == Direction.SE) return "SE"
return ""
}
}
/* Combines a sequence of directions into a path that is rooted at some point in the plane.
No restrictions are placed on Path objects which are immutable. */
class Path {
// static
static ADJ_LEN { 2.sqrt/2 - 1 }
// public constructor
construct new(startX, startY, directions) {
_originX = startX
_originY = startY
_directions = Lst.clone(directions)
_directionList = FrozenList.new(directions)
var endX = startX
var endY = startY
var diagonals = 0
for (direction in directions) {
endX = endX + direction.screenX
endY = endY + direction.screenY
if (direction.screenX != 0 && direction.screenY != 0) {
diagonals = diagonals + 1
}
}
_terminalX = endX
_terminalY = endY
_length = directions.count + diagonals * Path.ADJ_LEN
}
// private constructor
construct new_(that, deltaX, deltaY) {
_directions = that.directions
_directionList = that.directionList
_length = that.length
_originX = that.originX + deltaX
_originY = that.originY + deltaY
_terminalX = that.terminalX + deltaX
_terminalY = that.terminalY + deltaY
}
// property getters
directions { _directionList } // immutable list of directions that compose this path
originX { _originX } // x coordinate in the plane at which the path begins
originY { _originY } // y coordinate in the plane at which the path begins
terminalX { _terminalX } // x coordinate in the plane at which the path ends
terminalY { _terminalY } // y coordinate in the plane at which the path ends
length { _length } // length of the path using the standard Euclidean metric
// returns whether the path's point of origin is the same as its point of termination
isClosed { _originX == _terminalX && _originY == _terminalY }
// creates a new Path by translating this path in the plane.
translate(deltaX, deltaY) { Path.new_(this, deltaX, deltaY) }
// equals override
==(that) {
if (Object.same(this, that)) return true
if (!(that is Path)) return false
if (_originX != that.originX) return false
if (_originY != that.originY) return false
if (_terminalX != that.terminalX) return false
if (_terminalY != that.terminalY) return false
if (!Lst.areEqual(_directions, that.directions)) return false
return true
}
// string representation
toString { "X: %(originX), Y: %(originY), Path: %(_directions)" }
}
/* A simple implementation of the marching squares algorithm that can identify
perimeters in a supplied byte array. */
class MarchingSquares {
// constructor
construct new(width, height, data) {
_width = width
_height = height
_data = data // not copied but should not be changed
}
// property getters
width { _width } // width of the data matrix
height { _height } // height of the data matrix
data { _data } // data matrix
/* methods */
// finds the perimeter between a set of zero and non-zero values which
// begins at the specified data element - always returns a closed path
identifyPerimeter(initialX, initialY) {
if (initialX < 0) initialX = 0
if (initialX > _width) initialX = _width
if (initialY < 0) initialY = 0
if (initialY > _height) initialY = _height
var initialValue = value_(initialX, initialY)
if (initialValue == 0 || initialValue == 15) {
Fiber.abort("Supplied initial coordinates (%(initialX), %(initialY) " +
"do not lie on a perimeter.")
}
var directions = []
var x = initialX
var y = initialY
var previous = null
while (true) {
var direction
var v = value_(x, y)
if (v == 1 || v == 5 || v == 13) {
direction = Direction.N
} else if (v == 2 || v == 3 || v == 7) {
direction = Direction.E
} else if (v == 4 || v == 12 || v == 14) {
direction = Direction.W
} else if (v == 8 || v == 10 || v == 11) {
direction = Direction.S
} else if (v == 6) {
direction = (previous == Direction.N) ? Direction.W : Direction.E
} else if (v == 9) {
direction = (previous == Direction.E) ? Direction.N : Direction.S
} else {
Fiber.abort("Illegal state.")
}
directions.add(direction)
x = x + direction.screenX
y = y + direction.screenY
previous = direction
if (x == initialX && y == initialY) break
}
return Path.new(initialX, -initialY, directions)
}
// convenience version of above method to be used where no initial point is known
// returns null if there is no perimeter
identifyPerimeter() {
var size = width * height
for (i in 0...size) {
if (_data[i] != 0) return identifyPerimeter(i % _width, (i / _width).floor)
}
return null
}
// private utility methods
value_(x, y) {
var sum = 0
if (isSet_(x, y)) sum = sum | 1
if (isSet_(x+1, y)) sum = sum | 2
if (isSet_(x, y+1)) sum = sum | 4
if (isSet_(x+1, y+1)) sum = sum | 8
return sum
}
isSet_(x, y) {
return (x <= 0 || x > width || y <= 0 || y > height) ? false :
_data[(y - 1) * width + (x - 1)] != 0
}
}
var example = [
0, 0, 0, 0, 0,
0, 0, 0, 0, 0,
0, 0, 1, 1, 0,
0, 0, 1, 1, 0,
0, 0, 0, 1, 0,
0, 0, 0, 0, 0
]
var ms = MarchingSquares.new(5, 6, example)
var path = ms.identifyPerimeter()
System.print(path)
- Output:
X: 2, Y: -2, Path: [S, S, E, S, E, N, N, N, W, W]