Fibonacci Sampling Method (#1092)

* ✨ Implementing Fibonacci sampling method for disk, ball and sphere.

* 🧪 Adding tests for Fibonacci sampling.

* 🚚 Refactoring Fibonnaci Sampling. Added box method, and capacity to alter the ratio number. Also, sampling now preserves the lentype and numtype for the underlying geometry.

* 🐛 Fixed ustrip warning.

* 🧪 Updated tests for fibonacci sampling.

* 📝 Docs for fibonacci sampling.

* 📝 Moving Fibonacci Sampling to ContinuousSampling

* 🐛 Correcting allocations and fibonacci constructor.

* 🚚 Small formatting correction.

* 🚚 Refactor using the geom parametrization

* 🐛 removing comments and adding distortion for sphere.

* 🚚 Adding Ball{🌐} and Sphere{𝔼{3}}

* 🧪 Added tests for ArgumentError cases.

* 🐛 Fixed formatting

* Update src/sampling/fibonacci.jl

* Refactoring

* Adding license.

* Update src/sampling/fibonacci.jl

* 🧪 changing test for disk fibonacci sampling.

* Refactor tests

* Fix T in test

---------

Co-authored-by: Júlio Hoffimann <[email protected]>

docs/src/algorithms/sampling.md

@@ -113,4 +113,19 @@ sampler = MinDistanceSampling(3.0)
 points  = sample(grid, sampler) |> collect
 
 viz(points)
-```
+```
+
+### FibonacciSampling
+```@docs
+FibonacciSampling
+```
+
+```@example sampling
+sphere = Sphere((0.,0.,0.), 1.)
+
+# sample points using the Fibonacci lattice method
+sampler  = FibonacciSampling(100)
+points  = sample(sphere, sampler) |> collect
+
+viz(points)
+```

src/Meshes.jl

@@ -485,6 +485,7 @@ export
   RegularSampling,
   HomogeneousSampling,
   MinDistanceSampling,
+  FibonacciSampling,
   sampleinds,
   sample,
 

src/sampling.jl

@@ -65,6 +65,7 @@ sample(rng::AbstractRNG, g::Geometry, method::ContinuousSamplingMethod) = sample
 include("sampling/regular.jl")
 include("sampling/homogeneous.jl")
 include("sampling/mindistance.jl")
+include("sampling/fibonacci.jl")
 
 # ----------
 # UTILITIES

src/sampling/fibonacci.jl

@@ -0,0 +1,50 @@
+# ------------------------------------------------------------------
+# Licensed under the MIT License. See LICENSE in the project root.
+# ------------------------------------------------------------------
+
+"""
+	FibonacciSampling(n, ϕ = (1 + √5)/2)
+
+Generate `n` Fibonacci points with parameter `ϕ`.
+
+The golden ratio is used as the default value of `ϕ`,
+but other irrational numbers can be used.
+
+See <https://observablehq.com/@meetamit/fibonacci-lattices>
+and <https://www.johndcook.com/blog/2023/08/12/fibonacci-lattice>.
+"""
+struct FibonacciSampling{T<:Real} <: ContinuousSamplingMethod
+  n::Int
+  ϕ::T
+
+  function FibonacciSampling(n::Int, ϕ::T) where {T<:Real}
+    if n ≤ 0
+      throw(ArgumentError("Size must be positive"))
+    end
+    new{T}(n, ϕ)
+  end
+end
+
+FibonacciSampling(n::Int) = FibonacciSampling(n, (1 + √5) / 2)
+
+function sample(geom::Geometry, method::FibonacciSampling)
+  if paramdim(geom) != 2
+    throw(ArgumentError("Fibonacci sampling only defined for 2D geometries"))
+  end
+
+  fib = _fibmap(geom)
+
+  function point(i)
+    u = mod(i / method.ϕ, 1)
+    v = i / (method.n - 1)
+    geom(fib(u, v)...)
+  end
+
+  (point(i) for i in 0:(method.n - 1))
+end
+
+_fibmap(g) = (u, v) -> (u, v)
+_fibmap(d::Disk) = (u, v) -> (√u, v)
+_fibmap(b::Ball{𝔼{2}}) = (u, v) -> (√u, v)
+_fibmap(b::Ball{🌐}) = (u, v) -> (√u, v)
+_fibmap(s::Sphere{𝔼{3}}) = (u, v) -> (acos(1 - 2v) / π, u)

test/sampling.jl

@@ -387,6 +387,51 @@ end
   @test length(ps) > 0
 end
 
+@testitem "FibonacciSampling" setup = [Setup] begin
+  @test_throws ArgumentError sample(Box(cart(0, 0), cart(1, 1)), FibonacciSampling(-1))
+  @test_throws ArgumentError sample(Box(Point(0, 0, 0), Point(1, 1, 1)), FibonacciSampling(100))
+
+  box = Box(cart(1, 1), cart(4, 2))
+  ps = sample(box, FibonacciSampling(100)) |> collect
+  @test first(ps) isa Point
+  @test first(ps) ≈ cart(1, 1)
+  @test all(∈(box), ps)
+
+  box = Box(cart(0, 0), cart(1, 1))
+  ps = sample(box, FibonacciSampling(100, π)) |> collect
+  @test first(ps) isa Point
+  @test all(∈(box), ps)
+  @test ps[2] ≈ cart(mod(1 / π, 1), 1 / 99)
+
+  tbox = Box(cart(0, 0), cart(1, 1))
+  af = Affine(T[1 1; 0 1], T[2, 0])
+  tbox = af(tbox)
+  ps = sample(tbox, FibonacciSampling(100)) |> collect
+  @test first(ps) isa Point
+  @test first(ps) ≈ af(cart(0, 0))
+  @test all(∈(tbox), ps)
+
+  disk = Disk(Plane(cart(3, 0, 0), Vec(1, 0, 0)), T(2))
+  ps = sample(disk, FibonacciSampling(100)) |> collect
+  @test first(ps) isa Point
+  @test first(ps) ≈ centroid(disk)
+  @test all(p -> coords(p).x ≈ 3u"m", ps)
+  @test all(p -> -2u"m" < coords(p).y || coords(p).y < 2u"m" || isapprox(coords(p).y, 2u"m"; atol=1e-5u"m"), ps)
+  @test all(p -> -2u"m" < coords(p).z || coords(p).z < 2u"m" || isapprox(coords(p).z, 2u"m"; atol=1e-5u"m"), ps)
+
+  sphere = Sphere(cart(1, 1, 1), T(2))
+  ps = sample(sphere, FibonacciSampling(100)) |> collect
+  @test first(ps) isa Point
+  @test first(ps) ≈ cart(1, 1, 3)
+  @test all(∈(sphere), ps)
+
+  ball = Ball(cart(2, 1), T(0.1))
+  ps = sample(ball, FibonacciSampling(100)) |> collect
+  @test first(ps) isa Point
+  @test first(ps) ≈ centroid(ball)
+  @test all(∈(ball), ps)
+end
+
 @testitem "RNGs" setup = [Setup] begin
   dom = cartgrid(100, 100)
   for method in [UniformSampling(100), WeightedSampling(100), BallSampling(T(10))]