another convex hull / GrahamScan function added

src/hulls/graham.jl

@@ -62,3 +62,84 @@ function hull(pset::PointSet{2,T}, ::GrahamScan) where {T}
 
   PolyArea(c)
 end
+
+
+"""
+  hull_simple(p, h, s, v)
+
+calculate the convex hull for the matrix of points.  This is actually an implementation of the pseudocode from  [https://en.wikipedia.org/wiki/Graham_scan]  (https://en.wikipedia.org/wiki/Graham_scan).
+
+This is a fast and convinient method when point coordinates are in a matrix, and you have to create 10000000 hulls in a second.  
+
+parameters:   
+ * p is the matrix [y x w]  of the points; where 'y' is the point dimension and w is the number of points.  
+ * h is the preallocated vector of the indices of the points in p [w]
+ * s is the preallocated vector of the indices of the points in p, which makes a hull   [w]
+ * v preallocated vector of Folat64 [w]
+
+return the number of the points in a hull, and indices of the hull points in 's'
+"""
+function hull_simple(p, h, s, v)
+	j3, n = size(p)
+	p0 = 1
+
+	# detect p0
+	@inbounds for i=2:n
+		if p[2, i] < p[2, p0]   # minimum y
+			p0 = i
+		elseif p[2, i] == p[2, p0]
+			if p[1, i] < p[1, p0]
+				p0 = i
+			end
+		end
+	end
+	#@info "p0 = $p0"
+
+	@inbounds for i = 1:n
+		if i == p0
+			v[i] = 1.0
+			continue
+		end
+		v[i] = (((p[1, i] - p[1, p0])) / sqrt((p[1, i] - p[1, p0])^2 + (p[2, i] - p[2, p0])^2)) # this is the cosinus of the angle
+	end
+	#@info "not sorted v = $v"
+	sortperm!(h, v, rev = true)			
+	#@info "sorted v = $v"
+	u = 0
+	@inbounds for i = 1:n
+		# v = (p2[1] - p1[1])*(p3[2] - p1[2]) - (p2[2] - p1[2])*(p3[1] - p1[1])
+		while u > 1 &&  (p[1, s[u]] - p[1, s[u-1]])*(p[2, h[i]] - p[2, s[u-1]]) - (p[2, s[u]] - p[2, s[u-1]])*(p[1, h[i]] - p[1, s[u-1]]) <= 0.0
+			u -= 1
+		end
+		u += 1
+		s[u] = h[i]
+	end
+	return u
+end
+
+
+"""
+	zeroInside(p, s)
+
+check if zero point [0.0, 0.0] is inside the convex hull (hull created by hull_simple() function). 
+
+ * p is the matrix [y x w]  of the points; where 'y' is the point dimension and w is the number of points.  
+ * s is the vector of indices in p, of the points that we use for convex hull
+ * u number of the indices in 's'
+
+ return true if zero point is inside the hull.
+"""
+function zeroInside(p, s, u)
+	for i = 2:u
+		if (p[1, s[i]] - p[1, s[i-1]]) * (-p[2, s[i-1]]) - (p[2, s[i]] - p[2, s[i-1]])*(-p[1, s[i-1]]) < 0.0
+			return false
+		end
+	end
+	# the last segment:
+	if (p[1, s[1]] - p[1, s[u]])*(-p[2, s[u]]) - (p[2, s[1]] - p[2, s[u]])*(-p[1, s[u]]) < 0.0
+		return false
+	end
+	return true
+end
+
+export hull_simple, zeroInside

test/hulls.jl

@@ -35,4 +35,15 @@
     verts = vertices(boundary(chul))
     @test verts == P2[(0.5,-1),(1,0),(1,1),(0,1),(0,0)]
   end
+  @testset "GrahamScan_simple" begin
+    p = randn(2, 25)
+    #p .+= [2.0, 0.5]
+    _, n = size(p)
+    h = zeros(Int32, n)
+    v = zeros(Float64, n)
+    s = similar(h)
+    u = hull_simple(p, h, s, v)
+    #s = s[1:u]
+    @test zeroInside(p, s, u)
+  end
 end