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

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

残されたもう一つのﾊﾞｸﾞがわかりました。
今後の事を考え、ｵﾘｼﾞﾅﾙとは異なる処理をさせていた部分の
渡すべき引数を間違えていました･･･。
"符号なしInt を Long型で代用" とか "右ﾋﾞｯﾄｼﾌﾄ関数" とか
他人のせいにばかりしていましたが、自身の間違えでした。ｽｲﾏｾﾝ。

空間座標(-10,-10,-10) - (10,10,10)
分割ﾚﾍﾞﾙ　2
の時の
座標値(1,1,1) と座標値(9,9,9)
の、点のﾓｰﾄﾝ順序は
Lv2時の "56" と "63"
で、ﾎﾞﾘｭｰﾑとしての線形8分木ｲﾝﾃﾞｯｸｽNoは
"8" で正しいです。　やっと出来ました。
と思っていたのですが、別の座標値では正しくない結果が･･･。

空間座標と分割ﾚﾍﾞﾙは同じままで
座標値(9,9,9) と座標値(9.1,9,9)
を試すと

ﾎﾞﾘｭｰﾑとしての線形8分木ｲﾝﾃﾞｯｸｽNoが "8" と返ってきますが
泥臭い対応表から正しくは "72" が欲しいんです。
実はこの場合は、二つの点が同一の最小ﾚﾍﾞﾙ空間に収まった状態
なんです。

無い知識をフル回転してｵﾘｼﾞﾅﾙｺｰﾄﾞを見ても、同一の最小ﾚﾍﾞﾙ
空間に収まった状態に対しての処理が、書かれていません。
恐らく、全ての要素が最小ﾚﾍﾞﾙ空間を跨いでいる状態を前提として
書かれている感じがします。

その為、それに対応したｺｰﾄﾞとしてこのようにしてみました。

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

Private Const CLINER8TREEMANAGER_MAXLEVEL = 7   '有効空間分割最大ﾚﾍﾞﾙ
Private m_iPow&()
Private m_dwCellNum&                            '空間の数
Private m_CellAry&()                            '線形空間配列 ppCellAryの代わり

