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.

Task

Generate contours for a two-dimensional scalar field.

See: Marching squares

Julia

Uses the marching squares algorithm: see github.com/JuliaGeometry/Contour.jl/blob/master/src/Contour.jl <lang ruby>using Contour

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;

])

cl = first(levels(contours(1:6, 1:5, example))) xs, ys = coordinates(first(lines(cl)))

1. 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)] @show points

1. Plotting

using GLMakie plt = plot(example) lines!(xs, ys, linewidth = 8, color = :red) display(plt)

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

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.

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]