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

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

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]