曲線と戦ってみる1

ﾀｲﾄﾙを付けたもののｺﾞｰﾙを決めていない為、どんな結果になるのか?
何にたどり着くのか？自分でもわかっていません。

かなり以前の話ですが、"Unofficial CATIA User Forum"　で、
y4yamaさんが三次元曲線の円弧近似ﾏｸﾛを公開してくれていました。
手元にｺｰﾄﾞが残っていた為、転記させてもらいます。
(bunkatu と tole の値はｵﾘｼﾞﾅﾙでは無い可能性があります)

'vbs - Type1
Option Explicit
Dim icd(200, 2), pnt_onCrv(200)
Dim icd1(2), icd2(2), icd3(2)
Dim centre(2), ng_arcCnt, tole
Dim part1, Fact1, Hb1

Sub CATMain() 'catVBA
    Dim bunkatu
    bunkatu = 150: tole = 0.02
    ng_arcCnt = 0
    
    Dim Doc, aSel, ret1, crv1
    Set Doc = CATIA.ActiveDocument
    Set aSel = Doc.selection
    Set part1 = Doc.Part
    Set Fact1 = part1.HybridShapeFactory
    Set Hb1 = part1.HybridBodies.add
    Dim InputObjectType(0)
    InputObjectType(0) = "HybridShape"
    aSel.Clear
    ret1 = aSel.SelectElement2(InputObjectType, "sel CRV", True)
    If (ret1 <> "Normal") Then Exit Sub
    Set crv1 = aSel.Item(1).Value
    aSel.Clear

    Dim i As Integer, ratio As Double, j As Integer
    Dim rad1 As Double
    Dim pnt1, origin(2)
    For i = 0 To bunkatu
        ratio = i / bunkatu
        Set pnt1 = Fact1.AddNewPointOnCurveFromPercent(crv1, ratio, False)
        Hb1.AppendHybridShape pnt1
        part1.Update
        Call pnt1.GetCoordinates(origin)
        For j = 0 To 2: icd(i, j) = origin(j): Next
        'Fact1.DeleteObjectForDatum (pnt1)
        Set pnt_onCrv(i) = pnt1
    Next
    part1.Update    '以上でicd(0 To bunkatu, 2)に座標値がｾｯﾄできた
    
    Call gen_cir(0, CInt(bunkatu))
    Call MsgBox("円弧を作れなかった部分は " & ng_arcCnt & " 箇所です。")
End Sub

Sub gen_cir(ByVal i_st As Integer, ByVal i_ed As Integer)
    If (2 > i_ed - i_st) Then
        ng_arcCnt = ng_arcCnt + 1
        Exit Sub '3点ないので計算不可
    End If
    Dim i_mid As Integer, j As Integer
    i_mid = Int((i_ed + i_st) / 2)
    For j = 0 To 2
        icd1(j) = icd(i_st, j): icd2(j) = icd(i_mid, j): icd3(j) = icd(i_ed, j)
    Next
    Dim rad1 As Double, ng_flg
    rad1 = rad_3Pkeisan 'centreも確定した
    If (0 < rad1) Then
      ng_flg = 0
            '間の点がtoleを外れたらng
      For j = i_st + 1 To i_ed - 1
        If (tole < Math.Abs(rad1 - cal_leng(icd(j, 0) - centre(0) _
            , icd(j, 1) - centre(1), icd(j, 2) - centre(2)))) Then
            ng_flg = 1
            Exit For
        End If
      Next
    Else
      ng_flg = 1
    End If
    
    If (0 = ng_flg) Then
        Dim arc1 As HybridShapeCircle3Points
        Set arc1 = Fact1.AddNewCircle3Points(pnt_onCrv(i_st), pnt_onCrv(i_mid), pnt_onCrv(i_ed))
        Hb1.AppendHybridShape arc1
        part1.Update
    Else
        Call gen_cir(i_st, i_mid) '再帰呼出
        Call gen_cir(i_mid, i_ed)
    End If
