Marching squares: Difference between revisions

From Rosetta Code
Content added Content deleted
m (→‎{{header|J}}: not quite so lame, but still lame)
(Oval plot)
Line 75: Line 75:
=={{header|Julia}}==
=={{header|Julia}}==
Uses the marching squares algorithm: see github.com/JuliaGeometry/Contour.jl/blob/master/src/Contour.jl
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
<lang ruby>using Contour
import GLMakie as GM # GLMakie also defines Contour so import, not using


const example = Float64.([
const example = Float64.([
Line 86: Line 88:
])
])
cl = first(levels(contours(1:6, 1:5, example)))
const cl = first(levels(contours(1:6, 1:5, example)))
xs, ys = coordinates(first(lines(cl)))
const xs, ys = coordinates(first(lines(cl)))


# Showing the points of the contour as origin (0, 0) at bottom left
# Showing the points of the contour as origin (0, 0) at bottom left
points = [(Int(round(ys[i])) - 1, 6 - Int(round(xs[i]))) for i in eachindex(xs)]
const points = [(Int(round(ys[i])) - 1, 6 - Int(round(xs[i]))) for i in eachindex(xs)]
@show points
@show points


# oval of Cassini formula in z below, see formula at en.wikipedia.org/wiki/Cassini_oval#Equations
# Plotting
const xarray, yarray, a = -2.0:0.02:2.0, -2.0:0.02:2.0, 1.0
using GLMakie
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
plt = plot(example)
for x in xarray, y in yarray]
lines!(xs, ys, linewidth = 8, color = :red)
⚫
display(plt)


# 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>{{out}}
</lang>{{out}}
<pre>
<pre>
points = [(3, 4), (4, 3), (4, 2), (4, 1), (3, 0), (2, 1), (2, 1), (1, 2), (1, 3), (2, 4), (3, 4)]
points = [(3, 4), (4, 3), (4, 2), (4, 1), (3, 0), (2, 1), (2, 1), (1, 2), (1, 3), (2, 4), (3, 4)]
</pre>
</pre>
See http://alahonua.com/img/temp.png for a png of the above plot. Note due to column major ordering both figure and contour
are rotated compared to the rows and columns of the data given.


=={{header|Lua}}==
=={{header|Lua}}==

Revision as of 07:37, 2 July 2022

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


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)))

  1. 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

  1. 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]
  1. The first (and pehaps only significant) level is the 0 <-> 1 transition border
  2. 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]

])

  1. 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

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