Adding a new type for storing point patterns

src/PointPatterns.jl

@@ -10,10 +10,14 @@ using Distributions
 
 import Random
 
+include("patterns.jl")
 include("processes.jl")
 include("thinning.jl")
 
 export
+  # patterns
+  PointPattern,
+
   # point processes
   PointProcess,
   BinomialProcess,

src/patterns.jl

@@ -0,0 +1,65 @@
+"""
+    PointPattern(points, region, marks)
+    PointPattern(points, region)
+
+Constructs a point pattern with the `points::PointSet` and `region::Geometry`. 
+If `marks` are provided, they must be an `AbstractVector` of the same length as points. 
+If `marks` are not provided or are `nothing`, the point pattern is unmarked.
+"""
+struct PointPattern{T,D,G,M}
+  points::PointSet{T,D}
+  region::G
+  marks::M
+  function PointPattern(points::PointSet{T,D}, region::G, marks::M) where {T,D,G<:Meshes.GeometryOrDomain,M}
+    check_marks(points, marks)
+    check_region(points, region)
+    new{T,D,G,M}(points, region, marks)
+  end
+  function PointPattern(points::PointSet{T,D}, region::G) where {T,D,G<:Meshes.GeometryOrDomain}
+    PointPattern(points, region, nothing)
+  end
+end
+
+"""
+    check_region(points, region)
+
+Checks that all points are contained in the region.
+"""
+function check_region(points::PointSet, region)
+  @assert embeddim(points) == embeddim(region)
+  @assert all(∈(region), points)
+  nothing
+end
+
+
+"""
+    check_marks(points::PointSet{T,D}, marks) where {T,D,M}
+
+Checks that the marks are `nothing` if the points are unmarked or an `AbstractVector` of the same length as points if marks are provided.
+"""
+check_marks(points::PointSet{T,D}, marks::Nothing) where {T,D} = nothing
+function check_marks(points::PointSet{T,D}, marks::AbstractVector{M}) where {T,D,M}
+  @assert length(points) == length(marks)
+  nothing
+end
+function Base.vcat(pp₁::PointPattern, pp₂::PointPattern)
+  pp₁.region == pp₂.region || throw(ArgumentError("Point patterns must have the same region."))
+  points = PointSet(vcat(collect(pp₁.points), collect(pp₂.points)))
+  marks = join_marks(pp₁.marks, pp₂.marks)
+  PointPattern(points, pp₁.region, marks)
+end
+
+"""
+    join_marks(marks₁, marks₂)
+
+Concatenate two collections of marks. Will error if exactly one of the processes is marked.
+"""
+join_marks(marks₁::Nothing, marks₂::Nothing) = nothing
+join_marks(marks₁::AbstractArray, marks₂::AbstractArray) = vcat(marks₁, marks₂)
+join_marks(::Nothing, ::AbstractArray) =
+  throw(ArgumentError("Cannot concatenate marked and unmarked point patterns."))
+join_marks(::AbstractArray, ::Nothing) =
+  throw(ArgumentError("Cannot concatenate marked and unmarked point patterns."))
+
+Meshes.nelements(pp::PointPattern) = nelements(pp.points)
+Base.length(pp::PointPattern) = nelements(pp)

src/processes/binomial.jl

@@ -13,4 +13,4 @@ end
 
 ishomogeneous(p::BinomialProcess) = true
 
-randsingle(rng::Random.AbstractRNG, p::BinomialProcess, g) = PointSet(sample(rng, g, HomogeneousSampling(p.n)))
+randsingle(rng::Random.AbstractRNG, p::BinomialProcess, g) = PointPattern(PointSet(sample(rng, g, HomogeneousSampling(p.n))), g)

src/processes/cluster.jl

@@ -31,11 +31,15 @@ function randsingle(rng::Random.AbstractRNG, p::ClusterProcess, g)
   parents = rand(rng, p.proc, g)
 
   # generate offsprings
-  offsprings = filter(!isnothing, p.ofun.(parents))
+  offsprings = filter(!isnothing, p.ofun.(parents.points))
 
   # intersect with geometry
-  intersects = filter(!isnothing, [o ∩ g for o in offsprings])
+  intersects = if eltype(offsprings) <: PointPattern # this is a temporary fix
+    filter(!isnothing, [o.points ∩ g for o in offsprings])
+  else
+    filter(!isnothing, [o ∩ g for o in offsprings])
+  end
 
   # combine offsprings into single set
-  isempty(intersects) ? nothing : PointSet(mapreduce(collect, vcat, intersects))
+  isempty(intersects) ? nothing : PointPattern(PointSet(mapreduce(collect, vcat, intersects)), g)
 end

