[feature] diamond-square and MPD

Project.toml

@@ -1,9 +1,10 @@
 name = "NeutralLandscapes"
 uuid = "71847384-8354-4223-ac08-659a5128069f"
-authors = ["Timothée Poisot <[email protected]>", "Michael Krabbe Borregaard <[email protected]>"]
+authors = ["Timothée Poisot <[email protected]>", "Michael Krabbe Borregaard <[email protected]>", "Michael David Catchen <[email protected]>"]
 version = "0.0.1"
 
 [deps]
+Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
 NaNMath = "77ba4419-2d1f-58cd-9bb1-8ffee604a2e3"
 NearestNeighbors = "b8a86587-4115-5ab1-83bc-aa920d37bbce"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"

docs/src/gallery.md

@@ -32,7 +32,6 @@ demolandscape(WaveSurface(35, 3))
 ```
 
 ## Rectangular cluster
-
 ```@example gallery
 demolandscape(RectangularCluster())
 ```
@@ -56,3 +55,13 @@ demolandscape(PerlinNoise())
 sources = unique(rand(1:40000, 50))
 heatmap(NeutralLandscapes.classify!(rand(DistanceGradient(sources), (200, 200)), [0.5, 1, 1, 0.5]))
 ```
+
+## Diamond Square
+```@example gallery
+demolandscape(DiamondSquare())
+```
+
+## Midpoint Displacement
+```@example gallery
+demolandscape(MPD())
+```

src/NeutralLandscapes.jl

@@ -2,7 +2,8 @@ module NeutralLandscapes
 
 import NaNMath
 using Random: rand!
-using Statistics: quantile
+using Distributions: Normal
+using Statistics: quantile, mean
 using NearestNeighbors: KDTree, nn
 
 abstract type NeutralLandscapeMaker end
@@ -20,6 +21,9 @@ export PlanarGradient
 include(joinpath("algorithms", "edgegradient.jl"))
 export EdgeGradient
 
+include(joinpath("algorithms", "diamondsquare.jl"))
+export DiamondSquare, MPD
+
 include(joinpath("algorithms", "wavesurface.jl"))
 export WaveSurface
 

src/algorithms/diamondsquare.jl

