Merge pull request #82 from ArneSpang/main

Added inPolygon with tests

Author: boriskaus

src/utils.jl

@@ -7,6 +7,7 @@ export RotateTranslateScale
 export LithostaticPressure!
 export FlattenCrossSection
 export AddField, RemoveField
+export inPolyPoint, inPolyPointF, inPolygon!
 
 using NearestNeighbors
 
@@ -1640,4 +1641,106 @@ function LithostaticPressure!(Plithos::Array{T,N}, Density::Array{T,N}, dz::Numb
     Plithos[:] = reverse!(cumsum(reverse!(Plithos),dims=N))
 
     return nothing
+end
+
+"""
+    inPolygon!(INSIDE::Matrix, PolyX::Vector, PolyY::Vector, X::Matrix, Y::Matrix; fast=false)
+
+Checks if points given by matrices `X` and `Y` are in or on (both cases return true) a polygon given by `PolyX` and `PolyY`. Boolean `fast` will trigger faster version that may miss points that are exactly on the edge of the polygon. Speedup is a factor of 3.
+
+"""
+function inPolygon!(INSIDE::Matrix{Bool}, PolyX::Vector{T}, PolyY::Vector{T}, X::Matrix{T}, Y::Matrix{T}; fast=false) where T <: Real
+    iSteps = collect(eachindex(PolyX))
+    jSteps = [length(PolyX); collect(1:length(PolyX)-1)]
+
+    if fast
+        for j = 1 : size(X, 2)
+            for i = 1 : size(X, 1)
+                INSIDE[i,j] = inPolyPointF(PolyX, PolyY, X[i,j], Y[i,j], iSteps, jSteps)
+            end
+        end
+    else
+        for j = 1 : size(X, 2)
+            for i = 1 : size(X, 1)
+                INSIDE[i,j] = (inPolyPoint(PolyX, PolyY, X[i,j], Y[i,j], iSteps, jSteps) || inPolyPoint(PolyY, PolyX, Y[i,j], X[i,j], iSteps, jSteps))
+            end
+        end
+    end
+end
+
+"""
+    inPolygon!(inside::Vector, PolyX::Vector, PolyY::Vector, x::Vector, y::Vector; fast=false)
+
+Same as above but `inside`, `X` and `Y` and are vectors.
+
+"""
+function inPolygon!(inside::Vector{Bool}, PolyX::Vector{T}, PolyY::Vector{T}, x::Vector{T}, y::Vector{T}; fast=false) where T <: Real
+    iSteps = collect(eachindex(PolyX))
+    jSteps = [length(PolyX); collect(1:length(PolyX)-1)]
+
+    if fast
+        for i = eachindex(x)
+            inside[i] = inPolyPointF(PolyX, PolyY, x[i], y[i], iSteps, jSteps)
+        end
+    else
+        for i = eachindex(x)
+            inside[i] = (inPolyPoint(PolyX, PolyY, x[i], y[i], iSteps, jSteps) || inPolyPoint(PolyY, PolyX, y[i], x[i], iSteps, jSteps))
+        end
+    end
+end
+
+"""
+    inPolyPoint(PolyX::Vector, PolyY::Vector, x::Number, y::Number, iSteps::Vector, jSteps::)
+
+Checks if a point given by x and y is in or on (both cases return true) a polygon given by PolyX and PolyY, iSteps and jSteps provide the connectivity between the polygon edges. This function should be used through inPolygon!().
+
+"""
+function inPolyPoint(PolyX::Vector{T}, PolyY::Vector{T}, x::T, y::T, iSteps::Vector{Int64}, jSteps::Vector{Int64}) where T <: Real
+    inside1, inside2, inside3, inside4 = false, false, false, false
+    for (i,j) in zip(iSteps, jSteps)
+        xi = PolyX[i]
+        yi = PolyY[i]
+        xj = PolyX[j]
+        yj = PolyY[j]
+
+        con1 = ((yi > y) != (yj > y))
+        con2 = ((yi >= y) != (yj >= y))
+        if con1 && (x > (xj - xi) * (y - yi) / (yj - yi + eps()) + xi)
+            inside1 = !inside1
+        end
+
+        if con1 && (x >= (xj - xi) * (y - yi) / (yj - yi + eps()) + xi)
+            inside2 = !inside2
+        end
+
+        if con2 && (x > (xj - xi) * (y - yi) / (yj - yi + eps()) + xi)
+            inside3 = !inside3
+        end
+
+        if con2 && (x >= (xj - xi) * (y - yi) / (yj - yi + eps()) + xi)
+            inside4 = !inside4
+        end
+    end
+    return ((inside1 || inside2) || (inside3 || inside4))
+end
+
+"""
+    inPolyPointF(PolyX::Vector, PolyY::Vector, x::Number, y::Number, iSteps::Vector, jSteps::)
+
+Faster version of inPolyPoint() but will miss some points that are on the edge of the polygon.
+
+"""
+function inPolyPointF(PolyX::Vector{T}, PolyY::Vector{T}, x::T, y::T, iSteps::Vector{Int64}, jSteps::Vector{Int64}) where T <: Real
+    inside = false
+    for (i,j) in zip(iSteps, jSteps)
+        xi = PolyX[i]
+        yi = PolyY[i]
+        xj = PolyX[j]
+        yj = PolyY[j]
+
+        if ((yi > y) != (yj > y)) && (x > (xj - xi) * (y - yi) / (yj - yi + eps()) + xi)
+            inside = !inside
+        end
+    end
+    return inside
 end

test/test_utils.jl

@@ -239,3 +239,15 @@ cross_tmp = CrossSection(Data_EQ,Lon_level=16.4,section_width=10km)
 cross_tmp = CrossSection(Data_EQ,Start=(15.0,35.0),End=(17.0,37.0),section_width=10km)
 @test cross_tmp.fields.lon_proj[20] ==15.314329874961091 
 @test cross_tmp.fields.lat_proj[20] == 35.323420618580585
+
+# test inPolygon
+PolyX = [-2.,-1,0,1,2,1,3,3,8,3,3,1,2,1,0,-1,-2,-1,-3,-3,-8,-3,-3,-1,-2]
+PolyY = [3.,3,8.01,3,3,1,2,1,0,-1,-2,-1,-3,-3,-8,-3,-3,-1,-2,-1,0,1,2,1,3]
+xvec  = collect(-9:0.5:9); yvec = collect(-9:0.5:9); zvec = collect(1.:1.);
+X,Y,Z = meshgrid(xvec, yvec, zvec)
+X, Y  = X[:,:,1], Y[:,:,1]
+yN    = zeros(Bool, size(X))
+inPolygon!(yN, PolyX, PolyY, X, Y, fast=true)
+@test sum(yN) == 194
+inPolygon!(yN, PolyX, PolyY, X, Y)
+@test sum(yN) == 217