commit 3a15d9e39ecaeda0fe5433300bca2290d9af51a9
parent e5fe4e54a81b423af217a9eaa8ddd442a6f3ba4a
Author: Anders Damsgaard <andersd@riseup.net>
Date: Fri, 22 Sep 2017 13:23:55 -0400
begin implementing poisson-disk packing
Diffstat:
| M | src/grid.jl | | | 72 | +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++- |
| M | src/packing.jl | | | 102 | +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ |
| M | test/runtests.jl | | | 1 | + |
3 files changed, 174 insertions(+), 1 deletion(-)
diff --git a/src/grid.jl b/src/grid.jl
@@ -326,11 +326,81 @@ function isPointInCell(grid::Any, i::Int, j::Int, point::Vector{Float64},
end
end
+export isPointInGrid
+"""
+Check if a 2d point is contained inside the grid. The function uses either an
+area-based approach (`method = "Area"`), or a conformal mapping approach
+(`method = "Conformal"`). The area-based approach is more robust. This
+function returns `true` or `false`.
+"""
+function isPointInGrid(grid::Any, point::Vector{Float64},
+ sw::Vector{Float64} = Vector{Float64}(2),
+ se::Vector{Float64} = Vector{Float64}(2),
+ ne::Vector{Float64} = Vector{Float64}(2),
+ nw::Vector{Float64} = Vector{Float64}(2);
+ method::String="Conformal")
+
+ #sw, se, ne, nw = getCellCornerCoordinates(grid.xq, grid.yq, i, j)
+ nx, ny = size(grid.xq)
+ @views sw[1] = grid.xq[ 1, 1]
+ @views sw[2] = grid.yq[ 1, 1]
+ @views se[1] = grid.xq[ nx, 1]
+ @views se[2] = grid.yq[ nx, 1]
+ @views ne[1] = grid.xq[ nx, ny]
+ @views ne[2] = grid.yq[ nx, ny]
+ @views nw[1] = grid.xq[ 1, ny]
+ @views nw[2] = grid.yq[ 1, ny]
+
+ if method == "Area"
+ if areaOfQuadrilateral(sw, se, ne, nw) ≈
+ areaOfTriangle(point, sw, se) +
+ areaOfTriangle(point, se, ne) +
+ areaOfTriangle(point, ne, nw) +
+ areaOfTriangle(point, nw, sw)
+ return true
+ else
+ return false
+ end
+
+ elseif method == "Conformal"
+ x_tilde, y_tilde = conformalQuadrilateralCoordinates(sw, se, ne, nw,
+ point)
+ if x_tilde >= 0. && x_tilde <= 1. && y_tilde >= 0. && y_tilde <= 1.
+ return true
+ else
+ return false
+ end
+ else
+ error("method not understood")
+ end
+end
+
+export getGridCornerCoordinates
+"""
+ getGridCornerCoordinates(xq, yq)
+
+Returns grid corner coordinates in the following order (south-west corner,
+south-east corner, north-east corner, north-west corner).
+
+# Arguments
+* `xq::Array{Float64, 2}`: nominal longitude of q-points [degrees_E]
+* `yq::Array{Float64, 2}`: nominal latitude of q-points [degrees_N]
+"""
+@inline function getGridCornerCoordinates(xq::Array{Float64, 2},
+ yq::Array{Float64, 2})
+ nx, ny = size(xq)
+ @inbounds return Float64[xq[ 1, 1], yq[ 1, 1]],
+ Float64[xq[ nx, 1], yq[ nx, 1]],
+ Float64[xq[ nx, ny], yq[ nx, ny]],
+ Float64[xq[ 1, ny], yq[ 1, ny]]
+end
+
+
export getCellCornerCoordinates
"""
getCellCornerCoordinates(xq, yq, i, j)
-Returns grid corner coordinates in the following order (south-west corner,
+Returns grid-cell corner coordinates in the following order (south-west corner,
south-east corner, north-east corner, north-west corner).
# Arguments
diff --git a/src/packing.jl b/src/packing.jl
@@ -0,0 +1,102 @@
+## Functions for creating ice-floe packings
+"""
+Return random point in spherical annulus between `(r_i + r_j)` and `2.*(r_i +
+r_j)` around `x_i`. Note: there is slightly higher point density towards (r_i +
+r_j).
+"""
+function generateNeighboringPoint(x_i::Vector, r_i::Real, r_j::Real,
+ max_padding_factor::Real)
+
+ R = rand()*(r_i + r_j)*max_padding_factor + 2.*(r_i + r_j)
+ T = rand()*2.*pi
+ return [x_i[1] + R*sin(T), x_i[2] + R*cos(T)]
+end
+
+export poissonDiscSampling
+"""
+Generate disc packing in 2D using Poisson disc sampling with O(N) complexity, as
+described by [Robert Bridson (2007)](http://www.cs.ubc.ca/~rbridson/docs/bridson-siggraph07-poissondisk.pdf).
+
+# Arguments
+* `simulation::Simulation`: simulation object where ice floes are inserted.
+* `radius_max::Real`: largest ice-floe radius to use.
+* `radius_min::Real`: smallest ice-floe radius to use.
+* `sample_limit::Integer=30`: number of points to sample around each ice floe
+ before giving up.
+* `max_padding_factor::Real=2.`: this factor scales the padding to use during ice
+ floe generation in numbers of grain diameters.
+* `verbose::Bool=true`: show diagnostic information to stdout.
+
+"""
+function poissonDiscSampling(simulation::Simulation;
+ ;radius_max::Real=.1,
+ radius_min::Real=.1,
+ sample_limit::Integer=30,
+ max_padding_factor::Real=2.,
+ thickness::Real=1.,
+ verbose::Bool=true)
+
+ active_list = Int[] # list of points to originate search from
+
+ # Step 0: Use existing `grid` (ocean or atmosphere) for contact-search grid
+ if typeof(simulation.ocean.input_file) != Bool
+ grid = simulation.ocean
+ elseif typeof(simulation.atmosphere.input_file) != Bool
+ grid = simulation.atmosphere
+ else
+ error("poissonDiscSampling requires an ocean or atmosphere grid")
+ end
+ # save grid boundaries
+ sw, se, ne, nw = getGridCornerCoordinates(grid.xq, grid.yq)
+ width_x = se[1] - sw[1] # assume regular grid
+ width_y = nw[2] - sw[2] # assume regular grid
+
+ # Step 1: If grid is empty: select random initial sample and save its index
+ # to the background grid. Else: Make all ice floes active for search
+ if isempty(simulation.ice_floes)
+ r = rand()*(radius_max - radius_min) + radius_min
+ x0 = rand(2).*[width_x, width_y] + sw
+ addIceFloeCylindrical!(simulation, x0, r, thickness)
+ sortIceFloesInGrid!(simulation, grid)
+ push!(active_list, 1)
+ else
+ for i=1:length(simulation.ice_floes)
+ push!(active_list, i)
+ end
+ end
+
+ # Step 2: While the active list is not empty, choose a random index `i` from
+ # it. Generate up to `k` points chosen uniformly from the spherical annulus
+ # between radius `2*(r_i+r_j)` and `max_padding_factor*2*(r_i+r_j)` around
+ # `x_i`.
+ i = 0; j = 0; x_i = zeros(2); x_j = zeros(2); r_i = 0.; r_j = 0.
+ while !isempty(active_list)
+
+ i = rand(1:length(active_list))
+ deleteat!(active_list, i)
+
+ x_i = simulation.ice_floes[i].lin_pos
+ r_i = simulation.ice_floes[i].contact_radius
+ r_j = rand()*(radius_max - radius_min) + radius_min
+
+ for j=1:sample_limit
+
+ x_j = generateNeighboringPoint(x_i, r_i, r_j, max_padding_factor)
+ if !(isPointInGrid(grid, x_j))
+ continue
+ end
+
+ # if it doesn't collide, add it to the active list
+ if !checkCollisions(grid, x_j, dx, position_list, radius_list)
+ push!(position_list, x_j)
+ push!(radius_list, r_j)
+ push!(active_list, length(radius_list))
+ end
+ end
+ println("Active points: $(length(active_list))")
+ end
+ if verbose
+ info("Generated $(length(radius_list)) points")
+ end
+ return position_list, radius_list
+end
diff --git a/test/runtests.jl b/test/runtests.jl
@@ -2,6 +2,7 @@ import SeaIce
using Base.Test
include("icefloe.jl")
+include("packing.jl")
include("util.jl")
include("temporal.jl")
include("contact-search-and-geometry.jl")