@@ -0,0 +1,258 @@
+"""
+    DiamondSquare
+
+This type generates a neutral landscape using the diamond-sqaures
+algorithm, which produces fractals with variable spatialautocorrelation.
+
+https://en.wikipedia.org/wiki/Diamond-square_algorithm
+
+The algorithm is named diamond-square because it is an iterative procedure of
+"diamond" and "square" steps.
+
+
+
+The degree of spatial autocorrelation is controlled by a parameter `H`,
+which varies from 0.0 (low autocorrelation) to 1.0 (high autocorrelation) --- note this is non-inclusive and H = 0 and H = 1 will not behavive as expected.
+The result of the diamond-square algorithm is a fractal with dimension D = 2 + H.
+
+A similar algorithm, midpoint displacement *TODO link to mpd in docs*, almost
+identically, except that in DiamondSquare, the ed
+"""
+struct DiamondSquare <: NeutralLandscapeMaker
+    H::Float64
+    function DiamondSquare(H::T) where {T <: Real}
+        @assert 0 < H && H < 1
+        new(H)
+    end
+    DiamondSquare() = new(0.5)
+
+end
+
+"""
+    MPD()
+
+    Creates a midpoint-displacement algorithm object `MPD`. The degree of spatial autocorrelation is controlled by a parameter `H`,
+    which varies from 0.0 (low autocorrelation) to 1.0 (high autocorrelation) --- note this is non-inclusive and H = 0 and H = 1 will not behavive as expected.
+"""
+struct MPD <: NeutralLandscapeMaker
+    H::Float64
+    function MPD(H::T) where {T <: Real}
+        @assert 0 < H && H < 1
+        new(H)
+    end
+    MPD() = new(0.5)
+end
+
+
+"""
+    _landscape!(mat, alg::Union{DiamondSquare, MPD}; kw...)
+
+    Check if `mat` is the right size and already initialized.
+    If mat is not the correct size (DiamondSquare can only run on a lattice of size NxN where N = (2^n)+1 for integer n),
+    allocates the smallest lattice large enough to contain `mat` that can run DiamondSquare.
+    Will initialize `mat` to all zeros before running DiamondSquare if it is not initialized already.
+"""
+function _landscape!(mat, alg::Union{DiamondSquare, MPD}; kw...) where {IT <: Integer}
+
+    rightSize::Bool = _isPowerOfTwo(size(mat)[1]-1) && _isPowerOfTwo(size(mat)[2]-1)
+    latticeSize::Int = size(mat)[1]
+
+    dsMat = mat
+    if !rightSize
+        dim1, dim2 = size(mat)
+        smallestContainingLattice::Int = 2^ceil(log2(max(dim1, dim2))) + 1
+        @warn "$alg cannot be run on the input dimensions ($dim1 x $dim2),
+                and will instead run on the next smallest valid size ($smallestContainingLattice x $smallestContainingLattice).
+                This can slow performance as it involves additional memory allocation."
+        dsMat = zeros(smallestContainingLattice, smallestContainingLattice)
+    end
+    _diamondsquare!(dsMat, alg)
+
+    mat .= dsMat[1:size(mat)[1], 1:size(mat)[2]]
+end
+
+"""
+    _diamondsquare!(mat, alg)
+
+    Runs the diamond-square algorithm on a matrix `mat` of size
+    `NxN`, where `N=(2^n)+1` for some integer `n`, i.e (N=5,9,17,33,65)
+
+    Diamond-square is an iterative procedure, where the lattice is divided
+    into subsquares in subsequent rounds. At each round, the subsquares shrink in size,
+    as previously uninitialized values in the lattice are interpolated as a mean of nearby points plus random displacement.
+    As the rounds increase, the magnitude of this displacement decreases. This creates spatioautocorrelation, which is controlled
+    by a single parameter `H` which varies between `0` (no autocorrelation) and `1` (high autocorrelation)
+
+"""
+function _diamondsquare!(mat, alg)
+    latticeSize::Int = size(mat)[1]
+    numberOfRounds::Int = log2(latticeSize-1)
+    _initializeDiamondSquare!(mat, alg)
+
+    for round in 0:(numberOfRounds-1)    # counting from 0 saves us a headache later
+        subsquareSideLength::Int = 2^(numberOfRounds-(round))
+        numberOfSubsquaresPerAxis::Int = ((latticeSize-1) / subsquareSideLength)-1
+
+        for x in 0:numberOfSubsquaresPerAxis         # iterate over the subsquares within the lattice at this side length
+            for y in 0:numberOfSubsquaresPerAxis
+                subsquareCorners = _subsquareCornerCoordinates(x,y,subsquareSideLength)
+
+                _diamond!(mat, alg, round, subsquareCorners)
+                _square!(mat, alg, round, subsquareCorners)
+            end
+        end
+    end
+end
+
+"""
+    _initializeDiamondSquare!(mat, alg)
+
+    Initialize's the `DiamondSquare` algorithm by displacing the four corners of the
+    lattice using `displace`, scaled by the algorithm's autocorrelation `H`.
+"""
+function _initializeDiamondSquare!(mat, alg)
+    latticeSize = size(mat)[1]
+    corners = _subsquareCornerCoordinates(0,0, latticeSize-1)
+    for mp in corners
+        mat[mp...] = _displace(alg.H, 1)
+    end
+end
+
+"""
+    _subsquareCornerCoordinates(x::Int, y::Int, sideLength::Int)
+
+    Returns the coordinates for the corners of the subsquare (x,y) given a side-length `sideLength`.
+"""
+function _subsquareCornerCoordinates(x::Int, y::Int, sideLength::Int)
+    corners = [(1+sideLength*x, 1+sideLength*y), (1+sideLength*(x+1), 1+sideLength*y), (1+sideLength*x, sideLength*(y+1)+1), (1+sideLength*(x+1), 1+sideLength*(y+1))]
+end
+
+"""
+     _diamond!(mat, alg::DiamondSquare, round::Int, corners::AbstractVector{Tuple{Int, Int}})
+
+     Runs the diamond step of the `DiamondSquare` algorithm on the square defined by
+     `corners` on the matrix `mat`. The center of the square is interpolated from the
+     four corners, and is displaced. The displacement is drawn according to `alg.H` and round using `displace`
+"""
+function _diamond!(mat, alg, round::Int, corners::AbstractVector{Tuple{Int, Int}})
+    centerPt = _centerCoordinate(corners)
+    mat[centerPt...] = _interpolate(mat, corners) + _displace(alg.H, round)
+end
+
+"""
+    _square!(mat, alg::DiamondSquare, round::Int, corners::AbstractVector{Tuple{Int,Int}})
+
+    Runs the square step of the `DiamondSquare` algorithm on the square defined
+    by `corners` on the matrix `mat`. The midpoint of each edge of this square is interpolated
+    by computing the mean value of the two corners on the edge and the center of the square, and the
+    displacing it. The displacement is drawn according to `alg.H` and round using `displace`
+
+"""
+function _square!(mat, alg::DiamondSquare, round::Int, corners::AbstractVector{Tuple{Int, Int}})
+    bottomLeft,bottomRight,topLeft,topRight = corners
+    leftEdge, bottomEdge, topEdge, rightEdge = _edgeMidpointCoordinates(corners)
+    centerPoint = _centerCoordinate(corners)
+
+    mat[leftEdge...] = _interpolate(mat, [topLeft,bottomLeft,centerPoint]) + _displace(alg.H, round)
+    mat[bottomEdge...] = _interpolate(mat, [bottomLeft,bottomRight,centerPoint]) + _displace(alg.H, round)
+    mat[topEdge...] = _interpolate(mat, [topLeft,topRight,centerPoint]) + _displace(alg.H, round)
+    mat[rightEdge...] = _interpolate(mat, [topRight,bottomRight,centerPoint]) + _displace(alg.H, round)
+end
+
+"""
+    _square!(mat, alg::MPD, round::Int, corners::AbstractVector{Tuple{Int,Int}})
+
+    Runs the square step of the `MPD` algorithm on the square defined
+    by `corners` on the matrix `mat`. The midpoint of each edge of this square is interpolated
+    by computing the mean value of the two corners on the edge and the center of the square, and the
+    displacing it. The displacement is drawn according to `alg.H` and round using `displace`
+"""
+function _square!(mat, alg::MPD, round::Int, corners::AbstractVector{Tuple{Int, Int}})
+    bottomLeft,bottomRight,topLeft,topRight = corners
+    leftEdge, bottomEdge, topEdge, rightEdge = _edgeMidpointCoordinates(corners)
+    mat[leftEdge...] = _interpolate(mat, [topLeft,bottomLeft]) + _displace(alg.H, round)
+    mat[bottomEdge...] = _interpolate(mat, [bottomLeft,bottomRight]) + _displace(alg.H, round)
+    mat[topEdge...] = _interpolate(mat, [topLeft,topRight]) + _displace(alg.H, round)
+    mat[rightEdge...] = _interpolate(mat, [topRight,bottomRight]) + _displace(alg.H, round)
+end
+
+"""
+    _interpolate(mat, points::AbstractVector{Tuple{Int,Int}})
+
+    Computes the mean of a set of points, represented as a list of indecies to a matrix `mat`.
+"""
+function _interpolate(mat, points::AbstractVector{Tuple{Int,Int}})
+    return mean(collect(mat[pt...] for pt in points))
+end
+
+"""
+    _displace(H::Float64, round::Int)
+
+    `displace` produces a random value as a function of  `H`, which is the
+    autocorrelation parameter used in `DiamondSquare` and must be between `0`
+    and `1`, and `round` which describes the current tiling size for the
+    DiamondSquare() algorithm.
+
+    Random value are drawn from a Gaussian distribution using `Distribution.Normal`
+    The standard deviation of this Gaussian, σ, is set to (1/2)^(round*H), which will
+    move from 1.0 to 0 as `round` increases.
+
+"""
+function _displace(H::Float64, round::Int)
+    σ = (0.5)^(round*H)
+    return(rand(Normal(0, σ)))
+end
+
+"""
+    _centerCoordinate(corners::AbstractVector{Tuple{Int,Int}})
+
+    Returns the center coordinate for a square defined by `corners` for the
+    `DiamondSquare` algorithm.
+"""
+function _centerCoordinate(corners::AbstractVector{Tuple{Int,Int}})
+    bottomLeft,bottomRight,topLeft,topRight = corners
+    centerX::Int = (_xcoord(bottomLeft)+ _xcoord(bottomRight))/2
+    centerY::Int = (_ycoord(topRight)+ _ycoord(bottomRight))/2
+    return (centerX, centerY)
+end
+
+"""
+    _xcoord(pt::Tuple{Int,Int})
+
+    Returns the x-coordinate from a lattice coordinate `pt`.
+"""
+ _xcoord(pt::Tuple{Int,Int}) = pt[1]
+
+ """
+     _ycoord(pt::Tuple{Int,Int})
+
+     Returns the y-coordinate from a lattice coordinate `pt`.
+ """
+ _ycoord(pt::Tuple{Int,Int}) = pt[2]
+
+"""
+    _edgeMidpointCoordinates(corners::AbstractVector{Tuple{Int,Int}})
+
+    Returns an array of midpoints for a square defined by `corners` for the `DiamondSquare` algorithm.
+"""
+function _edgeMidpointCoordinates(corners::AbstractVector{Tuple{Int,Int}})
+    # bottom left, bottom right, top left, top right
+    bottomLeft,bottomRight,topLeft,topRight = corners
+
+    leftEdgeMidpoint::Tuple{Int,Int} = (_xcoord(bottomLeft), ( _ycoord(bottomLeft)+_ycoord(topLeft) )/2 )
+    bottomEdgeMidpoint::Tuple{Int,Int} = ( (_xcoord(bottomLeft)+ _xcoord(bottomRight))/2, _ycoord(bottomLeft) )
+    topEdgeMidpoint::Tuple{Int,Int} = ( (_xcoord(topLeft)+_xcoord(topRight))/2, _ycoord(topLeft))
+    rightEdgeMidpoint::Tuple{Int,Int} = ( _xcoord(bottomRight), (_ycoord(bottomRight)+_ycoord(topRight))/2)
+
+    edgeMidpoints = [leftEdgeMidpoint, bottomEdgeMidpoint, topEdgeMidpoint, rightEdgeMidpoint]
+    return edgeMidpoints
+end
+
+"""
+     _isPowerOfTwo(x::IT) where {IT <: Integer}
+
+     Determines if `x`, an integer, can be expressed as `2^n`, where `n` is also an integer.
+"""
+function _isPowerOfTwo(x::IT) where {IT <: Integer}
+    return log2(x) == floor(log2(x))
+end