Sub CATMain()
   Dim Min: Min = Array(-10#, -10#, -10#)   '最小位置
   Dim Max: Max = Array(10#, 10#, 10#)      '最大位置
   Dim lv&: lv = 2                          '分割ﾚﾍﾞﾙ
   Dim pnt: pnt = Array(9, 9, 9)            '調べたい座標値
   Dim pnt2: pnt2 = Array(9.1, 9, 9)        '調べたい座標値
   
   Dim w: w = ArySub(Max, Min)
   Dim unit: unit = AryDiv(w, CDbl(sl(1, lv)))
   Dim morton&: morton = GetPointElem(pnt, Min, unit)
   Dim morton2&: morton2 = GetPointElem(pnt2, Min, unit)
   
   If Not Init(lv) Then Exit Sub
   Dim spaNo: spaNo = GetMortonNumber(pnt2, pnt, Min, unit, lv)
   Dim msg$
   msg = "最小位置 : " + Join(Min, ",") + vbNewLine + _
         "最大位置 : " + Join(Max, ",") + vbNewLine + _
         "分割ﾚﾍﾞﾙ : " + CStr(lv) + vbNewLine + _
         "の8分木空間で" + vbNewLine + _
         "座標値1  : " + Join(pnt, ",") + vbNewLine + _
         "のﾓｰﾄﾝ順序は [ " + CStr(morton) + " ] " + vbNewLine + _
         "座標値2  : " + Join(pnt2, ",") + vbNewLine + _
         "のﾓｰﾄﾝ順序は [ " + CStr(morton2) + " ]" + vbNewLine + _
         "線形8分木ｲﾝﾃﾞｯｸｽNo: " + CStr(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

'小数点以下切り捨て
'http://www.openreference.org/articles/view/27
Private Function RoundDown(ByVal n#) As Double
    If n > 0 Then
        RoundDown = IIf(n - Int(n) = 0, n, Int(n))
    Else
        RoundDown = IIf(n - Fix(n) = 0, n, Fix(n))
    End If
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
    
    '各レベルでの空間数を算出
    ReDim m_iPow(CLINER8TREEMANAGER_MAXLEVEL + 1)
    m_iPow(0) = 1
    Dim i&
    For i = 1 To UBound(m_iPow)
        m_iPow(i) = m_iPow(i - 1) * 8
    Next
    
    'Levelレベル（0基点）の配列作成
    m_dwCellNum = (m_iPow(Level + 1) - 1) / 3
    ReDim m_CellAry(m_dwCellNum)
    For i = 0 To m_dwCellNum
        m_CellAry(i) = -1
    Next
    Init = True
End Function

'座標から線形8分木ｲﾝﾃﾞｯｸｽNoを取得
''' @param :p1-array(Double)-座標値1
'''        :p2-array(Double)-座標値2
'''        :Min-array(Double)-空間最小座標値
'''        :unit-array(Double)-空間単位ｻｲｽﾞ
'''        :Level-array(Double)-分割ﾚﾍﾞﾙ
''' @return:Long
Private Function GetMortonNumber(ByVal p1, ByVal p2, ByVal Min, ByVal unit, ByVal Level&) As Long
    '最小ﾚﾍﾞﾙにおける空間番号を算出
    Dim LT&: LT = GetPointElem(p1, Min, unit)
    Dim RB&: RB = GetPointElem(p2, Min, unit)
    
    '空間番号を排他的論理和して
    '最上位区切りから所属レベルを算出
    Dim Def&: Def = RB Xor LT
    Dim HiLevel&
    If Def = 0 Then
        HiLevel = 0 '最小ﾚﾍﾞﾙで同一空間番号の場合
    Else
        HiLevel = 1
        Dim i&, Check&
        For i = 0 To Level
            Check = (sr(Def, (i * 3))) And &H7
            If Not Check = 0 Then HiLevel = i + 1
        Next
    End If
    
    'ﾎﾞﾘｭｰﾑの空間番号とﾚﾍﾞﾙから線形8分木ｲﾝﾃﾞｯｸｽNoを算出
    Dim SpaceNum&: SpaceNum = sr(RB, (HiLevel * 3))
    Dim AddNum&: AddNum = (m_iPow(Level - HiLevel) - 1) / 7
    SpaceNum = SpaceNum + AddNum
    
    '不正ﾃﾞｰﾀﾁｪｯｸ
    If SpaceNum > m_dwCellNum Then
        GetMortonNumber = &HFFFFFFFF
    End If
    GetMortonNumber = SpaceNum
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

'座標→線形8分木要素番号変換関数
''' @param :p-array(Double)-座標値
'''        :Min-array(Double)-空間最小座標値
'''        :unit-array(Double)-空間単位ｻｲｽﾞ
''' @return:Long
Private Function GetPointElem(ByVal p As Variant, _
                              ByVal Min As Variant, _
                              ByVal unit As Variant) As Long
    GetPointElem = Get3DMortonNumber( _
            RoundDown((p(0) - Min(0)) / unit(0)), _
            RoundDown((p(1) - Min(1)) / unit(1)), _
            RoundDown((p(2) - Min(2)) / unit(2)))
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

' 右シフト(算術(>>)ではなく論理(>>>)シフトに相当)
Private Function sr(ByVal x&, ByVal n&) As Long
    If n = 0 Then
        sr = x
    Else
        Dim y: y = x And &H7FFFFFFF
        Dim z
        If n = 32 - 1 Then
            z = 0
        Else
            z = y \ CLng(2 ^ n) 'ひょっとしたらCLng要らないかも
        End If
        If y <> x Then z = z Or CLng(2 ^ (32 - n - 1))
        sr = z
    End If
End Function

2点が "同一の最小ﾚﾍﾞﾙ空間に収まった状態" は、こちらのXorした際の "Def" が
"0" の時です。

    ･･･
Private Function GetMortonNumber(ByVal p1, ByVal p2, ByVal Min, ByVal unit, ByVal Level&) As Long
    '最小ﾚﾍﾞﾙにおける空間番号を算出
    Dim LT&: LT = GetPointElem(p1, Min, unit)
    Dim RB&: RB = GetPointElem(p2, Min, unit)
    
    '空間番号を排他的論理和して
    '最上位区切りから所属レベルを算出
    Dim Def&: Def = RB Xor LT
    ･･･

(排他的論理和　って日本語考えた人、もっと短くして欲しかったです)

こちらを実行すると

線形8分木ｲﾝﾃﾞｯｸｽNoは "72" が得られました。　他にも幾つかの座標値で
試したところ正しい値が返って来ました。
やっと入り口にたどり着きました。まだまだこれから･･･。
(このﾍﾟｰｽだと、年内に実装するのは無理っぽいです)