読者です 読者をやめる 読者になる 読者になる

C#ATIA

↑タイトル詐欺 主にCATIA V5 の VBA

端点一致を探る、組み合わせテスト11

こちらの続きです。
端点一致を探る、組み合わせテスト10 - C#ATIA


前回の物を実装する為に書き換えました。

'vba
'モートン順序を利用した8分木空間分割の為のクズコード4
Option Explicit

Private Const CLINER8TREEMANAGER_MAXLEVEL = 7   '有効空間分割最大レベル

Private m_Level&                                '分割レベル
Private m_Tolerance#                            '一致トレランス
Private m_MaxCount&                             '同一空間内最大数(目安)
Private m_AxisCount&                            '空間分割時の各軸の最大数
Private m_MinPos                                '空間最小座標
Private m_Unit                                  '空間単位サイズ
Private m_ToleranceRatio                        '空間単位サイズに対してのトレランス比率
Private m_CellCount&()                          'レベル毎の空間数

Sub CATMain()
    'テスト座標
    Dim Poss As Collection: Set Poss = New Collection
    With Poss
        .add Array(-40#, -40#, -40#)
        .add Array(40#, 40#, 40#)
        .add Array(0, 0, 0)
        .add Array(0.1, 0, 0)
        .add Array(0.1, 0.001, 0)
        .add Array(-39, -39, -39)
        .add Array(-39, -39.1, -39)
        .add Array(-39, -39, -39.1)
    End With
    
    Dim lv&: lv = 2                  '分割レベル
    Dim Tol#: Tol = 0.001            '端点一致トレランス
    Dim MaxCount&: MaxCount = 3      '空間内最大数
    
    If Not Init(lv, Tol, MaxCount) Then Exit Sub
    
    '分割した空間毎に組み合わせる端点のインデックスを取得
    Dim List As Collection: Set List = GetLinerOctreeList(Poss)
    
    Dim msg$
    msg = "分割レベル : " & lv & vbNewLine + _
          "端点一致トレランス : " & Tol & vbNewLine + _
          "空間内最大数 : " & MaxCount & vbNewLine + _
          "空間数 : " & List.Count & vbNewLine + _
          "以下、組み合わせるリストです" & vbNewLine
    Dim msg2$, tmp: msg2 = ""
    For Each tmp In List
        msg2 = msg2 & List2Str(tmp) & vbNewLine
    Next
    MsgBox msg & msg2
End Sub

'テスト用
Private Function List2Str(ByVal List) As String
    Dim tmp, msg$: msg = ""
    For Each tmp In List
        msg = msg & tmp & ","
    Next
    List2Str = Mid(msg, 1, Len(msg) - 1)
End Function


'*** Octree ***
'http://marupeke296.com/COL_3D_No15_Octree.html

'線形8分木準備
''' @param :Level-long-分割レベル
''' @param :Tolerance-Double-一致トレランス
''' @param :MaxCount-long-同一空間内最大数(目安)
''' @return:Boolean
Function Init(ByVal Level&, ByVal Tolerance#, ByVal MaxCount&) As Boolean
    Init = False
    If Level >= CLINER8TREEMANAGER_MAXLEVEL Then Exit Function
    If Tolerance <= 0 Then Exit Function
    
    m_Level = Level
    m_Tolerance = Tolerance
    m_MaxCount = MaxCount
    m_ToleranceRatio = InitRangeAry(2, 0)
    m_AxisCount = sl(1, m_Level)
    
    ReDim m_CellCount(CLINER8TREEMANAGER_MAXLEVEL + 1)
    m_CellCount(0) = 1
    Dim i&
    For i = 1 To UBound(m_CellCount)
        m_CellCount(i) = m_CellCount(i - 1) * 8
    Next
    Init = True
End Function

'線形8分木リスト取得
''' @param :Pnts-Collection(array(Double))-座標値郡
''' @return:Collection(Collection(long))-空間別点Idx郡
Function GetLinerOctreeList(ByVal Pnts As Collection)
    '空間内最大数以下の場合、そのまま返す
    Dim DecidedList As Collection: Set DecidedList = New Collection
    If Pnts.Count < m_MaxCount Then
        Call DecidedList.add(Pnts)
        GoTo FuncEnd
    End If
    
    '座標値郡のIdxList作成
    Dim PntIdxList As Collection: Set PntIdxList = InitRangeList(Pnts.Count)
    
    'IdxListをルート空間として登録
    Dim CheckList As Collection: Set CheckList = New Collection
    Call CheckList.add(PntIdxList)
    
    'レベルの応じたIdx用配列
    Dim SpeceEnum, TempList As Collection, ReCheckList As Collection
    Dim Space, TempSpace, SpaAry, i&, Idx
    Do
        Set ReCheckList = New Collection
        '配置
        For Each Space In CheckList
            SpeceEnum = InitRangeAry(m_CellCount(m_Level), -1)
            Set TempList = New Collection
            
            If Not SetSpaceInfo(Pnts, Space) Then '空間が小さすぎる
                Call DecidedList.add(Space)
                GoTo Continue
            End If
            
            For Each Idx In Space
                SpaAry = GetMortonNum(Pnts(Idx))
                For i = 0 To UBound(SpaAry)
                    If SpeceEnum(SpaAry(i)) < 0 Then
                        Call TempList.add(InitSpace())
                        SpeceEnum(SpaAry(i)) = TempList.Count
                    End If
                    Call TempList.Item(SpeceEnum(SpaAry(i))).add(Idx)
                Next
            Next
            
            '空間毎の数が多いものを再配置
            For Each TempSpace In TempList
                Select Case True
                    Case TempSpace.Count < 2
                        '空間に1個しかないものは検証しない
                    Case TempSpace.Count > m_MaxCount
                        Call ReCheckList.add(TempSpace)
                    Case Else
                        Call DecidedList.add(TempSpace)
                End Select
            Next
Continue:
        Next
        If ReCheckList.Count < 1 Then Exit Do
        Set CheckList = ReCheckList
    Loop
    
FuncEnd:
    Set GetLinerOctreeList = DecidedList
End Function

'空のコレクション
Private Function InitSpace() As Collection
    Set InitSpace = New Collection
End Function

'空間情報設定
''' @param :Pnts-Collection(array(Double))-座標値郡
''' @param :Idxs-Collection(long)-座標値郡Idx
''' @return:Boolean
Private Function SetSpaceInfo(ByVal Pnts, ByVal Idxs As Collection) As Boolean
    SetSpaceInfo = False
    Dim SpSize: SpSize = GetSpaceSize_Idx(Pnts, Idxs)
    m_MinPos = AryAdd(SpSize(0), m_Tolerance * -1)
    Dim w: w = ArySub(AryAdd(SpSize(1), m_Tolerance), m_MinPos)
    m_Unit = AryDiv(w, m_AxisCount)
    
    Dim MinUnit#: MinUnit = m_Tolerance * 3
    Dim i&
    For i = 0 To 2
         If m_Unit(i) < MinUnit Then Exit Function
    Next
    
    For i = 0 To 2
        m_ToleranceRatio(i) = m_ToleranceRatio(i) / m_Unit(i)
    Next
    SetSpaceInfo = True
End Function

'座標→最小レベル空間番号郡取得
''' @param :Pos-array(Double)-座標値
''' @return:array(Long)
Private Function GetMortonNum(ByVal Pos As Variant) As Variant
    Dim Ratio#(2), Inte&(2), Dec#(2)
    Dim i&
    For i = 0 To 2
        Ratio(i) = (Pos(i) - m_MinPos(i)) / m_Unit(i)
        Inte(i) = Fix(Ratio(i))
        Dec(i) = Ratio(i) - Inte(i)
    Next
    
    Dim Axis(2) As Variant
    Dim AxisNums As Collection
    For i = 0 To 2
        Set AxisNums = New Collection
        AxisNums.add Inte(i)
        If Dec(i) <= m_ToleranceRatio(i) And Inte(i) > 0 Then AxisNums.add Inte(i) - 1
        If 1# - Dec(i) <= m_ToleranceRatio(i) And Inte(i) < m_AxisCount - 1 Then AxisNums.add Inte(i) + 1
        Set Axis(i) = AxisNums
    Next
    
    Dim x, y, z
    Dim SpaNo(): ReDim SpaNo(Axis(0).Count * Axis(1).Count * Axis(2).Count)
    Dim cnt&: cnt = -1
    For Each x In Axis(0)
        For Each y In Axis(1)
            For Each z In Axis(2)
                cnt = cnt + 1
                SpaNo(cnt) = Get3DMortonNumber(x, y, z)
                If SpaNo(cnt) < 0 And SpaNo(cnt) >= m_CellCount(m_Level) Then cnt = cnt - 1
    Next: Next: Next
    ReDim Preserve SpaNo(cnt)
    GetMortonNum = SpaNo
End Function

'ビット分割関数
''' @param :n-long
''' @return:long
Private Function BitSeparateFor3D(ByVal n&) As Long
    Dim s As Long: s = n
    s = (s Or sl(s, 8)) And &HF00F
    s = (s Or sl(s, 4)) And &HC30C3
    s = (s Or sl(s, 2)) And &H249249
    BitSeparateFor3D = s
End Function

'8分木モートン順序算出関数
''' @param :x-long
''' @param :y-long
''' @param :z-long
''' @return:long
Private Function Get3DMortonNumber(ByVal x&, ByVal y&, ByVal z&) As Long
   Get3DMortonNumber = BitSeparateFor3D(x) Or _
                       sl(BitSeparateFor3D(y), 1) Or _
                       sl(BitSeparateFor3D(z), 2)
End Function

'座標値郡から空間サイズ取得
''' @param :EndPnts-Collection(array(Double))-座標値郡
''' @return:array(array(Double))-0:最小値 1:最大値
Private Function GetSpaceSize_Idx(ByVal EndPnts As Collection, ByVal IdxList As Collection) As Variant
    Dim Min: Min = InitRangeAry(2, 1.79769313486231E+308)
    Dim Max: Max = InitRangeAry(2, -1.79769313486231E+308)
    
    Dim Idx, i&
    For Each Idx In IdxList
        For i = 0 To 2
            If Min(i) > EndPnts.Item(Idx)(i) Then Min(i) = EndPnts.Item(Idx)(i)
            If Max(i) < EndPnts.Item(Idx)(i) Then Max(i) = EndPnts.Item(Idx)(i)
        Next
    Next
    GetSpaceSize_Idx = Array(Min, Max)
End Function

'*** BitShift ***
'http://www.geocities.co.jp/SiliconValley/4334/unibon/asp/bitshift2.html
' 左シフト
Private Function sl(ByVal x&, ByVal n&) As Long
    If n = 0 Then
        sl = x
    Else
        Dim k: k = CLng(2 ^ (32 - n - 1))
        Dim d: d = x And (k - 1)
        Dim c: c = d * CLng(2 ^ n)
        If x And k Then c = c Or &H80000000
        sl = c
    End If
End Function

'*** VBA不足関数 ***
'配列同士の引き算-細かいチェック無し
Private Function ArySub(ByVal a, ByVal b) As Variant
    ArySub = Array(a(0) - b(0), a(1) - b(1), a(2) - b(2))
End Function

'配列と実数の足し算-細かいチェック無し
Private Function AryAdd(ByVal a, ByVal b#) As Variant
    AryAdd = Array(a(0) + b, a(1) + b, a(2) + b)
End Function

'配列と実数の割り算-細かいチェック無し
Private Function AryDiv(ByVal a, ByVal b#) As Variant
    AryDiv = Array(a(0) / b, a(1) / b, a(2) / b)
End Function

'初期化済み配列生成 - オブジェクトNG
Private Function InitRangeAry(ByVal Count&, ByVal Value As Variant)
    Dim ary() As Variant: ReDim ary(Count)
    Dim i&
        For i = 0 To Count
            ary(i) = Value
        Next
    InitRangeAry = ary
End Function

'初期化済みコレクション生成
Private Function InitRangeList(ByVal Count&) As Collection
    Dim List As Collection: Set List = New Collection
    Dim i&
    For i = 1 To Count
        List.add i
    Next
    Set InitRangeList = List
End Function

Init関数で線形8分木を作成する為の準備を行い、GetLinerOctreeList関数に
座標値のコレクションを渡す事で、空間分割された点群のインデックスを返します。

もう、GetLinerOctreeList関数の酷さ(3重ループ)と言ったら、言葉も出ない程
なのですが、他に思いつきませんでした・・・。

実際に実行した結果はこんな感じです。
f:id:kandennti:20161227175157p:plain
他の方には意味不明かと思いますが、実はオリジナルより効率が悪い部分が有ります。
両方のアルゴリズムを上手く利用できると効率良く処理できるのですが・・・。


テスト段階ですが、こちらのコードに今回の物を実装したのですが

単独な3D曲線の取得する1 - C#ATIA

劇的に速くなりました。 線形8分木の素晴らしさを感じています。
過去に作ったマクロで遅くて話にならないものがあるのですが、
これを利用すれば使い道が出てきそうな気がしてきました。

トレランスを導入した事で余計に混乱したのですが、以下は自分用に覚書です。
・分割した空間内に1つしか要素が無い場合、リストに入れる必要なし
・分割した空間サイズがトレランスに満たない場合、正しく処理されない可能性あり