Marching squares
- Task
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.
<lang J>step1=: Template:2 2 step2=: Template:(($y) step2a=: Template:If. 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
}}</lang>
Task example:
<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│ └───┘</lang>
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).
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 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[2])[2]) const ovalxs3, ovalys3 = coordinates(lines(figeight[2])[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)
</lang>
- 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. <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)
</lang>
- 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, }
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
<lang 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]), ']')
</lang>
- 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) ]
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. <lang ecmascript>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)</lang>
- Output:
X: 2, Y: -2, Path: [S, S, E, S, E, N, N, N, W, W]