End Sub

Function rad_3Pkeisan() As Double
    Dim V1(2) As Double, V2(2) As Double
    Dim n3(2) As Double, n2(2) As Double, kp(2) As Double
    Dim mid1(2) As Double, mid2(2) As Double, i As Integer
    
    For i = 0 To 2
        V1(i) = icd2(i) - icd1(i)
        V2(i) = icd3(i) - icd2(i)
        mid1(i) = (icd2(i) + icd1(i)) / 2#
        mid2(i) = (icd2(i) + icd3(i)) / 2#
    Next
    Call Unit_Vec(V1): Call Unit_Vec(V2)
    Call gaiseki(V1, V2, n3): Call gaiseki(V2, n3, n2)
    Call Unit_Vec(n2)
    Dim nv  As Double, na As Double, d As Double, t As Double
    nv = V1(0) * n2(0) + V1(1) * n2(1) + V1(2) * n2(2)
    If 0.000000000001 > Math.Abs(nv) Then
        rad_3Pkeisan = -1
    Else
        d = V1(0) * mid1(0) + V1(1) * mid1(1) + V1(2) * mid1(2)
        na = V1(0) * mid2(0) + V1(1) * mid2(1) + V1(2) * mid2(2)
        t = (d - na) / nv
        
        For i = 0 To 2: kp(i) = mid2(i) + t * n2(i): centre(i) = kp(i): Next
        rad_3Pkeisan = cal_leng(icd1(0) - kp(0), icd1(1) - kp(1), icd1(2) - kp(2))
    End If
End Function

Sub gaiseki(vec1() As Double, vec2() As Double, resVec() As Double)
    resVec(0) = vec1(1) * vec2(2) - vec1(2) * vec2(1)
    resVec(1) = vec1(2) * vec2(0) - vec1(0) * vec2(2)
    resVec(2) = vec1(0) * vec2(1) - vec1(1) * vec2(0)
End Sub

Sub Unit_Vec(ve1() As Double)
    Dim leng1  As Double, i
    leng1 = cal_leng(ve1(0), ve1(1), ve1(2))
    For i = 0 To 2: ve1(i) = ve1(i) / leng1: Next
End Sub

Function cal_leng(X As Double, Y As Double, z As Double)
    cal_leng = Math.Sqr(X * X + Y * Y + z * z)
End Function

当時、これらの処理が全く理解できませんでした。
(今でも数学的な部分の一部が、理解できていません・・・)

自分が理解した範囲で簡単に説明しておきます。
・まず、該当する曲線上に比率による点を多数作成します。
f:id:kandennti:20160615122441p:plain

・始点(左側)と終点、中間点を元に、3点通過円弧を演算によって
　算出します。
f:id:kandennti:20160615122448p:plain

・円弧と、3点通過円弧の始点-終点間の全ての点の距離が
　ﾄﾚﾗﾝｽ以内に収まっているかどうかを評価します。
　全ての点がﾄﾚﾗﾝｽ以内であれば、円弧は確定します。

・ﾄﾚﾗﾝｽ以内に収まっていない点が存在する場合は、
　先程の中間点を終点とし同様の処理し、残りの部分は中間点を
　始点とし同様の処理を繰り返し(再帰)て行きます。
f:id:kandennti:20160615122453p:plain
　これは二分探索木と呼ばれる手法ですね。

二分探索木 - Wikipedia

続いて、ここなさんが「このような処理の方が良いのでは？」と
提案されたｺｰﾄﾞがこちらです。

'vbs - Type2
Dim icd(200, 2), pnt_onCrv(200)
Dim icd1(2), icd2(2), icd3(2)
Dim centre(2), ng_arcCnt
Dim part1, Fact1, Hb1
Const bunkatu = 150
Const tole = 0.02

