Changeset 26709

Timestamp:

12/06/21 06:56:35 (3 years ago)

Author:

Mathieu Morlighem

Message:

CHG: added ContourToNodes

File:

• issm/trunk-jpl/src/jl/md/exp.jl (modified) (1 diff)

issm/trunk-jpl/src/jl/md/exp.jl

@@ -136 +136 @@
         return contours
 end# }}}
+#ContourToNodes{{{
+"""
+        ContourToNodes - Flag points that are in contour
+
+   More doc to come later....
+
+        Usage:
+                exp=expread(filename)
+
+# Examples:
+```julia-repl
+julia> exp=expread('domainoutline.exp')
+```
+
+# Arguments:
+- filename: the ARGUS file to read
+"""
+function ContourToNodes(x::Vector{Float64},y::Vector{Float64},filename::String,edgevalue::Float64)
+
+        #Read input file
+        contours = expread(filename)
+
+        #Initialize output
+        nbpts = length(x)
+        flags = zeros(Bool,nbpts)
+
+        #Loop over contours
+        for c in 1:length(contours)
+
+                #Get current contours
+                contour = contours[c]
+                xp      = contour.x
+                yp      = contour.y
+
+                #Check that we are within box
+                xmin = minimum(xp); xmax = maximum(xp)
+                ymin = minimum(yp); ymax = maximum(yp)
+
+                #Loop over all points provided
+                for ii in 1:nbpts
+
+                        #If this node is already within one of the contours, do not change it
+                        if(flags[ii]) continue end
+
+                        #Are we within bounds?
+                        if(x[ii]<xmin || x[ii]>xmax || y[ii]<ymin || y[ii]>ymax) continue end
+
+                        #we are potentially inside... perform pnpoly test
+                        flags[ii] = pnpoly(xp, yp, x[ii], y[ii], edgevalue)
+                end
+        end
+        
+        return flags
+end# }}}
+
+function pnpoly(xp::Vector{Float64},yp::Vector{Float64},x::Float64,y::Float64,edgevalue::Float64) #{{{
+
+        npol = length(xp)
+
+        #Do we need to test for colinearity?
+        if(edgevalue!=2)
+                i = 1
+                j = npol
+                while(i<=npol)
+
+                        n1 = (yp[i]-yp[j])^2 + (xp[i]-xp[j])^2
+                        n2 = (y-yp[j])^2 + (x-xp[j])^2
+
+                        normp=sqrt(n1*n2)
+                        scalar=(yp[i]-yp[j])*(y-yp[j])+(xp[i]-xp[j])*(x-xp[j])
+
+                        if (scalar == normp)
+                                if (n2<=n1)
+                                        return edgevalue
+                                end
+                        end
+
+                        j =  i
+                        i += 1
+                end
+        end
+
+        #second test : point is neither on a vertex, nor on a side, where is it ?
+        i = 1
+        j = npol
+        c = false
+        while(i<=npol)
+                if (((yp[i]<=y && y<yp[j]) || (yp[j]<=y && y<yp[i])) &&
+            (x < (xp[j] - xp[i]) * (y - yp[i]) / (yp[j] - yp[i]) + xp[i]))
+         c = !c
+                end
+
+                j =  i
+                i += 1
+        end
+        return c
+end# }}}