C#ATIA

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

流用して、4分木る1

こちらで遅かった為の改善策を模索中です。
2D曲線の折れ線化を利用し、重複線の選択1 - C#ATIA

曲線の折れ線化については、後回しにして4分木に挑戦です。
8分木の計算式変更して、次元を一つ落とせば良いのだろうと
思っているのですが。

'vba test_Quadtree Ver0.0.1  using-'KCL0.09'
'モートン順序を利用した4分木空間分割テスト用クズコード

Option Explicit

Private Const MAXLEVEL = 2                      '有効空間分割最大レベル

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 Pnts As Collection: Set Pnts = New Collection
    With Pnts
        .Add Array(-40#, -40#)
        .Add Array(40#, 40#)
        .Add Array(0#, 0#)
        .Add Array(25#, -25#)
        .Add Array(1#, 1#)
    End With
    
    '初期設定
    Dim Tol#: Tol = 0.001            '端点一致トレランス
    Dim MaxCount&: MaxCount = 2      '空間内最大数
    
    If Not SetStart(Tol, MaxCount) Then
        MsgBox "設定値が不正です"
        Exit Sub
    End If
    
    '座標値郡のIdxList作成
    Dim PntIdxList As Collection: Set PntIdxList = InitRangeList(Pnts.Count)
    Dim Id: Id = SetSpaceInfo(Pnts, PntIdxList)
    
    '確認
    Call DumpQuadIdx(Pnts)
End Sub

'確認用
Private Sub DumpQuadIdx(ByVal PosList)
    Dim Pos
    Debug.Print " **** "
    For Each Pos In PosList
        Debug.Print "Pos:" & Pos(0) & "," & Pos(1), _
                    "MotonNo:" & GetPointElem(Pos)
    Next
End Sub

'線形4分木準備
''' @param :Tolerance-Double-一致トレランス
''' @param :MaxCount-long-同一空間内最大数(目安)
''' @return:Boolean
Private Function SetStart(ByVal Tolerance#, ByVal MaxCount&) As Boolean
    SetStart = False
    If Tolerance <= 0 Then Exit Function
    
    m_Tolerance = Tolerance
    m_MaxCount = MaxCount
    m_ToleranceRatio = InitRangeAry(2, 0)
    
    ReDim m_CellCount(MAXLEVEL + 1)
    m_CellCount(0) = 1
    Dim i&
    For i = 1 To UBound(m_CellCount)
        m_CellCount(i) = m_CellCount(i - 1) * 4
    Next
    SetStart = True
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)
    If Not SetLevel(W) Then Exit Function
    m_Unit = AryDiv(W, m_AxisCount)
    
    Dim i&
    For i = 0 To 1
        m_ToleranceRatio(i) = m_ToleranceRatio(i) / m_Unit(i)
    Next
    SetSpaceInfo = True
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 1
            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

'空間サイズとトレランスからレベル算出し設定
Private Function SetLevel(ByVal W) As Boolean
    SetLevel = False
    Dim Min#: Min = 1.79769313486231E+308
    Dim i&
    For i = 0 To 1
         If Min > W(i) Then Min = W(i)
    Next
    Dim TmpLv&: TmpLv = Fix(Log_n((Min / (m_Tolerance * 2 + 0.002)), 2))
    If TmpLv > MAXLEVEL Then
        m_Level = MAXLEVEL
    Else
        m_Level = TmpLv
    End If
    
    If m_Level < 1 Then Exit Function
    m_AxisCount = sl(1, m_Level)
    SetLevel = True
End Function

'ビット分割関数2D
''' @param :n-long
''' @return:long
Private Function BitSeparateFor2D(ByVal n&) As Long
    Dim S&: S = n
    S = (S Or sl(S, 8)) And &HFF00FF
    S = (S Or sl(S, 4)) And &HF0F0F0F
    S = (S Or sl(S, 2)) And &H33333333
    BitSeparateFor2D = (S Or sl(S, 1)) And &H55555555
End Function

'4分木モートン順序算出関数
''' @param :x-long
''' @param :y-long
''' @return:long
Private Function Get2DMortonNumber(ByVal X&, ByVal Y&) As Long
   Get2DMortonNumber = BitSeparateFor2D(X) Or _
                       sl(BitSeparateFor2D(Y), 1)
End Function

'線形4分木インデックス取得関数
''' @param :Pos-array(Double)
''' @return:long
Private Function GetPointElem(ByVal Pos As Variant) As Long
   GetPointElem = Get2DMortonNumber(Fix((Pos(0) - m_MinPos(0)) / m_Unit(0)), _
                                    Fix((Pos(1) - m_MinPos(1)) / m_Unit(1)))
End Function


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

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

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

'nを底とする対数
Private Function Log_n(X, n)
    Log_n = Log(X) / Log(n)
End Function

' 左シフト
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

'初期化済み配列生成 - オブジェクト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

実行してみた感じはこちら。

Pos:-40,-40   MotonNo:0
Pos:40,40     MotonNo:15
Pos:0,0       MotonNo:12
Pos:25,-25    MotonNo:5
Pos:1,1       MotonNo:12

てっきり上手く行くと思ったのに・・・・。

こちらの画像を参考にすると
その8 4分木空間分割を最適化する!
モートン番号は 10,5,6,15,6 の順で出て欲しいのに。
Y方向が逆になってます。 何故? 何処間違えたのだろう?