Sub CATMain() 'catVBA
    ng_arcCnt = 0
    
    Dim Doc, aSel, ret1, crv1
    Set Doc = CATIA.ActiveDocument
    Set aSel = Doc.selection
    Set part1 = Doc.Part
    Set Fact1 = part1.HybridShapeFactory
    Set Hb1 = part1.HybridBodies.add
    Dim InputObjectType(0)
    InputObjectType(0) = "HybridShape"
    aSel.Clear
    ret1 = aSel.SelectElement2(InputObjectType, "sel CRV", True)
    If (ret1 <> "Normal") Then Exit Sub
    Set crv1 = aSel.Item(1).Value
    aSel.Clear

    Dim i As Integer, ratio As Double, j As Integer
    Dim rad1 As Double
    Dim pnt1, origin(2)
    For i = 0 To bunkatu
        ratio = i / bunkatu
        Set pnt1 = Fact1.AddNewPointOnCurveFromPercent(crv1, ratio, False)
        Hb1.AppendHybridShape pnt1
        part1.Update
        Call pnt1.GetCoordinates(origin)
        For j = 0 To 2: icd(i, j) = origin(j): Next
        'Fact1.DeleteObjectForDatum (pnt1)
        Set pnt_onCrv(i) = pnt1
    Next
    part1.Update    '以上でicd(0 To bunkatu, 2)に座標値がｾｯﾄできた
    
    Call gen_cir(0, CInt(bunkatu))
    Call MsgBox("円弧を作れなかった部分は " & ng_arcCnt & " 箇所です。")
End Sub

Sub gen_cir(ByVal i_st As Integer, ByVal i_ed As Integer)
    If (2 > i_ed - i_st) Then
        ng_arcCnt = ng_arcCnt + 1
        Exit Sub '3点ないので計算不可
    End If
    Dim i_mid As Integer, j As Integer
    i_mid = Int((i_ed + i_st) / 2)
    For j = 0 To 2
        icd1(j) = icd(i_st, j): icd2(j) = icd(i_mid, j): icd3(j) = icd(i_ed, j)
    Next
    Dim rad1 As Double, ng_flg
    rad1 = rad_3Pkeisan 'centreも確定した
    If (0 < rad1) Then
      ng_flg = 0
            '間の点がtoleを外れたらng
      For j = i_st + 1 To i_ed - 1
        If (tole < Math.Abs(rad1 - cal_leng(icd(j, 0) - centre(0) _
            , icd(j, 1) - centre(1), icd(j, 2) - centre(2)))) Then
            ng_flg = 1
            Exit For
        End If
      Next
    Else
      ng_flg = 1
    End If
    
    If (0 = ng_flg) Then
        Dim arc1 As HybridShapeCircle3Points
        Set arc1 = Fact1.AddNewCircle3Points(pnt_onCrv(i_st), pnt_onCrv(i_mid), pnt_onCrv(i_ed))
        Hb1.AppendHybridShape arc1
        part1.Update
            i_st = i_ed
        Call gen_cir(i_st, bunkatu) '再帰呼出
    Else
        i_ed = i_ed - 1
        Call gen_cir(i_st, i_ed) '再帰呼出
    End If
End Sub

'以下の部分は Type1と同じです

変数を定数化したり、ｽｺｰﾌﾟを変更したりしている部分もありますが、
gen_cir関数の最後の方の、再帰の呼び出し方を変更しています。

・円弧の評価する部分までは、Type1と同様です。

・ﾄﾚﾗﾝｽを外れている場合、Type1は中間点を終点とし再帰させてい
　ますが、こちらは終点の1個手前を終点とし再帰させています。
f:id:kandennti:20160615122502p:plain

確か「こちらの方法の方が、生成される円弧の数が少なくなる
可能性がある」と言った内容だったような記憶があります。

C#ATIA

↑ﾀｲﾄﾙ詐欺　主にFusion360API　偶にCATIA V5 VBA(絶賛ネタ切れ中)

曲線と戦ってみる1