umegusa's blog

備忘録

3D convex hull

3次元凸包アルゴリズムを組んでみました。
こちらを参考に実装してみました。
一応メモ

ちなみにRubyなのでかなり遅いです。
ぶっちゃけ計算向きではないです。

入力は3次元の座標。
出力は凸包を形成する3点になります。

convexhull.rb

require './point3D'
require './face'

## incremental convexhull

# 3D convex hull
# param: points the Point3D coord
# return: a series of indices to create the faces of the hull

def convexhull(points)
	if points.length < 3 then return nil end
	
	face = Face.new(points)
	face.vertices(0,1,2)
	
	# 最初に4面体を作り、計算をする。
	# 頂点が4つ未満の場合、nullを返す
	v = centroid(points, 3, face)
	
	if face.isVisible(v) then face.flip end
	
	face0 = Face.new(points)
	face0.vertices(3, face.i0, face.i1)
	if face0.isVisible(v) then face0.flip end
	
	face1 = Face.new(points)
	face1.vertices(3, face.i1, face.i2)
	if face1.isVisible(v) then face1.flip end

	face2 = Face.new(points)
	face2.vertices(3, face.i2, face.i0)
	if face2.isVisible(v) then face2.flip end
	
	# 4面体を形成する
	valid_faces = [face,face0,face1,face2]
	visible_faces = []
	tmp_faces = []
	result = Hash::new
	
	# すべての頂点に処理を掛ける。最初の4つの点は省く
	(4..points.length-1).each{ |i|
		v = points[i]
		
		visible_faces = []
		
		valid_faces.each{ |face|
			if face.isVisible(v) then
				visible_faces << face
			end
		}
		
		# 4面体の頂点の中にあるため、次の点で比較
		if visible_faces.length == 0 then next end
		
		# 頂点が凸砲の外側にある場合、リストから表示されているすべての面を削除
		visible_faces.each{ |face|
			delete_id = valid_faces.index(face)
			valid_faces.delete_at(delete_id)
		}
		
		# 片面だけの場合、3面で構成しても大丈夫なので3面で構成
		if visible_faces.length == 1 then
			face = visible_faces[0]
			input_face = Face.new(points)
			input_face.vertices(i, face.i0, face.i1)
			valid_faces << input_face
			
			input_face = Face.new(points)
			input_face.vertices(i, face.i1, face.i2)
			valid_faces << input_face
			
			input_face = Face.new(points)
			input_face.vertices(i, face.i2, face.i0)
			valid_faces << input_face
			next
		end
		
		# すべての頂点が既存のものより外側にあった場合、すべての面から新しく作る
		tmp_faces = []
		count = 0
		
		# 新しいvisible_facesを作成する
		visible_faces.each{ |face|
			input_face = Face.new(points)
			
			input_face.vertices(i, face.i0, face.i1)
			tmp_faces << input_face
			
			input_face = Face.new(points)
			input_face.vertices(i, face.i1, face.i2)
			tmp_faces << input_face
			
			input_face = Face.new(points)
			input_face.vertices(i, face.i2, face.i0)
			tmp_faces << input_face
			count = count + 1
		}
		
		
		tmp_faces.each{ |face|
			flag = true
			
			# 面の前に点がないかを検索する
			tmp_faces.each{ |other|
				if flag then
					if face != other then
						if face.isVisible(other.centroid) then
							face = nil
							break
							flag = false
						end
					end
				end
			}
			
			# 面の前に点が無い場合格納する
			if !face.nil? then
				valid_faces << face
			end
		}
		
	}
	
	h_count = 0
	valid_faces.each{ |face|
		input_array = [face.i0, face.i1, face.i2]
		result[h_count] = input_array
		h_count = h_count + 1
	}
	
	# 3D covex hullの結果を返す
	return result

end

def centroid(points, index, face)
	p = points[index]
	p0 = points[face.i0]
	p1 = points[face.i1]
	p2 = points[face.i2]
	return Point3D.new( (p.x + p0.x + p1.x + p2.x) / 4, (p.y + p0.y + p1.y + p2.y) / 4, (p.z + p0.z + p1.z + p2.z) / 4)
end


point = []

