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

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

昨日の悩みが全て解消するような、良い方法を思い付きました。
時間が無い為、先にｺｰﾄﾞです。

'vba
'ﾓｰﾄﾝ順序を利用した8分木空間分割の為のｸｽﾞｺｰﾄﾞ3
Option Explicit

Private Const CLINER8TREEMANAGER_MAXLEVEL = 7   '有効空間分割最大ﾚﾍﾞﾙ

Sub CATMain()
   Dim min: min = Array(-40#, -40#, -40#)   '最小位置
   Dim Max: Max = Array(40#, 40#, 40#)      '最大位置
   Dim lv&: lv = 2                          '分割ﾚﾍﾞﾙ
   Dim Pnt: Pnt = Array(0, 0, 0)            '調べたい座標値
   Dim Tolerance#: Tolerance = 0.001        '端点一致ﾄﾚﾗﾝｽ
   
   Dim w: w = ArySub(Max, min)
   Dim unit: unit = AryDiv(w, CDbl(sl(1, lv)))
   
   If Not Init(lv) Then Exit Sub
   Dim SpaNo As Variant: SpaNo = GetMortonNum(Pnt, min, unit, Tolerance)
   Dim msg$
   msg = "最小位置 : " + Join(min, ",") + vbNewLine + _
         "最大位置 : " + Join(Max, ",") + vbNewLine + _
         "分割ﾚﾍﾞﾙ : " + CStr(lv) + vbNewLine + _
         "の8分木空間で" + vbNewLine + _
         "座標値  : " + Join(Pnt, ",") + vbNewLine + _
         "線形8分木ｲﾝﾃﾞｯｸｽNo: " + Join(SpaNo, " , ")
   MsgBox msg
End Sub

'*** 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 AryDiv(ByVal a, ByVal b#) As Variant
    AryDiv = Array(a(0) / b, a(1) / b, a(2) / b)
End Function


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

'線形8分木配列を構築
Private Function Init(ByVal Level&) As Boolean
    Init = False
    
    '設定最高レベル以上の空間は作れない
    If Level >= CLINER8TREEMANAGER_MAXLEVEL Then Exit Function
    
    Init = True
End Function

'ビット分割関数
Private Function BitSeparateFor3D(ByVal n As Byte) 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分木モートン順序算出関数
Private Function Get3DMortonNumber(ByVal x As Byte, ByVal y As Byte, ByVal z As Byte) As Long
   Get3DMortonNumber = BitSeparateFor3D(x) Or _
                       sl(BitSeparateFor3D(y), 1) Or _
                       sl(BitSeparateFor3D(z), 2)
End Function

'座標→最小ﾚﾍﾞﾙ空間番号郡取得
''' @param :p-array(Double)-座標値
'''        :Min-array(Double)-空間最小座標値
'''        :unit-array(Double)-空間単位ｻｲｽﾞ
'''        :tol-Double-ﾄﾚﾗﾝｽ
''' @return:array(Long)
Private Function GetMortonNum(ByVal p As Variant, _
                              ByVal min As Variant, _
                              ByVal unit As Variant, _
                              ByVal tol As Double) As Variant
    '最小位置から各軸の比率算出
    Dim ratio#(2), inte&(2), dec#(2) '比率,整数,小数
    Dim tolratio#(2) '各軸のﾄﾚﾗﾝｽの比率
    Dim i&
    For i = 0 To 2
        ratio(i) = (p(i) - min(i)) / unit(i)
        inte(i) = Fix(ratio(i))
        dec(i) = ratio(i) - inte(i)
        tolratio(i) = tol / unit(i)
    Next
    
    'ﾄﾚﾗﾝｽを考慮
    Dim axis(2) As Variant
    Dim spapos As Collection
    For i = 0 To 2
        Set spapos = New Collection
        spapos.add inte(i)
        If dec(i) <= tolratio(i) Then spapos.add inte(i) - 1
        If 1# - dec(i) <= tolratio(i) Then spapos.add inte(i) + 1
        Set axis(i) = spapos
    Next
    
    '空間番号算出
    Dim x, y, z
    Dim SpaNo(): ReDim SpaNo(axis(0).Count * axis(1).Count * axis(2).Count)
    Dim j&, k&, 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)
    Next: Next: Next
    ReDim Preserve SpaNo(cnt)
    GetMortonNum = SpaNo
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

細かなﾁｪｯｸ、ｾﾝｽの無い変数・関数名は見逃してください。

今までは点の座標値とﾄﾚﾗﾝｽを考慮し、ﾎﾞﾘｭｰﾑとして線形8分木のｲﾝﾃﾞｯｸｽを
求めていましたが、ﾄﾚﾗﾝｽを考慮して点として最小ﾚﾍﾞﾙのﾓｰﾄﾝ順序の空間番号
を求めるようにしました。(わかり難い表現です･･･)

上記のｺｰﾄﾞを実行した結果はこちら

空間の重心位置(要はど真ん中)の場合、以前ではﾙｰﾄ空間に属する事に
なっていました。　この場合、他に生成された全ての空間と端点一致の
処理が必要でした。　それに対し今回は、最低ﾚﾍﾞﾙの空間の8ヶ所のみの
ﾁｪｯｸで済みます。　この方法であれば
・空間ﾎﾟｲﾝﾀ配列(ｵﾘｼﾞﾅﾙのppCellAry)が、最低ﾚﾍﾞﾙ分のみ
・衝突ﾘｽﾄの巡回が、再帰・非再帰とかではなく単純なﾙｰﾌﾟ1回のみ
・同一空間での要素数が多い際、再度空間分割することが容易
　（極端に離れた要素がある場合への対応が楽）
他にも、ｵﾘｼﾞﾅﾙｺｰﾄﾞではﾁｪｯｸ漏れを起こす可能性がある場合に気が付き
それにも対応したつもりです。

ｵﾘｼﾞﾅﾙｺｰﾄﾞをVBA化する事ばかり考えていたので、計算式の意味を
把握していませんでした。　VBAの商が切り下げじゃないと騒いだ時に
計算式の意味を良く考えれば良かったです。