Marching squares: Difference between revisions

m
 
(21 intermediate revisions by 9 users not shown)
Line 1:
{{draft task|Image processing}}
 
;Task:
Line 6:
See: [https://en.wikipedia.org/wiki/Marching_squares Marching squares]
 
 
=={{header|11l}}==
{{trans|Wren}}
 
<syntaxhighlight lang="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())
</syntaxhighlight>
 
{{out}}
<pre>
X: 2, Y: -2, Path: [S, S, E, S, E, N, N, N, W, W]
</pre>
 
=={{header|J}}==
 
This is a partial implementation, see the talk page for some discussion of the untouched issues.
 
<syntaxhighlight lang="j">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
}}</syntaxhighlight>
 
Task example:
 
<syntaxhighlight lang="j"> 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│
└───┘</syntaxhighlight>
 
Here, <code>img</code> 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:
 
<syntaxhighlight lang="j">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</syntaxhighlight>
 
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 <code>require'viewmat'</code>) with <code>viewmat Z</code> and <code>viewmat 1 (<"1;C)} 200 200$0</code>, 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.)
 
=={{header|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
<lang ruby>using Contour
<syntaxhighlight lang="julia">using Contour
import GLMakie as GM # GLMakie also defines Contour so import, not using
 
const example = Float64.([
Line 20 ⟶ 212:
])
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
 
</lang>{{out}}
# 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)
</syntaxhighlight>{{out}}
<pre>
points = [(3, 4), (4, 3), (4, 2), (4, 1), (3, 0), (2, 1), (2, 1), (1, 2), (1, 3), (2, 4), (3, 4)]
</pre>
 
=={{header|Lua}}==
Based on the Phix and Wren solutions.
<syntaxhighlight lang="lua">-- 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)
</syntaxhighlight>{{out}}
<pre>
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, }
</pre>
 
=={{header|Perl}}==
{{trans|Raku}}
<syntaxhighlight lang="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);</syntaxhighlight>
{{out}}
<pre>X: 2, Y: -2, Path: SSESENNNWW</pre>
 
=={{header|Phix}}==
Based on the same code as the Wren example.
<!--<lang Phix>(phixonline)-->
<!--<syntaxhighlight lang="phix">(phixonline)-->
<span style="color: #008080;">with</span> <span style="color: #008080;">javascript_semantics</span>
<span style="color: #008080;">enum</span> <span style="color: #000000;">E</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">N</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">W</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">S</span>
Line 56 ⟶ 420:
<span style="color: #008080;">end</span> <span style="color: #008080;">if</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span>
<span style="color: #008080;">switch</span> <span style="color: #000000;">mask</span> <span style="color: #008080;">do</span>
<span style="color: #008080;">case</span> <span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #000000;">5</span><span style="color: #0000FF;">,</span><span style="color: #000000;">13</span> <span style="color: #0000FF;">:</span> <span style="color: #000000;">direction</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">N</span>
<span style="color: #008080;">case</span> <span style="color: #000000;">2</span><span style="color: #0000FF;">,</span><span style="color: #000000;">3</span><span style="color: #0000FF;">,</span><span style="color: #000000;">7</span> <span style="color: #0000FF;">:</span> <span style="color: #000000;">direction</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">E</span>
Line 70 ⟶ 434:
<span style="color: #008080;">if</span> <span style="color: #000000;">cx</span><span style="color: #0000FF;">=</span><span style="color: #000000;">x</span> <span style="color: #008080;">and</span> <span style="color: #000000;">cy</span><span style="color: #0000FF;">=</span><span style="color: #000000;">y</span> <span style="color: #008080;">then</span> <span style="color: #008080;">exit</span> <span style="color: #008080;">end</span> <span style="color: #008080;">if</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">while</span>
<span style="color: #000080;font-style:italic;">-- return 0-based indexes to match other entries</span>
<span style="color: #008080;">return</span> <span style="color: #0000FF;">{</span><span style="color: #000000;">x</span><span style="color: #0000FF;">-</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span> <span style="color: #0000FF;">-(</span><span style="color: #000000;">y</span><span style="color: #0000FF;">-</span><span style="color: #000000;">1</span><span style="color: #0000FF;">),</span> <span style="color: #000000;">directions</span><span style="color: #0000FF;">}</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">if</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span>
<span style="color: #008080;">return</span> <span style="color: #0000FF;">{-</span><span style="color: #000000;">01</span><span style="color: #0000FF;">,-</span><span style="color: #000000;">01</span><span style="color: #0000FF;">,</span><span style="color: #008000;">"Not found!"</span><span style="color: #0000FF;">}</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">function</span>
Line 85 ⟶ 450:
<span style="color: #7060A8;">printf</span><span style="color: #0000FF;">(</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #008000;">"X: %d, Y: %d, Path: %s\n"</span><span style="color: #0000FF;">,</span><span style="color: #000000;">identifyPerimeter</span><span style="color: #0000FF;">(</span><span style="color: #000000;">example</span><span style="color: #0000FF;">))</span>
<!--</langsyntaxhighlight>-->
{{out}}
<pre>
Line 92 ⟶ 457:
 
=={{header|Python}}==
<langsyntaxhighlight lang="python">from numpy import array, round
from skimage import measure
 
Line 107 ⟶ 472:
contours = round(measure.find_contours(example, 0.1))[0]
print('[', ', '.join([str((p[1], 5 - p[0])) for p in contours]), ']')
</langsyntaxhighlight>{{out}}
<pre>
[ (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) ]
</pre>
 
=={{header|Raku}}==
{{trans|Phix}}
<syntaxhighlight lang="raku" line># 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!"
}</syntaxhighlight>
Output is the same as the Phix entry.
 
=={{header|Wren}}==
{{libheader|Wren-seq}}
This is a translation of the [http://www.tomgibara.com/computer-vision/marching-squares 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.
<langsyntaxhighlight ecmascriptlang="wren">import "./seq" for Lst, FrozenList
 
/* A direction in the plane. */
Line 332 ⟶ 744:
var ms = MarchingSquares.new(5, 6, example)
var path = ms.identifyPerimeter()
System.print(path)</langsyntaxhighlight>
 
{{out}}
1,480

edits