I'm working on modernizing Rosetta Code's infrastructure. Starting with communications. Please accept this time-limited open invite to RC's Slack.. --Michael Mol (talk) 20:59, 30 May 2020 (UTC)

# Marching squares

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.

Generate contours for a two-dimensional scalar field.

See: Marching squares

## J

This is a partial implementation, see the talk page for some discussion of the untouched issues.

`step1=: {{2 2 #[email protected](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}}`

`   img=: 4~:i.3 2   img1 11 10 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=: 1X=: |:Y=:201#"0]0.02*i:100Z0=: (*:(*:X)+(*:Y)) + (_2*(*:a)*X -&*: Y) + *:aZ=: (Z0>:0.8)*Z0<:1.2C=: 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 Contourimport 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 leftconst 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#Equationsconst xarray, yarray, a = -2.0:0.02:2.0, -2.0:0.02:2.0, 1.0const 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 holesconst figeight = levels(contours(1:size(z, 1), 1:size(z, 2), z))const ovalxs, ovalys = coordinates(lines(figeight))const ovalxs2, ovalys2 = coordinates(lines(figeight))const ovalxs3, ovalys3 = coordinates(lines(figeight)) 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, clockwiselocal 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, d, d			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	endend local function findFirstOnLayer (iLayer, data)	for y = 1, #data do -- from 1 to hight		for x = 1, #data do -- from 1 to width			local value = data[y][x]			if value == iLayer then				return x, y -- only one contour			end		end	endend 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}, \$#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)
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
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
end if
end for
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
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, roundfrom 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 foundcontours = round(measure.find_contours(example, 0.1))print('[', ', '.join([str((p, 5 - p)) 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]
```