src/processes/inhibition.jl

@@ -11,4 +11,4 @@ struct InhibitionProcess{D<:Real} <: PointProcess
   δ::D
 end
 
-randsingle(rng::Random.AbstractRNG, p::InhibitionProcess, g) = PointSet(sample(rng, g, MinDistanceSampling(p.δ)))
+randsingle(rng::Random.AbstractRNG, p::InhibitionProcess, g) = PointPattern(PointSet(sample(rng, g, MinDistanceSampling(p.δ))), g)

src/processes/poisson.jl

@@ -32,7 +32,7 @@ function randsingle(rng::Random.AbstractRNG, p::PoissonProcess{<:Real}, g)
   n = rand(rng, Poisson(λ * V))
 
   # simulate n points
-  iszero(n) ? nothing : PointSet(sample(rng, g, HomogeneousSampling(n)))
+  iszero(n) ? nothing : PointPattern(PointSet(sample(rng, g, HomogeneousSampling(n))), g)
 end
 
 #--------------------
@@ -56,7 +56,7 @@ function randsingle(rng::Random.AbstractRNG, p::PoissonProcess{<:AbstractVector}
   n = rand(rng, Poisson(sum(λ)))
 
   # simulate point pattern
-  iszero(n) ? nothing : PointSet(sample(rng, d, HomogeneousSampling(n, λ)))
+  iszero(n) ? nothing : PointPattern(PointSet(sample(rng, d, HomogeneousSampling(n, λ))), d)
 end
 
 function maxintensity(p::PoissonProcess{<:Function}, g)

src/processes/union.jl

@@ -17,7 +17,7 @@ ishomogeneous(p::UnionProcess) = ishomogeneous(p.p₁) && ishomogeneous(p.p₂)
 function randsingle(rng::Random.AbstractRNG, p::UnionProcess, g)
   pp₁ = rand(rng, p.p₁, g)
   pp₂ = rand(rng, p.p₂, g)
-  PointSet([collect(pp₁); collect(pp₂)])
+  vcat(pp₁, pp₂) # note this is unmarked as PointProcess is not marked currently.
 end
 
 # ----------

src/thinning/random.jl

@@ -36,21 +36,35 @@ thin(p::PoissonProcess{<:Function}, t::RandomThinning{<:Real}) = PoissonProcess(
 # THINNING POINT PATTERN
 # -----------------------
 
-function thin(pp::PointSet, t::RandomThinning{<:Real})
-  draws = rand(Bernoulli(t.p), nelements(pp))
-  inds = findall(isequal(1), draws)
+function thin(pp::PointPattern, t::ThinningMethod)
+  inds = thinning_index(pp.points, pp.marks, t)
+   points = collect(view(pp.points, inds))
+  marks = if pp.marks === nothing
+    nothing
+  else
+    collect(view(pp.marks, inds))
+  end
+  PointPattern(PointSet(points), pp.region, marks)
+end
+
+function thin(pp::PointSet, t::ThinningMethod)
+  inds = thinning_index(pp, nothing, t)
   points = collect(view(pp, inds))
   PointSet(points)
 end
 
-function thin(pp::PointSet, t::RandomThinning{<:Function})
+function thinning_index(pp::PointSet, marks::Nothing, t::RandomThinning{<:Real})
+  draws = rand(Bernoulli(t.p), nelements(pp))
+  findall(isequal(1), draws)
+end
+
+function thinning_index(pp::PointSet, marks::Nothing, t::RandomThinning{<:Function})
   inds = Vector{Int}()
   for j in 1:nelements(pp)
     x = element(pp, j)
     if rand(Bernoulli(t.p(x)))
       push!(inds, j)
     end
   end
-  points = collect(view(pp, inds))
-  PointSet(points)
-end
+  return inds
+end

test/patterns.jl

@@ -0,0 +1,57 @@
+import PointPatterns: check_marks, check_region
+@testset "Patterns" begin
+
+    @testset "check_region" begin
+        seg = Segment((0.0, 0.0), (11.3, 11.3))
+        tri = Triangle((0.0, 0.0), (5.65, 0.0), (5.65, 5.65))
+        quad = Quadrangle((0.0, 0.0), (0.0, 4.0), (4.0, 4.0), (4.0, 0.0))
+        box = Box((0.0, 0.0), (4.0, 4.0))
+        ball = Ball((1.0, 1.0), 2.25)
+        outer = [(0, -4), (4, -1), (4, 1.5), (0, 3)]
+        hole1 = [(0.2, -0.2), (1.4, -0.2), (1.4, 0.6), (0.2, 0.6)]
+        hole2 = [(2, -0.2), (3, -0.2), (3, 0.4), (2, 0.4)]
+        poly = PolyArea([outer, hole1, hole2])
+        regions = [seg, tri, quad, box, ball, poly]
+        points = PointSet([(100,100)])
+        for r in regions
+            @test_throws AssertionError check_region(points, r)
+        end
+    end
+
+    @testset "check_marks" begin
+        points = PointSet([(0,1), (0,3), (0,5), (0,7), (0,9)])
+        @test isnothing(check_marks(points, [1, 2, 3, 4, 5]))
+        @test isnothing(check_marks(points, nothing))
+        @test_throws AssertionError check_marks(points, [1, 2, 3, 4])
+    end
+
+    @testset "vcat" begin
+        region = Box((0,0), (1,10))
+        points₁ = PointSet([(0,1), (0,3), (0,5), (0,7), (0,9)])
+        marks₁ = [1, 2, 3, 4, 5]
+        points₂ = PointSet([(1,1), (1,3), (1,5), (1,7), (1,9)])
+        marks₂ = [6, 7, 8, 9, 10]
+
+        pp₁ = PointPattern(points₁, region, marks₁)
+        pp₂ = PointPattern(points₂, region, marks₂)
+        pp₃ = PointPattern(points₁, region, nothing)
+        pp₄ = PointPattern(points₂, region, nothing)
+        pp₅ = PointPattern(points₂, Box((0,0), (1,11)), marks₂)
+
+        pp₁₂ = vcat(pp₁, pp₂)
+        pp₃₄ = vcat(pp₃, pp₄)
+        @test pp₁₂.points == PointSet(vcat(parent(points₁), parent(points₂)))
+        @test pp₁₂.marks == vcat(marks₁, marks₂)
+        @test pp₁₂.region == pp₁.region
+        @test pp₃₄.points == PointSet(vcat(parent(points₁), parent(points₂)))
+        @test isnothing(pp₃₄.marks)
+        @test pp₃₄.region == pp₃.region
+
+        @test_throws ArgumentError vcat(pp₁, pp₃)
+        @test_throws ArgumentError vcat(pp₂, pp₄)
+        @test_throws ArgumentError vcat(pp₁, pp₅)
+
+        @test length(pp₁) == 5
+        @test PointPatterns.nelements(pp₁) == 5
+    end
+end

test/processes.jl

@@ -26,7 +26,8 @@
   @testset "Basic" begin
     for p in procs, g in geoms
       pp = rand(p, g)
-      @test all(∈(g), pp)
+      @test g == pp.region
+      @test all(∈(g), pp.points)
     end
   end
 
@@ -43,7 +44,7 @@
     for g in [grid, mesh]
       p = PoissonProcess(λ.(centroid.(g)))
       pp = rand(p, g)
-      @test all(∈(g), pp)
+      @test all(∈(g), pp.points)
     end
 
     # empty pointsets
@@ -78,7 +79,7 @@
       cp = ClusterProcess(p, ofun)
       pp = rand(cp, g)
       if !isnothing(pp)
-        @test all(∈(g), pp)
+        @test all(∈(g), pp.points)
       end
     end
   end
@@ -91,8 +92,8 @@
 
     s = rand(p, b, 2)
     @test length(s) == 2
-    @test s[1] isa PointSet
-    @test s[2] isa PointSet
+    @test s[1] isa PointPattern
+    @test s[2] isa PointPattern
     @test nelements.(s) == [100, 100]
   end
 end

test/runtests.jl

@@ -3,7 +3,7 @@ using Meshes
 using Test, Random
 
 # list of tests
-testfiles = ["processes.jl", "thinning.jl"]
+testfiles = ["patterns.jl" ,"processes.jl", "thinning.jl"]
 
 @testset "PointPatterns.jl" begin
   for testfile in testfiles

test/thinning.jl

@@ -5,14 +5,14 @@
     pp = rand(p, q)
     tp = thin(pp, RandomThinning(0.3))
     @test length(tp) ≤ length(pp)
-    xs = coordinates.(tp)
+    xs = coordinates.(tp.points)
     @test all(0 .≤ first.(xs) .≤ 4.0)
     @test all(0 .≤ last.(xs) .≤ 4.0)
 
     λ(s::Point2) = sum(coordinates(s) .^ 2)
     tp = thin(pp, RandomThinning(s -> λ(s) / λ(Point(4.0, 4.0))))
     @test length(tp) ≤ length(pp)
-    xs = coordinates.(tp)
+    xs = coordinates.(tp.points)
     @test all(0 .≤ first.(xs) .≤ 4.0)
     @test all(0 .≤ last.(xs) .≤ 4.0)
   end