Add isinside functionality (#96)

Author: dkarrasch

.github/workflows/ci.yml

@@ -30,7 +30,7 @@ jobs:
       - uses: julia-actions/julia-buildpkg@v1
       - uses: julia-actions/julia-runtest@v1
       - uses: julia-actions/julia-processcoverage@v1
-      - uses: codecov/codecov-action@v2
+      - uses: codecov/codecov-action@v3
         if: ${{ matrix.version == '1' && matrix.os == 'ubuntu-latest' }}
         with:
           file: lcov.info

Project.toml

@@ -1,6 +1,6 @@
 name = "CoherentStructures"
 uuid = "0c1513b4-3a13-56f1-9cd2-8312eaec88c4"
-version = "0.4.10"
+version = "0.4.11"
 
 [deps]
 ArnoldiMethod = "ec485272-7323-5ecc-a04f-4719b315124d"

docs/Manifest.toml

@@ -147,8 +147,8 @@ compute_closed_orbits
 First of all, [`Singularity`](@ref), [`EllipticBarrier`](@ref) and [`EllipticVortex`](@ref)
 objects can be passed as initial conditions to the [`flow`](@ref) function, which returns
 a (time-resolved) vector of corresponding objects. Moreover, we have advection with
-adaptive point insertion (known in oceanography as "Dritschel advection") and an function
-to compute the area of [`EllipticBarrier`](@ref) and [`EllipticVortex`](@ref) objects.
+adaptive point insertion (known in oceanography as "Dritschel advection") and some more
+convenience functions.
 
 ```@docs
 FlowGrowParams
@@ -157,4 +157,5 @@ area
 centroid
 clockwise
 Base.extrema(::EllipticBarrier)
+isinside
 ```

docs/src/elliptic.md

@@ -1755,3 +1755,39 @@ function clockwise(barrier::EllipticBarrier, velo::ODEFunction{false}, t0; p=not
 end
 clockwise(vortex::EllipticVortex, velo, t0; p=nothing) =
     clockwise(last(vortex.barriers), velo, t0; p=p)
+
+"""
+    isinside(p, polygon)
+
+Determine whether the point `p` (or [`Singularity`](@ref)) lies inside `polygon`, which can
+be of type `Vector{SVector{2}}`, `EllipticBarrier` or `EllipticVortex`. In the latter case,
+the outermost (i.e., the last `EllipticBarrier` in the `barriers` field) is used.
+"""
+isinside(p::SVector{2}, poly::Vector{<:SVector{2}}) = !iszero(winding_number(p, poly))
+isinside(p::SVector{2}, barrier::EllipticBarrier) = isinside(p, barrier.curve)
+isinside(p::SVector{2}, vortex::EllipticVortex) = isinside(p, last(vortex.barriers))
+isinside(s::Singularity, curve) = isinside(s.coords, curve)
+
+function quadrant(p::SVector{2})
+    if p[1] > 0
+        p[2] >= 0 && return 1
+        return 4
+    else
+        p[2] >= 0 && return 2
+        return 3
+    end
+end
+function winding_number(q::SVector{2}, poly::Vector{<:SVector{2}})
+    qprev = qfirst = quadrant(poly[1] - q)
+    qnext = quadrant(poly[2] - q)
+    Δq = qnext - qprev
+    @inbounds for i in 3:length(poly)
+        qnext = quadrant(poly[i] - q)
+        qdiff = qnext - qprev
+        Δq += qdiff == 3 ? -1 : qdiff == -3 ? 1 : qdiff
+        qprev = qnext
+    end
+    qdiff = qfirst - qnext
+    Δq += qdiff == 3 ? -1 : qdiff == -3 ? 1 : qdiff
+    return Δq ÷ 4
+end

src/ellipticLCS.jl

@@ -31,6 +31,7 @@ export
 	area,
 	centroid,
 	clockwise,
+	isinside,
 
 	#diffusion_operators.jl
 	gaussian,

src/exports.jl

@@ -173,16 +173,19 @@ end
     @test area(square) ≈ 4 rtol=5eps()
     @test centroid(square) ≈ ones(2) atol=5eps()
     barrier = CS.EllipticBarrier(circle, SVector{2}(0.0, 0.0), 1.0, true)
-    vortices = CS.EllipticVortex(SVector{2}(0.0, 0.0), [barrier])
+    vortex = CS.EllipticVortex(SVector{2}(0.0, 0.0), [barrier])
     LL, UR = extrema(barrier)
-    @test area(vortices) == area(barrier) == area(circle)
-    @test centroid(vortices) == centroid(barrier) == centroid(circle)
-    @test extrema(vortices) == extrema(barrier)
+    @test area(vortex) == area(barrier) == area(circle)
+    @test centroid(vortex) == centroid(barrier) == centroid(circle)
+    @test extrema(vortex) == extrema(barrier)
     @test LL ≈ -ones(2) rtol=1e-5
     @test UR ≈ ones(2) rtol=1e-6
     @test !clockwise(barrier, ODEFunction((u, p, t) -> SVector{2}(-u[2], u[1])), 0)
     @test clockwise(barrier, ODEFunction((u, p, t) -> SVector{2}(-u[2], u[1])), 0) ==
-        clockwise(vortices, ODEFunction((u, p, t) -> SVector{2}(-u[2], u[1])), 0)
+        clockwise(vortex, ODEFunction((u, p, t) -> SVector{2}(-u[2], u[1])), 0)
+    @test isinside(SVector(1-eps(),0), vortex)
+    @test isinside(Singularity(vortex.center, 1), vortex)
+    @test !isinside(SVector(1+eps(),0), vortex)
     barrier = CS.EllipticBarrier(square, SVector{2}(1.0, 1.0), 1.0, true)
     LL, UR = extrema(barrier)
     @test LL ≈ zeros(2)