C#ATIA

↑タイトル詐欺 主にFusion360API 偶にCATIA V5 VBA(絶賛ネタ切れ中)

単独な3D曲線の取得する2

こちらの続きです。
単独な3D曲線の取得する1 - C#ATIA


前回は総当りの非常に効率の悪い物でしたが、線形8分木による
スクリーニング処理を利用し単独な3D曲線を検索するように修正しました。

まずはCATMainを持つ標準モジュール "FindSingle3DCurve.bas" です。

'vba FindSingle3DCurve Ver0.0.2  using-'KCL0.09' 'FindSingle3DCurve_Octree0.0.1'
'単独な3D曲線の取得

Option Explicit

'*** 設定 ***
Private Const Tolerance = 0.001                         '端点一致トレランス
Private Const OctLv = 2&                                '8分木分割レベル(最大7)
Private Const InSpaceMaxCount = 200&                    '同一空間内最大数(目安)
Private Const SingleCurveFooter = " - SingleCurve!!"    '単独曲線名に追記する文字
'************

Sub CATMain()
    'ドキュメントのチェック
    If Not CanExecute(Array("PartDocument", "ProductDocument")) Then Exit Sub
    Dim Msg$
    
    '形状セットの選択
    Msg = "3D曲線の入っている形状セットを指定してください : ESCキー 終了"
    Dim Hb As HybridBody
    Set Hb = KCL.SelectItem(Msg, "HybridBody")
    If KCL.IsNothing(Hb) Then Exit Sub
    
    '線の取得
    Dim Lines As Collection: Set Lines = GetLines(Hb)
    If KCL.IsNothing(Lines) Then
        Msg = "指定した' " + Hb.Name + " 'には、3D曲線がありませんでした!"
        MsgBox Msg: Exit Sub
    End If
    
    '警告
    Msg = CStr(Lines.Count) + "本あります。実行しますか?"
    
    If MsgBox(Msg, vbYesNo + vbInformation) = vbNo Then Exit Sub
    Debug.Print "** start ** : " + CStr(Lines.Count) + "本" 'debug
    Debug.Print "OctLv : " + CStr(OctLv) + "  InSpaceMaxCount : " + CStr(InSpaceMaxCount) 'debug
    KCL.SW_Start 'debug
    
    '端点座標取得
    Dim EndPnts As Collection: Set EndPnts = GetEndPntCoordinates(Lines)
    Debug.Print "EndPnts : " + CStr(KCL.SW_GetTime) 'debug
    
    '組み合わせリスト取得
    Dim OctList As Collection
    If FindSingle3DCurve_Octree.Init(OctLv, Tolerance, InSpaceMaxCount) Then
        Set OctList = FindSingle3DCurve_Octree.GetLinerOctreeList(EndPnts)
    End If
    
    If KCL.IsNothing(OctList) Then
        Msg = "エラーです。マクロ内の設定値を修正して下さい(特に 'OctLv' )"
        MsgBox Msg: Exit Sub
    End If
    Debug.Print "OctList : " + CStr(KCL.SW_GetTime) 'debug
    
    '単独曲線インデックス取得
    Dim SingleList As Collection: Set SingleList = GetSingleIdx(EndPnts, OctList)
    Debug.Print "SingleList : " + CStr(SingleList.Count) + " : " + CStr(KCL.SW_GetTime)
    If SingleList.Count < 1 Then
        Msg = "指定した' " + Hb.Name + " 'には、単独な3D曲線が見つかりませんでした!" + vbNewLine + _
              "(トレランス " & Tolerance & " mm)"
        MsgBox Msg: Exit Sub
    End If
    
    '確認
    Msg = SingleList.Count & "本の単独な3D曲線が見つかりました!!" + vbNewLine + _
          "名前の最後に ' " + SingleCurveFooter + " ' を追記しますか?"
    If MsgBox(Msg, vbYesNo + vbInformation) = vbNo Then GoTo MacroEnd
    
    'リネーム
    Call AppendName(Lines, SingleList)
    
MacroEnd:
    'Call Dump(ByVal SingleList, ByVal Lines, ByVal EndPnts) 'debug
    MsgBox "終了!"
End Sub

'単独曲線のリネーム
Private Sub AppendName(ByVal Lines As Collection, ByVal IdxList As Collection)
    Dim i
    For Each i In IdxList
        With Lines.Item(i)
            .Name = .Name & SingleCurveFooter
        End With
    Next
