diff --git a/Project.toml b/Project.toml index f11312c..93beaee 100644 --- a/Project.toml +++ b/Project.toml @@ -1,9 +1,10 @@ name = "NeutralLandscapes" uuid = "71847384-8354-4223-ac08-659a5128069f" -authors = ["Timothée Poisot ", "Michael Krabbe Borregaard "] +authors = ["Timothée Poisot ", "Michael Krabbe Borregaard ", "Michael David Catchen "] 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" diff --git a/docs/src/gallery.md b/docs/src/gallery.md index 3d38fff..3215848 100644 --- a/docs/src/gallery.md +++ b/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()) +``` diff --git a/src/NeutralLandscapes.jl b/src/NeutralLandscapes.jl index 9dbbb18..0899eb6 100644 --- a/src/NeutralLandscapes.jl +++ b/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, MidpointDisplacement + include(joinpath("algorithms", "wavesurface.jl")) export WaveSurface diff --git a/src/algorithms/diamondsquare.jl b/src/algorithms/diamondsquare.jl new file mode 100644 index 0000000..4776cf0 --- /dev/null +++ b/src/algorithms/diamondsquare.jl @@ -0,0 +1,269 @@ +""" + DiamondSquare + +This type generates a neutral landscape using the diamond-squares +algorithm, which produces fractals with variable spatial autocorrelation. + +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 behave as expected. +The result of the diamond-square algorithm is a fractal with dimension D = 2 + H. + +A similar algorithm, midpoint-displacement, functions almost +identically, except that in DiamondSquare, the square step interpolates +edge midpoints from the nearest two corners and the square's center, where as +midpoint-displacement only intepolates from the nearest corners (see `MidpointDisplacement`). + +""" +struct DiamondSquare <: NeutralLandscapeMaker + H::Float64 + function DiamondSquare(H::T) where {T <: Real} + @assert 0 <= H < 1 + new(H) + end +end + +""" + DiamondSquare() + + Initializes DiamondSqaure with H = 0.5. +""" +DiamondSquare() = DiamondSquare(0.5) + +""" + MidpointDisplacement() + + Creates a midpoint-displacement algorithm object `MidpointDisplacement`. 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. + + A similar algorithm, diamond-square, functions almost + identically, except that in diamond-square, the square step interpolates + edge midpoints from the nearest two corners and the square's center, where as + `MidpointDisplacement` only intepolates from the nearest corners (see `DiamondSquare`). +""" +struct MidpointDisplacement <: NeutralLandscapeMaker + H::Float64 + function MidpointDisplacement(H::T) where {T <: Real} + @assert 0 <= H < 1 + new(H) + end +end +MidpointDisplacement() = MidpointDisplacement(0.5) + +""" + _landscape!(mat, alg::Union{DiamondSquare, MidpointDisplacement}; kw...) + + Check if `mat` is the right size. + 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. +""" +function _landscape!(mat, alg::Union{DiamondSquare, MidpointDisplacement}; 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 = 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.*i for i in [(x,y), (x+1, y), (x, y+1), (x+1, 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::MidpointDisplacement, round::Int, corners::AbstractVector{Tuple{Int,Int}}) + + Runs the square step of the `MidpointDisplacement` 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::MidpointDisplacement, 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(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 x & (x-1) == 0 +end