samples_2.4 => samples

share/grafx2/scripts/samples/picture/thomson/lib/convex_hull.lua

@@ -0,0 +1,219 @@
+-- convxhull.lua : support for computing the convex
+-- hull of a set of points.
+--
+-- inspired from: https://gist.github.com/anonymous/5184ba0bcab21d3dd19781efd3aae543
+--
+-- Version: 02-jan-2017
+--
+-- Copyright 2016-2017 by Samuel Devulder
+--
+-- This program is free software; you can redistribute
+-- it and/or modify it under the terms of the GNU
+-- General Public License as published by the Free
+-- Software Foundation; version 2 of the License.
+-- See <http://www.gnu.org/licenses/>
+
+if not ConvexHull then
+
+local function sub(u,v)
+	return {u[1]-v[1],u[2]-v[2],u[3]-v[3]}
+end
+
+local function mul(k,u)
+	return {k*u[1],k*u[2],k*u[3]}
+end
+
+local function cross(u,v)
+	return {u[2]*v[3] - u[3]*v[2],
+	        u[3]*v[1] - u[1]*v[3],
+			u[1]*v[2] - u[2]*v[1]}
+end
+
+local function dot(u,v)
+	return u[1]*v[1] + u[2]*v[2] + u[3]*v[3]
+end
+
+local function unit(u)
+	local d=dot(u,u)
+	return d==0 and u or mul(1/d^.5, u)
+end
+
+ConvexHull = {}
+
+function ConvexHull:new(coordFct)
+	local o = {
+		points={},
+		coord=coordFct
+	}
+	setmetatable(o, self)
+	self.__index = self
+	return o
+end
+
+function ConvexHull.coord(elt)
+	return {elt[1],elt[2],elt[3]}
+end
+
+function ConvexHull:vect(a,b)
+	return sub(self.coord(b),self.coord(a))
+end
+
+function ConvexHull:normal(face)
+	local u=self:vect(face[1],face[2])
+	local v=self:vect(face[1],face[3])
+	return cross(u,v)
+end
+
+function ConvexHull:printPoint(p)
+	return '('..table.concat(self.coord(p),',')..')'
+end
+
+function ConvexHull:printFace(F)
+	return '['..self:printPoint(F[1])..' '..
+	            self:printPoint(F[2])..' '..
+				self:printPoint(F[3])..']'
+end
+
+function ConvexHull:seen(face,p)
+	local N=self:normal(face)
+	local P=self:vect(face[1],p)
+	return dot(N,P)>=0
+end
+
+function ConvexHull:bdry(faces)
+	local code={n=0}
+	function code.encode(pt,...)
+		if pt then
+			local k = code[pt]
+			if not k then
+				k = code.n+1
+				code[k]  = pt
+				code[pt] = k
+				code.n   = k
+			end
+			local rest = code.encode(...)
+			return rest and (k..','..rest) or ""..k
+		end
+	end
+	function code.decode(str)
+		local i = str:find(',')
+		if i then
+			local k = str:sub(1,i-1)
+			return code[tonumber(k)],code.decode(str:sub(i+1))
+		else
+			return code[tonumber(str)]
+		end
+	end
+	local set = {}
+	local  function add(...)
+		set[code.encode(...)] = true
+	end
+	local function rem(...)
+		set[code.encode(...)] = nil
+	end
+	local function keys()
+		local r = {}
+		for k in pairs(set) do
+			r[{code.decode(k)}] = true
+		end
+		return r
+	end
+	for F in pairs(faces) do
+		add(F[1],F[2])
+		add(F[2],F[3])
+		add(F[3],F[1])
+	end
+	for F in pairs(faces) do
+		rem(F[1],F[3])
+		rem(F[3],F[2])
+		rem(F[2],F[1])
+	end
+	return keys()
+end
+
+function ConvexHull:addPoint(p)
+	-- first 3 points
+	if self.points then
+		if p==self.points[1] or p==self.points[2] then return end
+		table.insert(self.points,p)
+
+		if #self.points==3 then
+			self.hull={
+				{self.points[1],self.points[2],self.points[3]},
+				{self.points[1],self.points[3],self.points[2]}
+			}
+			self.points=nil
+		end
+	else
+		local seenF,n = {},0
+		for _,F in ipairs(self.hull) do
+			if F[1]==p or F[2]==p or F[3]==p then return end
+			if self:seen(F,p) then seenF[F]=true;n=n+1 end
+		end
+
+		if n==#self.hull then
+			-- if can see all faces, unsee ones looking "down"
+			local N
+			for F in pairs(seenF) do N=self:normal(F); break; end
+			for F in pairs(seenF) do
+				if dot(self:normal(F),N)<=0 then
+					seenF[F] = nil
+					n=n-1
+				end
+			end
+		end
+
+		-- remove (old) seen faces
+		local z=#self.hull
+		for i=#self.hull,1,-1 do
+			if seenF[self.hull[i]] then
+				table.remove(self.hull,i)
+			end
+		end
+
+		-- insert new boundaries with seen faces
+		for E in pairs(self:bdry(seenF)) do
+			table.insert(self.hull,{E[1],E[2],p})
+		end
+	end
+	return self
+end
+
+function ConvexHull:verticesSet()
+	local v = {}
+	if self.hull then
+		for _,F in ipairs(self.hull) do
+			v[F[1]] = true
+			v[F[2]] = true
+			v[F[3]] = true
+		end
+	end
+	return v
+end
+
+function ConvexHull:verticesSize()
+	local n = 0
+	for _ in pairs(self:verticesSet()) do n=n+1 end
+	return n
+end
+
+function ConvexHull:distToFace(F,pt)
+	local N=unit(self:normal(F))
+	local P=self:vect(F[1],pt)
+	return dot(N,P)
+end
+
+function ConvexHull:distToHull(pt)
+	local d
+	for _,F in ipairs(self.hull) do
+		local t = self:distToFace(F,pt)
+		d = d==nil and t or
+			(0<=t and t<d or
+			 0>=t and t>d) and t or
+			d
+		if d==0 then break end
+	end
+	return d
+end
+
+end -- ConvexHull