End Sub

'初期化済み配列生成 - オブジェクト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 GetSingleIdx(ByVal EndPnts As Collection, ByVal OctList As Collection) As Collection
    Dim i&, j&, i_crv_idx, j_crv_idx
    Dim CombAry As Variant: CombAry = InitRangeAry(EndPnts.Count * 0.5, False)
    Dim Space As Collection
    For Each Space In OctList
        For i = 1 To Space.Count
            i_crv_idx = EndPnts.Item(Space.Item(i))(3)
            For j = i + 1 To Space.Count
                j_crv_idx = EndPnts.Item(Space.Item(j))(3)
                If i_crv_idx = j_crv_idx Then GoTo Continue
                If CombAry(i_crv_idx) And CombAry(j_crv_idx) Then GoTo Continue
                If IsPosEqual(EndPnts.Item(Space.Item(i)), EndPnts.Item(Space.Item(j))) Then
                    CombAry(i_crv_idx) = True
                    CombAry(j_crv_idx) = True
                End If
Continue:
            Next
        Next
    Next
    
    Dim SingleIdxList As Collection: Set SingleIdxList = New Collection
    For i = 1 To UBound(CombAry)
        If CombAry(i) = False Then SingleIdxList.add i
    Next
    Set GetSingleIdx = SingleIdxList
End Function

'座標値一致
Private Function IsPosEqual(ByVal p1 As Variant, ByVal p2 As Variant) As Boolean
    IsPosEqual = False
    Dim i&
    For i = 0 To 2
        If Abs(p2(i) - p1(i)) > Tolerance Then Exit Function
    Next
    If Sqr(Abs((p2(0) - p1(0)) * (p2(0) - p1(0)) + _
               (p2(1) - p1(1)) * (p2(1) - p1(1)) + _
               (p2(2) - p1(2)) * (p2(2) - p1(2)))) < Tolerance Then
        IsPosEqual = True
    End If
End Function

'端点取得
Private Function GetEndPntCoordinates(ByRef Lines As Collection) As Collection
    Dim Doc As PartDocument: Set Doc = KCL.GetParent_Of_T(Lines.Item(1), "PartDocument")
    Dim SPAWB As Workbench: Set SPAWB = Doc.GetWorkbench("SPAWorkbench")
    Dim EndPnt As Collection: Set EndPnt = New Collection
    Dim Pos(8) As Variant
    Dim i&
    For i = 1 To Lines.Count
        Call SPAWB.GetMeasurable(Lines.Item(i)).GetPointsOnCurve(Pos)
        EndPnt.add KCL.JoinAry(KCL.GetRangeAry(Pos, 0, 2), Array(i))
        EndPnt.add KCL.JoinAry(KCL.GetRangeAry(Pos, 6, 8), Array(i))
    Next
    Set GetEndPntCoordinates = EndPnt
End Function

'線取得
Private Function GetLines(ByVal Hb As HybridBody) As Collection
    Dim Hss As HybridShapes: Set Hss = Hb.HybridShapes
    If Hss.Count < 1 Then Exit Function
    
    Dim Pt As Part: Set Pt = KCL.GetParent_Of_T(Hb, "PartDocument").Part
    Dim Fact As HybridShapeFactory: Set Fact = Pt.HybridShapeFactory
    Dim Lines As Collection: Set Lines = New Collection
    Dim Hs As HybridShape
    For Each Hs In Hss
        Select Case Fact.GetGeometricalFeatureType(Hs)
            Case 2, 3, 4
                Lines.add Hs
        End Select
    Next
    Set GetLines = Lines
End Function

'debug用
Private Sub Dump(ByVal SingleList, ByVal Lines, ByVal EndPnts)
    Dim ss
    For Each ss In SingleList
        Debug.Print ss & ":" & Lines.Item(EndPnts.Item(ss)(3)).Name
    Next
End Sub

'ver0.0.1 - 総当りのダメなやつ
'ver0.0.2 - 線形8分木を実装

設定を行う為の定数部分の説明です。
・"Tolerance"
 端点一致を判断するトレランスです。一般的に考え0.001以下は無意味な気がしています。
 又、データの状態次第ですが、空間分割する際にも影響がある為あまり
 大きな値(100.0とか)では、処理時間が長くなる可能性があります。

