Initial commit.

Author: JakeGrainger

.gitignore

@@ -0,0 +1,2 @@
+.vscode
+.DS_store

Fig1_introexample.jl

@@ -0,0 +1,89 @@
+using PointProcessFilters, CairoMakie, Meshes, ProgressMeter, RCall, LinearAlgebra
+
+## Simulation
+@info "Perform the simulation..."
+side_length = 200
+side_padding = 50
+box_x = [-side_padding, side_length + side_padding]
+box_y = box_x
+
+λ = 0.01
+kappa = 0.8λ
+scale = 3
+mu = λ / kappa
+
+@rput kappa scale mu box_x box_y
+R"""
+    library(spatstat)
+    set.seed(1234)
+    simulated_process = rThomas(kappa, scale, mu, win=owin(box_x, box_y))
+"""
+@rget simulated_process
+simulated_points = PointSet(Meshes.Point.(simulated_process[:x], simulated_process[:y]))
+@info "points simulated."
+
+## Filtering
+@info "Apply the filters..."
+interval = [0, 0.05, 0.1, 0.15]
+filtered_process = [
+    centeredcirclepass(simulated_points, interval[2]),
+    ringpass(simulated_points, interval[2], interval[3]),
+    ringpass(simulated_points, interval[3], interval[4]),
+]
+filter_resolution = 100
+filter_x = range(0, side_length, filter_resolution)
+filter_y = filter_x
+Y = @showprogress "evaluating filters... " [y.(filter_x, filter_y')
+                                            for y in filtered_process]
+Ylims = extrema(stack(Y))
+@info "filtered processes computed."
+
+## sdf of thomas process
+function thomassdf(k; kappa, scale, mu)
+    return kappa * mu * (1 + mu * exp(-scale^2 * (2π)^2 * norm(k)^2))
+end
+
+## plotting
+@info "Start plotting..."
+fig_res = (80, 80)
+
+## plot the spectra
+@info "- plotting spectra"
+figure = Figure(resolution = fig_res .* (1, 2 / 3), figure_padding = 1)
+ax = Axis(figure[1, 1], yscale = identity, xticks = [0], yticks = [0], xgridcolor = :black,
+          ygridcolor = :black)
+band!(ax, interval, zeros(length(interval)), fill(0.025, length(interval)),
+      color = (:lightgrey, 0.5))
+vlines!(interval, color = :black, linestyle = :dash, linewidth = 1 / 2)
+lines!(ax, range(0, 0.2, 1000), x -> thomassdf(x, kappa = kappa, scale = scale, mu = mu),
+       color = :black)
+
+ylims!(ax, (-0.001, 0.025))
+hidedecorations!(ax, grid = false)
+hidespines!(ax)
+resize_to_layout!(figure)
+
+save("fig/fig1/sdf.pdf", figure)
+figure
+
+## plot the points
+@info "- plotting points"
+figure = Figure(resolution = fig_res, figure_padding = 1)
+ax = Axis(figure[1, 1], limits = (0, side_length, 0, side_length))
+viz!(ax, simulated_points, color = :black, pointsize = 2)
+hidedecorations!(ax)
+resize_to_layout!(figure)
+save("fig/fig1/points.pdf", figure)
+figure
+
+## plot filters
+@info "- plotting filters"
+for i in eachindex(Y)
+    figure_filt = Figure(resolution = fig_res, figure_padding = 1)
+    ax_filt = Axis(figure_filt[1, 1])
+    contourf!(ax_filt, filter_x, filter_y, Y[i], colorrange = Ylims, colormap = :tempo)
+    hidedecorations!(ax_filt)
+    save("fig/fig1/filtered_$i.pdf", figure_filt)
+end
+
+@info "Finished plotting."

Fig2_lgcp.jl

@@ -0,0 +1,106 @@
+using CairoMakie, PointProcessFilters, ProgressMeter, GaussianRandomFields, RCall, Meshes
+import Meshes: Point, PointSet, Box
+include("src/iso_pgram.jl")
+import .IsoPP: periodogram_iso
+
+## set parameter
+side_length = 100
+side_padding = 25
+box_x = [-side_padding, side_length + side_padding]
+box_y = box_x
+nu = 3 / 2
+lengthscale = 10
+sigma = sqrt(2)
+lambda_noise = 0.05
+
+## simualate the random field (R package currently broken so cant use spatstat)
+import Random;
+Random.seed!(1234);
+cov = CovarianceFunction(2, Matern(lengthscale, nu, σ = sigma))
+grid_sides = range(box_x..., 250), range(box_y..., 250)
+grf = GaussianRandomField(cov, CirculantEmbedding(), grid_sides...)
+intensity = exp.(GaussianRandomFields.sample(grf) .- 4)
+heatmap(intensity)
+
+##
+@rput intensity box_x box_y lambda_noise
+R"""
+    library(spatstat)
+    set.seed(1234)
+    window = owin(box_x, box_y)
+    intensity_image = as.im(aperm(intensity), window)
+    signal_process = rpoispp(intensity_image, nsim=2)
+    noise_process = rpoispp(lambda_noise, win=window, nsim=2)
+"""
+@rget signal_process noise_process
+sim_names = collect(keys(signal_process))
+signal = Dict([k => PointSet(Point.(signal_process[k][:x], signal_process[k][:y]))
+               for k in sim_names])
+noisy = Dict([k => PointSet([Point.(signal_process[k][:x], signal_process[k][:y]);
+                             Point.(noise_process[k][:x], noise_process[k][:y])])
+              for k in sim_names])
+
+## apply filters
+@info "Filtering (may take a few minutes)..."
+area = diff(box_x)[1] * diff(box_y)[1]
+filtered = Dict([k => centeredcirclepass(noisy[k], 0.05).(grid_sides[1], grid_sides[2]') .-
+                      length(noisy[k]) / area for k in sim_names])
+@info "filtering completed."
+
+##
+clims = extrema(vcat(values(filtered)...))
+f = Figure()
+heatmap(f[1, 1], grid_sides..., intensity, colormap = :matter)
+for i in 1:2
+    viz(f[2, i], signal[sim_names[i]], pointsize = 4, color = :black,
+        axis = (limits = (0, 100, 0, 100),))
+    viz(f[3, i], noisy[sim_names[i]], pointsize = 4, color = :black,
+        axis = (limits = (0, 100, 0, 100),))
+    heatmap(f[4, i], grid_sides..., filtered[sim_names[i]], colormap = :tempo,
+            colorrange = clims)
+end
+f
+
+region = Box(Point(box_x[1],box_y[1]), Point(box_x[2],box_y[2]))
+pgram = Dict([s=>periodogram_iso(noisy[s], region, (1:20).*0.01) for s in sim_names])
+λ = Dict([s=>length(noisy[s]) / area for s in sim_names])
+
+##
+@info "Making figures..."
+clims = extrema(vcat(values(filtered)...))
+function make_figure()
+    f = Figure(resolution = (100, 100), figure_padding = 1)
+    ax = Axis(f[1, 1], limits = (0, 100, 0, 100), aspect = 1)
+    hidedecorations!(ax)
+    return f, ax
+end
+
+figure, ax = make_figure()
+contourf!(ax, grid_sides..., intensity, colormap = :matter)
+save("fig/fig2/intensity.pdf", figure)
+
+for i in 1:2
+    figure_p, ax_p = make_figure()
+    viz!(ax_p, signal[sim_names[i]], pointsize = 2, color = :black)
+    save("fig/fig2/cp$i.pdf", figure_p)
+
+    figure_p, ax_p = make_figure()
+    viz!(ax_p, noisy[sim_names[i]], pointsize = 2, color = :black)
+    save("fig/fig2/pts$i.pdf", figure_p)
+
+    figure_p, ax_p = make_figure()
+    contourf!(ax_p, grid_sides..., filtered[sim_names[i]], colormap = :tempo,
+              colorrange = clims)
+    save("fig/fig2/filt$i.pdf", figure_p)
+
+    figure_p = Figure(resolution = (100,60), fontsize=8, figure_padding=1)
+    ax_p = Axis(figure_p[1,1], xticks = [0], yticks=[0], xgridcolor=:black, ygridcolor=:black)
+    ylims!(ax_p, (-1,21))
+    band!(ax_p, [0,0.05], fill(-1,2), fill(21,2), color = (:lightgrey,0.5))
+    vlines!(ax_p, [0,0.05], color = :black, linestyle=:dash, linewidth=1/2)
+    lines!(ax_p, pgram[sim_names[i]].t, pgram[sim_names[i]].sdf./pgram[sim_names[i]].lambda.-1, color = :black)
+    hidedecorations!(ax_p, grid=false)
+    hidespines!(ax_p)
+    save("fig/fig2/sdf_$i.pdf", figure_p)
+end
+@info "figures finished."

Fig3_lansingwoods.jl

@@ -0,0 +1,55 @@
+using PointProcessFilters, CairoMakie, Meshes, ProgressMeter
+import Meshes: Point, PointSet, Box
+using RCall
+include("src/iso_pgram.jl")
+import .IsoPP: periodogram_iso
+
+@info "Begin analysis of spatstat data"
+
+R"""
+    library(spatstat)
+    treedata = lansing
+"""
+@rget treedata
+
+lansing_points = Dict([s=>PointSet(unique(Point.(treedata[:x],treedata[:y])[treedata[:marks].==s])) for s in unique(treedata[:marks])])
+region = Box(Point(0,0), Point(1,1))
+
+species = ["hickory", "maple"]
+λ = Dict([s=>length(lansing_points[s]) / measure(region) for s in species])
+ppfilt = Dict([s=>centeredcirclepass(lansing_points[s], 6) for s in species])
+
+lags = range(0,1,100)
+pp_w = Dict(@showprogress [s=>ppfilt[s].(lags,lags').-λ[s] for s in species])
+
+##
+pgram = Dict([s=>periodogram_iso(lansing_points[s], region, 1:20) for s in species])
+
+folder = "fig/fig3/"
+clims = extrema(stack([pp_w[sp] for sp in species]))
+for sp in species
+    figure_p = Figure(resolution = (100,100), fontsize=8, figure_padding=1)
+    ax_p = Axis(figure_p[1,1], limits=(0,1,0,1))
+    hidedecorations!(ax_p)
+    viz!(ax_p, lansing_points[sp], pointsize = 2, color = :black)
+    save(folder * "pts_$sp.pdf", figure_p)
+
+    figure_p = Figure(resolution = (100,60), fontsize=8, figure_padding=1)
+    ax_p = Axis(figure_p[1,1], xticks = [0], yticks=[0], xgridcolor=:black, ygridcolor=:black)
+    ylims!(ax_p, (-1,15))
+    xlims!(ax_p, (-1,21))
+    band!(ax_p, [0,6], fill(-1,2), fill(15,2), color = (:lightgrey,0.5))
+    vlines!(ax_p, [0,6], color = :black, linestyle=:dash, linewidth=1/2)
+    lines!(ax_p, pgram[sp].t, pgram[sp].sdf./pgram[sp].lambda.-1, color = :black)
+    hidedecorations!(ax_p, grid=false)
+    hidespines!(ax_p)
+    save(folder * "sdf_$sp.pdf", figure_p)
+
+    figure_p = Figure(resolution = (100,100), fontsize=8, figure_padding=1)
+    ax_p = Axis(figure_p[1,1])
+    hidedecorations!(ax_p)
+    cf=contourf!(ax_p, lags, lags, pp_w[sp], colorrange = clims, colormap=:tempo)
+    save(folder * "filt_$sp.pdf", figure_p)
+end
+
+@info "finished plotting."

Fig4_impulseresponses.jl

@@ -0,0 +1,47 @@
+using CairoMakie, PointProcessFilters, Meshes, ProgressMeter
+import CairoMakie: Point2f, Circle, Rect
+import PointProcessFilters: CenteredBox, CenteredBall
+import Meshes: Point, Box
+
+## set parameters
+r1 = 5
+r2 = 10
+rsmall = (r2 - r1) / 2
+defaults = (color = (:gray), strokewidth = 0)
+
+## make figure
+figure = Figure(resolution = 300 .* (1, 2 / 4), figure_padding = 1)
+ax = [Axis(figure[1, i], aspect = 1, limits = (r2 + 1) .* (-1, 1, -1, 1)) for i in 1:4]
+hidedecorations!.(ax)
+hidespines!.(ax)
+poly!(ax[1], Circle(Point2f(0, 0), r2); defaults...)
+poly!(ax[1], Circle(Point2f(0, 0), r1); color = :white)
+
+poly!(ax[2], Circle(Point2f(5, r1 + rsmall), rsmall); defaults...)
+poly!(ax[2], Circle(Point2f(-5, -r1 - rsmall), rsmall); defaults...)
+
+poly!(ax[3], Rect(-r2, -r2, 2r2, 2r2); defaults...)
+poly!(ax[3], Rect(-r1, -r1, 2r1, 2r1); color = :white)
+
+poly!(ax[4], Rect(5 - rsmall, r1, 2rsmall, 2rsmall); defaults...)
+poly!(ax[4], Rect(-5 + rsmall, -r1, -2rsmall, -2rsmall); defaults...)
+
+regions = (CenteredBall{2}(r2) - CenteredBall{2}(r1),
+           Ball(Point(5, r1 + rsmall), rsmall),
+           CenteredBox(Point(r2, r2)) - CenteredBox(Point(r1, r1)),
+           Box(Point(5 - rsmall, r1), Point(5 + rsmall, r1 + 2rsmall)))
+h = impulse_response.(RegionPassFilter.(regions))
+x = range(-0.5, 0.5, 100)
+h_evaluated = @showprogress "evaluating ir... " [Makie.pseudolog10.(h[i].(x, x'))
+                                                 for i in 1:4]
+
+ax = [Axis(f[2, i], aspect = 1) for i in 1:4]
+hidedecorations!.(ax)
+hidespines!.(ax)
+for i in 1:4
+    heatmap!(ax[i], x, x, h_evaluated[i], colormap = :tempo,
+             colorrange = extrema(stack(h_evaluated)), interpolate = true)
+end
+
+save("fig/fig4/transfer_functions.pdf", figure)
+figure

LICENSE

@@ -0,0 +1,21 @@
+MIT License
+
+Copyright (c) 2023 Jake P. Grainger
+
+Permission is hereby granted, free of charge, to any person obtaining a copy
+of this software and associated documentation files (the "Software"), to deal
+in the Software without restriction, including without limitation the rights
+to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+copies of the Software, and to permit persons to whom the Software is
+furnished to do so, subject to the following conditions:
+
+The above copyright notice and this permission notice shall be included in all
+copies or substantial portions of the Software.
+
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
+SOFTWARE.