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