Marching squares: Difference between revisions
(added Raku programming solution) |
m (→{{header|Raku}}: insignificant changes) |
||
Line 339:
for ^width X ^height -> (\x,\y) {
unless data[y;x] == 0 {
my
my $mask = 0;
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 so all
}
when * ∈ ( 1, 5, 13 ) { N }
when * ∈ ( 2, 3, 7 ) { E }
Line 354 ⟶ 353:
when * ∈ ( 6 ) { $previous == N ?? W !! E }
when * ∈ ( 9 ) { $previous == E ?? N !! S }
} {
▲ ($cx,$cy) <<+=<< |(@dx[.value], @dy[.value]) }
} until $cx==x and $cy==y ;
return
▲ }
}
}
|
Revision as of 16:04, 9 July 2022
- 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).
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:
<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</lang>
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 <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[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)
</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) ]
Raku
<lang perl6># 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) { unless data[y;x] == 0 { my ($cx,$cy,$directions,$previous) = x, y; repeat { my $mask = 0; 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 so all $mx>1, $my>1, data[$my-1;$mx-1] != 0
} 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!"
}</lang> 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. <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]