open("./input/test.txt",'r'){|strs|
  while line = strs.gets
  	str = line.split(" ")
  	coord = Point3D.new(str[0].to_f,str[1].to_f, str[2].to_f)
  	point << coord
  end
}

puts point.size

hull = convexhull(point)
points = []
(0..10).each{|i|
	points << point[i]
}


p hull

point3D.rb

class Point3D
  def initialize(x,y,z)
    @x = x
    @y = y
    @z = z
  end

  # return Point3D
  def distance(p)
    dx = @x - p.x
    dy = @y - p.y
    dz = @z - p.z
    return Math::sqrt(dx ** 2 + dy ** 2 + dz ** 2)
  end
  
  # 外積
  # return Point3D
  def cross(p)
  	a = @y * p.z - @z * p.y
  	b = @z * p.x - @x * p.z
  	c = @x * p.y - @y * p.x
  	return Point3D.new(a,b,c)
  end
  
  # 内積
  # a*b = axbx + ayby + azbz
  # return Double
  def dot(p)
	return @x * p.x + @y * p.y + @z * p.z
  end

  # ベクトルの大きさ
  # return Double
  def vectorLength
    return Math::sqrt(@x ** 2 + @y ** 2 + @z ** 2)
  end

  # 正規化 : 単位ベクトルに正規化する
  # return Point3D
  def norm
  	point = Point3D.new(@x,@y,@z)
  	len = 1.0 / point.vectorLength
  	@x = @x / len
  	@y = @y / len
  	@z = @z / len
  end
  
  def subtract(p)
    a = @x - p.x
    b = @y - p.y
    c = @z - p.z
    return Point3D.new(a,b,c)
  end
  
  # 右手法か判断
  # 正規化した後、計算する
  # z軸を中心にする。
  # z軸が正の時→x,yが正ならば1
  # z軸が負の時→x,yが負ならば1  を返す
  def systemComparison
    r_value = 0
#    if @x * @y == @z or 
#       @y * @z == @x or 
#       @z * @x == @y then r_value = 1 end
	if @z > 0 and @x > 0 and @y > 0 then return 0 end
	if @z < 0 and @x < 0 and @y < 0 then return 0 end
       
#    puts "#{@x} #{@y} #{@z} "
 
    return r_value
  end
  
  attr_reader :x, :y, :z
end


face.rb

require './point3d'

class Face
	def initialize(point)
		@points = point
	end
	
	def vertices(i0, i1, i2)
		@i0 = i0
		@i1 = i1
		@i2 = i2
		
		computePanel
	end
	
	def computePanel
		p1 = @points[@i0]
		p2 = @points[@i1]
		p3 = @points[@i2]
		
		@a = p1.y * (p2.z - p3.z) + p2.y * (p3.z - p1.z) + p3.y * (p1.z - p2.z)
		@b = p1.z * (p2.x - p3.x) + p2.z * (p3.x - p1.x) + p3.z * (p1.x - p2.x)
		@c = p1.x * (p2.y - p3.y) + p2.x * (p3.y - p1.y) + p3.x * (p1.y - p2.y)
		@d = -(p1.x * (p2.y * p3.z - p3.y * p2.z) + p2.x * (p3.y * p1.z - p1.y * p3.z) + p3.x * (p1.y * p2.z - p2.y * p1.z))
		
	end
	
	def isVisible(point)
		return (@a * point.x + @b * point.y + @c * point.z + @d) > 0
	end
	
	def outputtest(point)
		puts "a::#{@a} b::#{@b} c::#{@c} d::#{@d}"
		puts "x::#{point.x} y::#{point.y} z::#{point.z}"
		return @a * point.x + @b * point.y + @c * point.z + @d
	end
	
	def centroid
		p0 = @points[@i0]
		p1 = @points[@i1]
		p2 = @points[@i2]
		return Point3D.new( (p0.x + p1.x + p2.x)/3.0, (p0.y + p1.y + p2.y)/3.0, (p0.z + p1.z + p2.z)/3.0)
	end
	
	def flip
		t = @i0
		@i0 = @i1
		@i1 = t
		computePanel
	end
	attr_reader :points, :i0, :i1, :i2
end