Marching squares

From Rosetta Code
Marching squares is a draft programming task. It is not yet considered ready to be promoted as a complete task, for reasons that should be found in its talk page.
Task

Generate contours for a two-dimensional scalar field.

See: Marching squares


11l

Translation of: Wren
-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

Translation of: 11l
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

Translation of: Raku
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

Translation of: Phix
# 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

Library: Wren-seq

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]
Retrieved from "https://rosettacode.org/wiki/Marching_squares?oldid=369870"