・"OctLv"
 空間分割する際の分割レベルです。テストを行った感じでは大きくした方が処理時間が短く
 なる傾向にあります。 最大7です。

・"InSpaceMaxCount"
 空間分割した際、同一空間内の要素数の最大数ですが、あくまで目安です。
 こちらは値は、大小どちらにした方が効率良くなるのか正直把握しておりません。
 
・"SingleCurveFooter"
 前回同様、Hitした単独曲線名に追記する文字です。


続いて、組み合わせ効率を改善するための標準モジュール "FindSingle3DCurve_Octree.bas" です。
あまりに複雑な上、他にも流用できそうな気がしたので、単独なモジュールとしました。

'vba FindSingle3DCurve_Octree Ver0.0.1  using-'KCL0.08'
'モートン順序を利用した8分木空間分割ライブラリ
'FindSingle3DCurveマクロ専用です

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&()                          'レベル毎の空間数

'*** 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)
    '座標値郡のIdxList作成
    Dim PntIdxList As Collection: Set PntIdxList = InitRangeList(Pnts.Count)
    
    '空間内最大数以下の場合、そのまま返す
    Dim DecidedList As Collection: Set DecidedList = New Collection
    If Pnts.Count < m_MaxCount Then
        Call DecidedList.add(PntIdxList)
        GoTo FuncEnd
    End If
    
    '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個しかないものは検証しない
                        'Call DecidedList.add(TempSpace)
                    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
    Set PntIdxList = Nothing
    Set DecidedList = Nothing
    Set CheckList = Nothing
    Set ReCheckList = Nothing
    Set TempList = Nothing
    ReDim SpeceEnum(0)
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関数に、分割レベル,・トレランス,・空間内最大数 を渡して初期化し
GetLinerOctreeList関数に座標値の入ったコレクションを渡すことで、
組み合わせを行うための、座標値のインデックスが空間毎に入ったコレクションを返します。
引数-Collection(Array(x,y,z-Double)) の形で配列は3次元以上でもOKです。
戻り値-Collection(Collection(Long))です。返すのはあくまでも列挙する為の
インデックスで、座標値ではありません。

後は、標準モジュール KCL が必要になります。
f:id:kandennti:20161228121941p:plain
こんな感じで3つが同一プロジェクトに入っている必要があります。


これを、前回(Ver0.0.1)の物と比較しながら実際に処理した結果です。

Ver0.0.1
** start ** : 4100本
EndPnts : 1.73
SingleList : 1 : 2771.944

Ver0.0.2
** start ** : 4100本
OctLv : 2  InSpaceMaxCount : 200
EndPnts : 1.771
OctList : 3.607
SingleList : 1 : 9.34

Ver0.0.2
** start ** : 4100本
OctLv : 7  InSpaceMaxCount : 200
EndPnts : 2.214
OctList : 3.33
SingleList : 1 : 4.128

SingleList行の最後の値が、トータルの処理時間です。
4099本→4100本に1本増やしたのは、使用していたサンプルデータに単独な
3D曲線が無かった為です。
f:id:kandennti:20161228122006p:plain

4100本を処理した際、Ver0.0.1では45分以上の処理時間が必要だったのですが
Ver0.0.2では、分割レベル2で9秒 分割レベル7では4秒程まで短縮しました。
線形8分木の効果が絶大です。 これであれば十分耐えられる処理時間かと
思います。
但し、"FindSingle3DCurve.bas" 側ではもう少し効率化出来そうなアイデアが
幾つかあるので、今後の課題とします。


正直に書くと、線形8分木に取り組んだのは初めてでしたが、線形4分木に
ついては過去にDrawの重複線削除マクロで取り組んだことがあります。
ですが、線形4分木の衝突リストを作成する為の空間の巡回方法がどうしても
理解できず、情けない話断念した経緯があり、今回も実装までたどり着けるか
どうか当初は解りませんでした。
(このブログを始めた理由の一つは、重複線削除マクロをC#化し線形4分木を
実装することでした。 …全くやっていませんが)

線形8分木の実装は年内では間に合わないと思っていたのですが、こちらの
思い付きの効果が大きく
端点一致を探る、組み合わせテスト10 - C#ATIA
ギリギリですが、年内には出来上がったつもりです。(バグがあるかも…)


今年はこれでおしまいです。昨年も書きましたが、こんなにブログを続けると
思っていませんでしたよ。