Shade3D 公式

30 曲面に沿った掃引 [ Shade Labo ]


#1


**16 / 07 / 01**     一部のロジックを修正

     修正ロジックについては 28 水平掃引 [ fig 28 - 11 ] を参照
     script の修正箇所については 28 水平掃引28 - 6 Script 修正 の項を参照


関連記事 :
  • 26 滑らかな掃引(1/2)
  • 27 滑らかな掃引(2/2)
  • 28 水平掃引

30 - 1  カスタマイズされた掃引


**28 水平掃引** では **c_sweep( )** をベースにして、断面配置のルールを変更した掃引体を作成しています。

大まかな処理の流れをまとめてみます。

  1. 断面配置のルールに応じて掃引中心線の 分割ルール を定め分割

  2. 一定のルールの下に 断面を配置

  3. 交差方向 handle 長さの基準となる面法線 Vp を求め

  4. Vp を用いて交差方向 handle 長さ を定める

       断面配置で ねじれ を与えた場合、handle 方向の求め方は2つに分かれる

       直線 / 平坦性 部分を検出し、そこでは handle 方向を中心線のそれに合わせる

       それ以外では、隣接断面ポイント座標から求める



**c_sweep( )** による滑らかな掃引と、**h_sweep( )** による水平掃引の処理の本質的な違いは 2)の **断面配置ルール** です。

そして、それに従って 分割ルール直線 / 平坦性 の定義 が決定されます。


このことから、さらなるカスタマイズの可能性が見えてきます。

すなわち、任意の断面配置ルールを定め、それに応じて 分割ルールと直線 / 平坦性の定義をモディファイすれば、様々な掃引体を作ることが可能です。


ここでは、その一つとして 曲面に沿った掃引 を考えてみます。


#2

30 - 2  仕様


**曲面に沿った掃引 s_sweep( )** の仕様を以下に記します。
**< 機能 >**

自由曲面内の線形状を掃引中心線とし、中心線の接線ベクトルと自由曲面の面法線を基準 にして断面配置を行う

     

     [ fig 30 - 1 ]



**< 対象 >**

自由曲面内の線形状 を掃引中心線とする

面法線は指定線形状の 交差方向下流 の bezier patch から求める

ただし、交差方向が開いている自由曲面の終端の線形状が対象の場合は、上流の patch から求める

中心線に U 方向一点接触重点や重線があっても構わない

V 方向( 交差方向 )については、極のような一点接触重点は構わないが、重線は不可

     

     [ fig 30 - 2 ]

          [ fig 30 - 3 ]



**< 中心線分割 >**

中心線は 曲がり( 全曲率 )ベースによるものと、ねじれ ベースによるものを比較し、よりクリティカルな方を基準にして分割

          [ fig 30 - 4 ]



**< 平坦性の定義 >**

掃引中心線を基準とした基準面にねじれがない = 平坦 とみなす

          [ fig 30 - 5 ]


#3

30 - 3  断面配置


関連記事:
  • 21 Quaternion

断面配置は **21 Quaternion** で紹介した関数 **attitude_control_quaternion( )** を用います。
def **attitude_control_quaternion**( u1 as vec3, v1 as vec3, u2 as vec3, v2 as vec3, fit as integer = 0 ) as quaternion

          [ fig 30 - 6 ]


**s_sweep( )** の関数 **s_sweep_section_quaternion( )** を [ script 30 - 1 ] に記します。
[ script 30 - 1 ]
#  s_sweep_section_quaternion(bz as vec3 list, bp as vec3 matrix) as quaternion
#
#	bz, bp から掃引時の断面形状の姿勢を与える quaternion を返す
#	vn1, vn2 :	当該区間の 始端 / 終端 での面法線

def s_sweep_section_quaternion(bz, bp) :	
	u1 = labo.bezier_line_tangent(bz, 0)					#  outhandle 側接線ベクトル

	if u1.abs2() == 0 :
		return quaternion()									#  重点 ( 面積のない patch も重点と見なす )
	else :
		v1 = labo.bezier_surface_v_tangent_2(bp, 0., 0.)	#  始端での V 方向接線ベクトル
#		if (v1 == None) or (v1.abs2 == 0) :					#  V 方向につぶれた patch		( 注記-2 参照 )
#			return None										#  サポートされない自由曲面	( 注記-2 参照 )
			
	u2 = labo.bezier_line_tangent(bz, 1)					#  inhandle 側接線ベクトル
	v2 = labo.bezier_surface_v_tangent_2(bp, 0., 1.)		#   終端での V 方向接線ベクトル
#	if (v2 == None) or (v2.abs2 == 0) :						#  V 方向につぶれた patch		( 注記-2 参照 )
#		return None											#  サポートされない自由曲面	( 注記-2 参照 )
	
	vn1 = u1*v1			#  注記-1 参照
	vn2 = u2*v2			#  注記-1 参照
	return labo.attitude_control_quaternion(u1, vn1, u2, vn2, 0)	
		
#	注記-1:
#		vn1, vn2 を直接 bezier_surface_normal( ) あるいは bezier_surface_normal_2( ) で求めてはダメ
#
#	注記-2:
#		V 方向につぶれた patch を検出した場合、サポートされない自由曲面として None を返すが、
#		既に実行した関数 bezier_divide_by_total_curvature_or_torsion( ) でチェック済みであるので不要

#4

30 - 4  曲面の 曲がり と ねじれ


一般に曲面に対する **ねじれ** を表現する計量はなく、**曲がり** の計量である面曲率はここでは用いません。

曲面上の曲線を中心とした local な領域において 曲がりねじれ を次のように定義します。


  • 曲面上を走るある曲線を中心にして、曲面の曲線近傍における局所的な微小幅の面を取り出す

  • 曲線の接線方向を Vt、面の面法線を Vn とすれば、 Vt で前後方向、Vn で上下方向、Vt x Vn で左右方向をそれぞれ定める

  • この時、 pitching と yawing を合成したものを 曲がり、rolling を ねじれ と定義する



**< Bezier 曲面上の Bezier 曲線を中心とした 局所的な 曲がり と ねじれ >**
1)Bezier 曲面 bp の U 方向線形状 Lu を中心線とする
          [ fig 30 - 7 ]
2)中心線 Lu 上の点 **P1** に対し、接線ベクトルと面法線ベクトルを **U1**, **Vn1** とする
3)点 **P1** から中心線上を微小距離 Δs 進んだ点 **P2** での接線ベクトルと面法線ベクトルを **U2**, **Vn2** とし、面法線ベクトル **Vn1** を起点 **P1** から **P2** に曲がりに沿って回転移動したものを **W** とする
4)このとき、微小区間での面の **曲がり** と **ねじれ** を次のように表す

       面の 曲がり = P1 - P2 間での中心線の曲がり = U1U2 の傾き

       面の ねじれ = WVn2 との傾き


5)プログラムでは、連続する各微小区間での 曲がり, ねじれ それぞれの和で、当該区間全体での 曲がり量, ねじれ量 とみなす

     

     [ fig 30 - 8 ]


#5

30 - 5  中心線分割


中心線は **曲がり** ベースによるものと、**ねじれ** ベースによるものを比較し、よりクリティカルな方を基準にして分割します。

大まかな処理の流れは次のようになります。


1)当該 bezier 区間の 屈曲点( 変曲点 )を求め、その点で分割


2)1)で分割された bezier 区間を parameter 基準で 8分割 し、曲がり量( 全曲率 )と ねじれ量を 概算
       指定するそれぞれのしきい値と比較して、どちらがクリティカルであるかを判定


3)2)で判定された クリティカルなベース で bezier 区間を最適分割して 曲がり量 / ねじれ量 を計算


4)当該 Bezier 区間の 曲がり量 / ねじれ量 から 分割数 N を決定


5)曲がり量 ベースでの分割では曲がりの N 分割を与える bezier parameter を線形補間で求め 分割

       ねじれ量 ベースでの分割では 単純な bezier parameter の N 分割で分割



**5)**の処理で、ねじれ量がクリティカルな場合に、ねじれ量ベースで N 分割するのではなく、**単純な bezier parameter の N 分割** にしているのは以下の理由によります。
曲がりの大きな所とねじれの大きな所は一致しないことが多く、ねじれ量基準で分割すると曲がり分割の不均一さが目立ち、形状全体の滑らかさに障害となってしまいます。

また、基準となる自由曲面のレンダリング自体が bezier parameter の等分割( default 8 分割 )でなされており、これはある程度の曲がり具合の等分割に近い分割を与えますので、自由曲面と合わせたレンダリング結果の親和性の上でも、parameter の N 分割の方が好ましい結果が得られます。

     

     [ fig 30 - 9 ]



[ script 30 - 2 ] は中心線分割に係わる部分を抜粋したものです。

関数 bezier_divide_by_total_curvature_or_torsion( ) によって、上記 2)~ 5)の処理を行い、分割 parameter list を返します。

なお、上で述べた ねじれ量ベースで N 分割 するコードもコメントアウトした状態で記述してあります。


[ script 30 - 2 ]

#  bezier_divide_by_total_curvature_or_torsion( bz as vec3 list, bp, as vec3 matrix, threshold1 as float, threshold2 as float) as float list
#
#		bz : 		bezier control point 座標を格納した vec3 list
#		bp : 		bezier patch
#		threshold1 :	分割しきい値 ( 全曲率ベース )
#		threshold2 :	分割しきい値 ( ねじれベース )
#
#		bezier line bz を分割する bezier parameter を返す。
#		当該区間の全曲率あるいは bP のねじれの変化量の内、大なる方を分割基準とする
#		分割がない場合は、空 list が返される	
#		面法線が求まらない = ねじれが求まらない 場合は、None が返される

def bezier_divide_by_total_curvature_or_torsion(bz, bp, threshold1, threshold2) :
	from math import acos, pi

	paramL = []								#  分割点の bezier parameter を格納する list
	
	u1 = labo.bezier_line_tangent(bz, 0)
	if u1.abs2() == 0 :
		return paramL						#  重点による長さのない区間
	else :
		uu1 = u1
		
	u2 = labo.bezier_line_tangent(bz, 1)
	
	a = u1.dot(u2)
	if abs(a) > 0.99999994 :
		w = bz[3] - bz[0]
		if w.abs2() == 0 :
			return paramL					#  重点とみなす
		elif a < - 0.99999994 :
			w.norm()
			if abs(w.dot(u1)) > 0.99999994 :		#  判定 - 1
				return paramL				#  重点と同じ扱いにする  
		
	v1 = labo.bezier_surface_v_tangent_2(bp, 0., 0.)
	if (v1 == None) or (v1.abs2() == 0) :	
		return None							#  patch V 方向がつぶれている
	else :
		vn1 = u1*v1
		vn1.norm()
		vvn1 = vn1	
	

	#  8 分割で全曲率 theta1 と ねじれ theta2 を求める
	theta1, theta2 = 0, 0
	for i in range(1, 9) :
		t = i/8.	
		u2 = labo.bezier_line_tangent(bz, t)
		v2 = labo.bezier_surface_v_tangent_2(bp, 0., t)	
		if (v2 == None) or (v2.abs2() == 0) :	
			return None						#  patch V 方向がつぶれている あるいは V 方向自己交差
			
		cosT1 = max(-1, min(1, u1.dot(u2)))
		theta1 += acos(cosT1)				#  全曲率
		
		vn2 = u2*v2
		vn2.norm()
		Q = quaternion().slerp(u1, u2)		#  上記の 判定 - 1 より、Q != None
		w = Q*vn1	
		
		cosT2 = max(-1, min(1, w.dot(vn2)))
		theta2 += abs(acos(cosT2))			#  ねじれ
		
		u1 = u2
		vn1 = vn2
		

	#  theta1/threshold1 と theta2/threshold2 の比較から、分割基準のベースを選択して分割
	pL = [0.]								#  各区間を分割する bezer parameter を格納する list
	cL = [0.]								#  各区間までの全曲率を格納する list
	a1 = theta1/threshold1
	a2 = theta2/threshold2

	if a1 > 0.2 or a2 > 0.2 :
		if a1 > a2 :
			make_total_curvature_list(bz, 0., uu1, 1., u2, pL, cL, 8)	#  最適分割を行って cL[], pL[] にデータを格納する 再帰関数
			theta = cL[len(cL) - 1]									#  全曲率
			if theta > threshold1 :									#  当該 bezier 区間の全曲率がしきい値を越えている
				n = int((theta - 0.00001)/threshold1) + 1			#  分割数 ( 丸め誤差による分割数の不均一防止のため 0.00001 を減ずる )
				dc = theta/n				
				paramL = [0.]										#  分割点の bezier parameter を格納 (最初は 0. )
				j = 1
				for i in range(1, n) :
					c = i*dc
					while c > cL[j] :
						j += 1
					paramL.append(pL[j] - (pL[j] - pL[j - 1])*(cL[j] - c)/(cL[j ] - cL[j - 1]))	#  分割点の bezier parameter
				paramL.append(1.)													#  分割点の bezier parameter list の最後に 1 を追加
			
			
		else :				
			make_torsion_list(bp, 0., uu1, vvn1, 1., u2, vn2, pL, cL, 8)		#  最適分割を行って cL[], pL[] にデータを格納する 再帰関数
			theta = cL[len(cL) - 1]								#  ねじれ
		
			if theta > threshold2 :								#  当該 bezier 区間のねじれ量がしきい値を越えている
				n = int((theta - 0.00001)/threshold2) + 1		#  分割数 ( 丸め誤差による分割数の不均一防止のため 0.00001 を減ずる )
			
				#   純粋にねじれ量をベースにして分割すると、曲がりに対する分割がうまくいかない 
				#  ベースとなる自由曲面との親和性も項考慮すれば、両者を折衷して、単純な parameter ベースでの分割を採用する
			
#				dc = theta/n				
#				paramL = [0.]									#  分割点の bezier parameter を格納 (最初は 0. )
#				j = 1
#				for i in range(1, n) :
#					c = i*dc
#					while c > cL[j] :
#						j += 1
#					paramL.append(pL[j] - (pL[j] - pL[j - 1])*(cL[j] - c)/(cL[j ] - cL[j - 1]))	#  分割点の bezier parameter
#				paramL.append(1.)	
			
				paramL = [float(i)/n for i in range(n + 1)]		#  分割点の bezier parameter を格納 (最初は 0. 最後は 1.)
		
	return paramL
	
	
	

#  make_total_curvature_list(bz as vec3 list, t1 as float, v1 as vec3, t2 as float, v2 as vec3, pL as float list, cL as float list, n as int)
#
#		bz の parameter 区間 t1~t2 を最適サイズに分割して、cL[], pL[] にデータを格納する再帰関数
#
#		bz : 		bezier control point 座標を格納した vec3 list
#		t1, v1 :		計算対象区間の 始端 bezier parameter と その点での接線ベクトル	
#		t2, v2 :		計算対象区間の 終端 bezier parameter と その点での接線ベクトル
#		pL :			各区間を分割する bezer parameter を格納
#		cL :			分割された各区間までの全曲率を格納
#		n :			分割数 ( bezier parameter t1~t2 を n 分割する )

def make_total_curvature_list(bz, t1, v1, t2, v2, pL, cL, n) :
	from math import acos
	
	u = v1
	s1 = t1
	ds = (t2 - t1)/n								#  parameter  増分
	
	for i in range(1, n + 1) :						#  分割数 n で分割
		s2 = s1 + ds

		if i < n :
			v = labo.bezier_line_tangent(bz, s2)	# parameter t での接線ベクトル
		else :
			v = v2
			
		cosT = max(-1, min(1, u.dot(v)))
		
		if cosT >= 0.99966 :						#  0.99966 = cos(1.5度)
			pL.append(s2)
			cL.append(cL[len(cL) - 1] + acos(cosT))	
		else :
			theta =  acos(cosT)
			m = max(2, 1 + int(theta*180/pi/1.5))						#  分割数
			make_total_curvature_list(bz, s1, u, s2, v, pL, cL, m)		#  再分割
			
		u = v
		s1 = s2
		

		
		
#  make_torsion_list(bp as vec3 list, t1 as float, uu1 as vec3, vvn1 as bec3, t2 as float, uu2 as vec3, vvn2 as vec3, pL as float list, tL as float list, n as int)
#
#		bp の 区間 t1~t2 を最適サイズに分割して、cL[], pL[] にデータを格納する再帰関数
#
#		bp : 			bezier patch bp
#		t1, u1, vn1 :		計算対象区間の 始端 bezier parameter と その点での接線ベクトルと面法線	
#		t2, u2, vm2 :	計算対象区間の 終端 bezier parameter と その点での接線ベクトルと面法線
#		pL :				各区間を分割する bezer parameter を格納
#		tL :				分割された各区間までの全曲がり量を格納
#		n :				分割数 ( bezier parameter t1~t2 を n 分割する )

def make_torsion_list(bp, t1, uu1, vvn1, t2, uu2, vvn2, pL, tL, n) :
	from math import acos, pi

	u1 = uu1
	vn1 = vvn1
	s1 = t1
	ds = (t2 - t1)/n										#  parameter  増分

	for i in range(1, n + 1) :								#  分割数 n で分割
		s2 = s1 + ds
		if i < n :
			u2 = labo.bezier_surface_u_tangent(bp, 0., s2)
			v2 = labo.bezier_surface_v_tangent_2(bp, 0., s2)
			vn2 = u2*v2
			vn2.norm()
												
		else :
			u2 = uu2
			vn2 = vvn2
		
		Q = quaternion().slerp(u1, u2)
		w = Q*vn1
		cosT = max(-1, min(1, w.dot(vn2)))
				
		#  捻れ量ベースの N 分割は求めないので、曲がりの時よりも最小角を大きめにとる
		if cosT >= 0.998629 :								#  0.998629 = cos(3度)	
			pL.append(s2)
			tL.append(tL[len(tL) - 1] + acos(cosT))	
		else :
			theta = acos(cosT)
			m = max(2, 1 + int(theta*180/pi/3))				#  分割数
			make_torsion_list(bp, s1, u1, vn1, s2, u2, vn2, pL, tL, m)	#  再分割
	
		u1 = u2	
		vn1 = vn2
		s1 = s2

#6

30 - 6  平坦性


平坦性判定は **非ねじれ判定** になります。

平非ねじれ判定リストを返す関数 s_flatness_list( ) を [ script 30 - 3 ] に記します。


[ script 30 - 3 ]
#  s_flatness_list(bZs as vec3 list, bPs as vec3 matrix list, closed as bool) as bool list
#		bZs :	掃引中心線の control point 座標 list
#		bPs :	掃引中心線を含む bezier patch list
#		closed :	中心線の開閉情報
#
#		bPs の各 patch 区間の 非ねじれ list を返す
#
#		list の始端と終端には次なる項目を追加
#			中心線が閉じている場合、その隣接項目
#			中心線が開いている場合、False

def s_flatness_list(bZs, bPs, closed):
	flatL = []
	n = len(bPs)
	for i in range(n) :
		bz = bZs[i]
		bp = bPs[i]
		
		u1 = labo.bezier_line_tangent(bz, 0)	#  outhandle 側接線ベクトル
		u2 = labo.bezier_line_tangent(bz, 0.5)	#  中点での接線ベクトル
		
		if u1.abs2() == 0 :						#  重点
			flatL.append(False)
		
		else :
			v1 = labo.bezier_surface_v_tangent_2(bp, 0., 0.)	#  ここの関数が呼ばれる時点で、v1, v2  が求まることは既に担保されている
			v2 = labo.bezier_surface_v_tangent_2(bp, 0., 0.5)
			vn1 = u1*v1
			vn1.norm()
			vn2 = u2*v2
			vn2.norm()
			Q = quaternion().slerp(u1, u2)
			w = Q*vn1
			flatL.append(w.dot(vn2) >= 0.99985)	#   本来なら しきい値は 0.9999 とするが、ベースとなる自由曲面の持つ理想形状との誤差を考慮し、若干緩めに設定
		
	#   関数 h_sweep_handle( ) に与える引数の選択において、平坦性 list fL の隣接二値を比較するため、 list の先頭と終端に値を追加
	if not closed :
		flatL.insert(0, True)
		flatL.append(True)
	else :
		flatL.insert(0, flatL[len(flatL) - 1])
		flatL.append(flatL[1])
		
	return flatL

#7

30 - 7  Script


以上から水平掃引の関数 **h_sweep( )** をベースにして、曲面に沿った掃引を作成する機能に改編した関数 **s_sweep( )** を作成、しきい角度などを指定するユーザーインターフェースを追加して [ script 30 - 4 ] として下に記します。
**< 注記 >**

script 内で bezier patch list bPs( 選択された掃引中心線形状を含む自由曲面 ) と bezier list bZs( 掃引中心線 )の二つを併記しているが、これらには

     bPs[k][0] = bZs[k]

なる関係があり、bPs のみの handling にすることも可能


ここでは、U 方向( 中心線 )に関する処理 と V 方向( 自由曲面交差方向 )に関する処理を識別しやすくするために、あえて双方を併記し使い分けている

u = labo.bezier_line_tangent(bz, 1)

v = labo.bezier_surface_v_tangent_2(bp, 0., 1.)


[ script 30 - 4 ]
**16 / 07 / 01**   ロジック修正に伴い 関数 **get_vt( )** を変更

     変更箇所については 28 水平掃引 [ script 28 - 6 ] を参照


h_sweep( ) からの変更や追加部分には行末に ############ を付しています

掃引中心線としての自由曲面内の線形状 と 掃引断面 を選択して実行します。

import labo
from matrix import *
from quaternion import *

		

##  全面変更	  ########################################################################
		
#  s_sweep_section_quaternion(bz as vec3 list, bp as vec3 matrix) as quaternion
#
#	bz, bp から掃引時の断面形状の姿勢を与える quaternion を返す
#	vn1, vn2 :	当該区間の 始端 / 終端 での面法線

def s_sweep_section_quaternion(bz, bp) :	
	u1 = labo.bezier_line_tangent(bz, 0)					#  outhandle 側接線ベクトル

	if u1.abs2() == 0 :
		return quaternion()									#  重点 ( 面積のない patch も重点と見なす )
	else :
		v1 = labo.bezier_surface_v_tangent_2(bp, 0., 0.)	#  始端での V 方向接線ベクトル
#		if (v1 == None) or (v1.abs2 == 0) :					#  V 方向につぶれた patch		( 注記-2 参照 )
#			return None										#  サポートされない自由曲面	( 注記-2 参照 )
			
	u2 = labo.bezier_line_tangent(bz, 1)					#  inhandle 側接線ベクトル
	v2 = labo.bezier_surface_v_tangent_2(bp, 0., 1.)		#   終端での V 方向接線ベクトル
#	if (v2 == None) or (v2.abs2 == 0) :						#  V 方向につぶれた patch		( 注記-2 参照 )
#		return None											#  サポートされない自由曲面	( 注記-2 参照 )
	
	vn1 = u1*v1			#  注記-1 参照
	vn2 = u2*v2			#  注記-1 参照
	return labo.attitude_control_quaternion(u1, vn1, u2, vn2, 0)	
		
#	注記-1:
#		vn1, vn2 を直接 bezier_surface_normal( ) あるいは bezier_surface_normal_2( ) で求めてはダメ
#
#	注記-2:
#		V 方向につぶれた patch を検出した場合、サポートされない自由曲面として None を返すが、
#		既に実行した関数 bezier_divide_by_total_curvature_or_torsion( ) でチェック済みであるので不要

##################################################################################			
		
			
				
#  sweep_plane_normal(p0 as vec3, p1 as vec3, p2 as vec3, p3 as vec3) as vec3
#
#	掃引中心線の bezier data bz (p0, p1, p2) から、掃引体交差方向の handle 長さを求めるための面法線を返す
#	基準となる anchor point に handle  が出ていることが前提
#
#		p0, p1, p2, p3 :	outhandle 側を求める場合は	bz[0], bz[1], bz[2], bz[3]
#							inhandle 側を求める場合は		bz[3], bz[2], bz[1], bz[0]

def sweep_plane_normal(p0, p1, p2, p3) :
	from math import pi, acos, sin, atan
	
	v1 = p1 - p0					#  v1.abs() != 0 が前提
	v2 = p2 - p1
	v1.norm()
	v2.norm()
	
	if v2.abs2() == 0 :
		v3 = p3 - p1
		v3.norm()
		if v3.abs2() == 0 :
			return v1
		else :
			vn = v1 + v3
			vn.norm()
			vp = v1 + vn
			vp.norm()
			return vp
		
	else :
		a = v1.dot(v2)
		if abs(a) >= 0.99999994 :
			return v1
		
		elif a >= 0 :
			vn = v1 + v2
			vn.norm()
			
			sinA = a
			cosA = (1 - sinA**2)**0.5
			beta = (pi - acos(a))/2
			r = (1 - sinA)/3
			
			axis = v1*v2
			axis.norm()
			v = v1 - v2
			v.norm()
			u = v/sin(beta) + v1*r/cosA
			u.norm()
			phi = acos(u.dot(v))

		else :
			axis = v1*v2
			phi = atan(4./3) - pi/4
			Q = quaternion().rotation(axis, pi/4)
			vn = Q*v1

		v3 = p3 - p2
		v3.norm()
		if v3.abs2() == 0 :
			return vn
		else :
			d = v1.dot(v2)
			w = -(v1-2*d*v2)
			w.norm()
			e = v3.dot(w)
			Q = quaternion().rotation(axis, e*phi)
			vp = Q*vn
			vp.norm()
			return vp

	
	

#  h_sweep_handle(w as vec3, vr as vec3,  Q as vec3, r as float, p as vec3, vp as vec3) as vec3
#
#	handle 座標を返す

def h_sweep_handle(w, vr, r, p, vp) :	
	#     L :  点 w を通り、vr 方向に進む直線
	#	PL1 : 点 p を通り面法線 vp なる面
	#	精度確保のため r = 1 とはせず、r = bounding box size / 16 としてある

	d1 = vp.dot(w - p)			#  点 w と 面 PL1 との符号付き距離
	d2 = vp.dot(w + r*vr - p)	#  点 w + r*vr と 面 PL1 との符号付き距離
	d = r*d1/(d1 - d2)			#  点 w を通り vr 方向に進む直線と、面 PL1 との交点までの 符号付き距離
	if d < 0 :					#  食い込みが発生する
		d = -d/12
	return w + d*vr
	

		
#  bezier_inflection_point( bz as vec3 list) as float list
#
#		bz : bezier control point 座標を格納した vec3 list
#		bezier line bz の屈曲点位置の bezier parameter を返す。
#		屈曲点がない場合は、空 list が返される

def bezier_inflection_point(bz) :
	w1 = (bz[1] - bz[0])*(bz[2] - bz[3])
	w2 = (bz[2] - bz[1])*(bz[3] - bz[0])
	w3 = (bz[2] - bz[0])*(bz[1] - bz[3])
	r1 = w1.norm()
	r2 = w2.norm()
	r3 = w3.norm()
	L = [[r1, w1], [r2, w2], [r3, w3]]
	L.sort()
	vn = L[2][1]
	if vn.norm() == 0 :				#  bz は直線であり、変曲点は存在しない
		return []

	#  bz2 : 原点を通り面法線 vn なる面に bz を投影	
	bz2 = []
	for i in range(4) :
		d = vn.dot(bz[i])
		bz2.append(bz[i] - d*vn)
	
	#  Q : bz2 の階差表現
	bz2 = matrix(False, bz2)
	Md = labo.bezier_d()
	Q = bz2*Md
	
	#  こちらの方法でも可
#	bz2 = matrix(True, bz2)
#	Md = labo.bezier_d()
#	Md.transpose()
#	Q = Md*bz2			

	A = Q[2]*Q[3]
	B = Q[1]*Q[3]
	C = Q[1]*Q[2]
	
	bb = labo.vec3_bounding_box_size(Q)
	bb.pop()
	bb2 = sorted(bb)
	for i in range(3) :
		if bb2[0] == bb[i] :
			a = A[i]
			b = B[i]
			c = C[i]
			t = []
			if a != 0 :
				d = b**2 - 4*a*c
				if d > 0 :
					tt = (-b - d**0.5)/(2*a)
					if tt > 0 and tt < 1 :
						t.append(tt)
					tt = (-b + d**0.5)/(2*a)
					if tt >= 0 and tt <= 1 :
						t.append(tt)
				elif d == 0 :
					tt = -b/(2*a)
					if tt > 0 and tt < 1 :
						t.append(tt)
			elif b != 0 :
				tt = -c/b
				if tt > 0 and tt < 1 :
					t.append(tt)
			
			if len(t) == 1 :
				v1 = labo.bezier_line_tangent(bz, 0)
				v2 = labo.bezier_line_tangent(bz, 1)
				v0 = labo.bezier_line_tangent(bz, t[0])
				if v0.dot(v1) >= 0.9999 or  v0.dot(v2) >= 0.9999 :	#  屈曲点での接線ベクトルが、始端 / 終端での接線ベクトルと
					t.pop()											#  ほとんど同じであれば、その屈曲点では分割しない
					
			elif len(t) == 2 :
				t.sort()	
				v1 = labo.bezier_line_tangent(bz, 0)
				v2 = labo.bezier_line_tangent(bz, 1)
				v0 = labo.bezier_line_tangent(bz, t[1])
				if v0.dot(v2) >= 0.9999 :							#  屈曲点での接線ベクトルが、終端での接線ベクトルと
					del t[1]										#  ほとんど同じであれば、その屈曲点では分割しない
				v0 = labo.bezier_line_tangent(bz, t[0])
				if v0.dot(v1) >= 0.9999 :							#  屈曲点での接線ベクトルが、始端での接線ベクトルと
					del t[0]										#  ほとんど同じであれば、その屈曲点では分割しない
			return t
				
	

##  全面変更	  ########################################################################

#  bezier_divide_by_total_curvature_or_torsion( bz as vec3 list, bp, as vec3 matrix, threshold1 as float, threshold2 as float) as float list
#
#		bz : 		bezier control point 座標を格納した vec3 list
#		bp : 		bezier patch
#		threshold1 :	分割しきい値 ( 全曲率ベース )
#		threshold2 :	分割しきい値 ( ねじれベース )
#
#		bezier line bz を分割する bezier parameter を返す。
#		当該区間の全曲率あるいは bP のねじれの変化量の内、大なる方を分割基準とする
#		分割がない場合は、空 list が返される	
#		面法線が求まらない = ねじれが求まらない 場合は、None が返される

def bezier_divide_by_total_curvature_or_torsion(bz, bp, threshold1, threshold2) :
	from math import acos, pi

	paramL = []								#  分割点の bezier parameter を格納する list
	
	u1 = labo.bezier_line_tangent(bz, 0)
	if u1.abs2() == 0 :
		return paramL						#  重点による長さのない区間
	else :
		uu1 = u1
		
	u2 = labo.bezier_line_tangent(bz, 1)
	
	a = u1.dot(u2)
	if abs(a) > 0.99999994 :
		w = bz[3] - bz[0]
		if w.abs2() == 0 :
			return paramL					#  重点とみなす
		elif a < - 0.99999994 :
			w.norm()
			if abs(w.dot(u1)) > 0.99999994 :		#  判定 - 1
				return paramL				#  重点と同じ扱いにする  
		
	v1 = labo.bezier_surface_v_tangent_2(bp, 0., 0.)
	if (v1 == None) or (v1.abs2() == 0) :	
		return None							#  patch V 方向がつぶれている
	else :
		vn1 = u1*v1
		vn1.norm()
		vvn1 = vn1	
	

	#  8 分割で全曲率 theta1 と ねじれ theta2 を求める
	theta1, theta2 = 0, 0
	for i in range(1, 9) :
		t = i/8.	
		u2 = labo.bezier_line_tangent(bz, t)
		v2 = labo.bezier_surface_v_tangent_2(bp, 0., t)	
		if (v2 == None) or (v2.abs2() == 0) :	
			return None						#  patch V 方向がつぶれている あるいは V 方向自己交差
			
		cosT1 = max(-1, min(1, u1.dot(u2)))
		theta1 += acos(cosT1)				#  全曲率
		
		vn2 = u2*v2
		vn2.norm()
		Q = quaternion().slerp(u1, u2)		#  上記の 判定 - 1 より、Q != None
		w = Q*vn1	
		
		cosT2 = max(-1, min(1, w.dot(vn2)))
		theta2 += abs(acos(cosT2))			#  ねじれ
		
		u1 = u2
		vn1 = vn2
		

	#  theta1/threshold1 と theta2/threshold2 の比較から、分割基準のベースを選択して分割
	pL = [0.]								#  各区間を分割する bezer parameter を格納する list
	cL = [0.]								#  各区間までの全曲率を格納する list
	a1 = theta1/threshold1
	a2 = theta2/threshold2

	if a1 > 0.2 or a2 > 0.2 :
		if a1 > a2 :
			make_total_curvature_list(bz, 0., uu1, 1., u2, pL, cL, 8)	#  最適分割を行って cL[], pL[] にデータを格納する 再帰関数
			theta = cL[len(cL) - 1]									#  全曲率
			if theta > threshold1 :									#  当該 bezier 区間の全曲率がしきい値を越えている
				n = int((theta - 0.00001)/threshold1) + 1			#  分割数 ( 丸め誤差による分割数の不均一防止のため 0.00001 を減ずる )
				dc = theta/n				
				paramL = [0.]										#  分割点の bezier parameter を格納 (最初は 0. )
				j = 1
				for i in range(1, n) :
					c = i*dc
					while c > cL[j] :
						j += 1
					paramL.append(pL[j] - (pL[j] - pL[j - 1])*(cL[j] - c)/(cL[j ] - cL[j - 1]))	#  分割点の bezier parameter
				paramL.append(1.)												#  分割点の bezier parameter list の最後に 1 を追加
			
			
		else :				
			make_torsion_list(bp, 0., uu1, vvn1, 1., u2, vn2, pL, cL, 8)		#  最適分割を行って cL[], pL[] にデータを格納する 再帰関数
			theta = cL[len(cL) - 1]								#  ねじれ
		
			if theta > threshold2 :								#  当該 bezier 区間のねじれ量がしきい値を越えている
				n = int((theta - 0.00001)/threshold2) + 1		#  分割数 ( 丸め誤差による分割数の不均一防止のため 0.00001 を減ずる )
			
				#   純粋にねじれ量をベースにして分割すると、曲がりに対する分割がうまくいかない 
				#  ベースとなる自由曲面との親和性も項考慮すれば、両者を折衷して、単純な parameter ベースでの分割を採用する
			
#				dc = theta/n				
#				paramL = [0.]									#  分割点の bezier parameter を格納 (最初は 0. )
#				j = 1
#				for i in range(1, n) :
#					c = i*dc
#					while c > cL[j] :
#						j += 1
#					paramL.append(pL[j] - (pL[j] - pL[j - 1])*(cL[j] - c)/(cL[j ] - cL[j - 1]))	#  分割点の bezier parameter
#				paramL.append(1.)	
			
				paramL = [float(i)/n for i in range(n + 1)]		#  分割点の bezier parameter を格納 (最初は 0. 最後は 1.)
		
	return paramL
	
	
	

#  make_total_curvature_list(bz as vec3 list, t1 as float, v1 as vec3, t2 as float, v2 as vec3, pL as float list, cL as float list, n as int)
#
#		bz の parameter 区間 t1~t2 を最適サイズに分割して、cL[], pL[] にデータを格納する再帰関数
#
#		bz : 		bezier control point 座標を格納した vec3 list
#		t1, v1 :		計算対象区間の 始端 bezier parameter と その点での接線ベクトル	
#		t2, v2 :		計算対象区間の 終端 bezier parameter と その点での接線ベクトル
#		pL :			各区間を分割する bezer parameter を格納
#		cL :			分割された各区間までの全曲率を格納
#		n :			分割数 ( bezier parameter t1~t2 を n 分割する )

def make_total_curvature_list(bz, t1, v1, t2, v2, pL, cL, n) :
	from math import acos
	
	u = v1
	s1 = t1
	ds = (t2 - t1)/n								#  parameter  増分
	
	for i in range(1, n + 1) :						#  分割数 n で分割
		s2 = s1 + ds

		if i < n :
			v = labo.bezier_line_tangent(bz, s2)	# parameter t での接線ベクトル
		else :
			v = v2
			
		cosT = max(-1, min(1, u.dot(v)))
		
		if cosT >= 0.99966 :						#  0.99966 = cos(1.5度)
			pL.append(s2)
			cL.append(cL[len(cL) - 1] + acos(cosT))	
		else :
			theta =  acos(cosT)
			m = max(2, 1 + int(theta*180/pi/1.5))					#  分割数
			make_total_curvature_list(bz, s1, u, s2, v, pL, cL, m)	#  再分割
			
		u = v
		s1 = s2
		

		
		
#  make_torsion_list(bp as vec3 matrix, t1 as float, uu1 as vec3, vvn1 as bec3, t2 as float, uu2 as vec3, vvn2 as vec3, pL as float list, tL as float list, n as int)
#
#		bp の 区間 t1~t2 を最適サイズに分割して、cL[], pL[] にデータを格納する再帰関数
#
#		bp : 			bezier patch bp
#		t1, u1, vn1 :	計算対象区間の 始端 bezier parameter と その点での接線ベクトルと面法線	
#		t2, u2, vm2 :	計算対象区間の 終端 bezier parameter と その点での接線ベクトルと面法線
#		pL :			各区間を分割する bezer parameter を格納
#		tL :			分割された各区間までの全曲がり量を格納
#		n :				分割数 ( bezier parameter t1~t2 を n 分割する )

def make_torsion_list(bp, t1, uu1, vvn1, t2, uu2, vvn2, pL, tL, n) :
	from math import acos, pi

	u1 = uu1
	vn1 = vvn1
	s1 = t1
	ds = (t2 - t1)/n										#  parameter  増分

	for i in range(1, n + 1) :								#  分割数 n で分割
		s2 = s1 + ds
		if i < n :
			u2 = labo.bezier_surface_u_tangent(bp, 0., s2)
			v2 = labo.bezier_surface_v_tangent_2(bp, 0., s2)
			vn2 = u2*v2
			vn2.norm()
												
		else :
			u2 = uu2
			vn2 = vvn2
		
		Q = quaternion().slerp(u1, u2)
		w = Q*vn1
		cosT = max(-1, min(1, w.dot(vn2)))
				
		#  捻れ量ベースの N 分割は求めないので、曲がりの時よりも最小角を大きめにとる
		if cosT >= 0.998629 :								#  0.998629 = cos(3度)	
			pL.append(s2)
			tL.append(tL[len(tL) - 1] + acos(cosT))	
		else :
			theta = acos(cosT)
			m = max(2, 1 + int(theta*180/pi/3))				#  分割数
			make_torsion_list(bp, s1, u1, vn1, s2, u2, vn2, pL, tL, m)	#  再分割
	
		u1 = u2	
		vn1 = vn2
		s1 = s2
		
##################################################################################


	

#  get_smooth_handle(p1 as vec3, p2 as vec3, p3 as vec3) as bool
#
#	in handle, anchor point, out handle が一直線に並んでいるかを返す
#	p1 :		in handle 座標
#	p2 :		anchor point 座標
#	p3 :		out handle 座標	
	
def get_smooth_handle(p1, p2, p3) :
	v1 = p2 - p1				#  ここに入る前に v1, v2 に長さがあることは確認済み
	v2 = p3 - p2
	v1.norm()
	v2.norm()
	return v1.dot(v2) >= 0.99999994
	
	
	
	
#  get_poisition_list(ob as shape, nop as int, closed as bool, mx2 as matrix) as vec3 list
#		ob :		掃引体交差方向線形状
#		nop :	ob の ポイント数
#		closed :	ob の開閉情報
#		mx2 :	ob の xshade.scene().local_to_world_matrix
#
#		ob の 隣接する3つの anchor point 座標の lisit の list [[vec3, vec3, vec3], [vec3, vec3, vec3}, ‥]  を返す

def get_poisition_list(ob, nop, closed, mx2) :
	pL = []

	p2 = vec3(ob.control_point(1).position)*mx2
	p1 =  vec3(ob.control_point(0).position)*mx2
	if not closed :	
		p0 = p1	
	else :														
		p0 = vec3(ob.control_point(nop - 1).position)*mx2
	pL.append([p0, p1, p2])
	
	for i in range(2, nop) :
		p0 = p1
		p1 = p2
		p2 = vec3(ob.control_point(i).position)*mx2
		pL.append([p0, p1, p2])
		
	p0 = p1
	p1 = p2
	if not closed :
		p2 = p1
	else :
		p2 = pL[0][1]
	pL.append([p0, p1, p2])
	
	return pL
	
	
	
	
#  s_flatness_list(bZs as vec3 list, bPs as vec3 matrix list, closed as bool) as bool list
#		bZs :		掃引中心線の control point 座標 list
#		bPs :		掃引中心線を含む bezier patch list
#		closed :	中心線の開閉情報
#
#		bPs の各 patch 区間の 非ねじれ list を返す
#
#		list の始端と終端には次なる項目を追加
#			中心線が閉じている場合、その隣接項目
#			中心線が開いている場合、False

def s_flatness_list(bZs, bPs, closed):
	flatL = []
	n = len(bPs)
	for i in range(n) :
		bz = bZs[i]
		bp = bPs[i]
		
		u1 = labo.bezier_line_tangent(bz, 0)	#  outhandle 側接線ベクトル
		u2 = labo.bezier_line_tangent(bz, 0.5)	#  中点での接線ベクトル
		
		if u1.abs2() == 0 :						#  重点
			flatL.append(False)
		
		else :
			v1 = labo.bezier_surface_v_tangent_2(bp, 0., 0.)	#  ここの関数が呼ばれる時点で、v1, v2  が求まることは既に担保されている
			v2 = labo.bezier_surface_v_tangent_2(bp, 0., 0.5)
			vn1 = u1*v1
			vn1.norm()
			vn2 = u2*v2
			vn2.norm()
			Q = quaternion().slerp(u1, u2)
			w = Q*vn1
			flatL.append(w.dot(vn2) >= 0.99985)	#   本来なら しきい値は 0.9999 とするが、ベースとなる自由曲面の持つ理想形状との誤差を考慮し、若干緩めに設定
		
	#   関数 h_sweep_handle( ) に与える引数の選択において、平坦性 list fL の隣接二値を比較するため、 list の先頭と終端に値を追加
	if not closed :
		flatL.insert(0, True)
		flatL.append(True)
	else :
		flatL.insert(0, flatL[len(flatL) - 1])
		flatL.append(flatL[1])
		
	return flatL
	
	

#  get_vt(w as vec3 list, q as vec3, r as float) as vec3
#
#		交差線形状の handle の方向を与える vector を返す
#		w[0] :	交差線形状の当該 point の一つ手前の point 座標
#		w[1] :	交差線形状の当該 point の座標
#		w[2] :	交差線形状の当該 point の一つ先の point 座標
#		q :		 中心線 point 座標
#		r :		中心線 boinding box size / 16

def get_vt(w, q, r) :
	#  w[0] = w[1] = w[2] となる三重点のケースではここに入らない
	v1 = w[0] - w[1]
	v2 = w[2] - w[1]
	v1.norm()
	v2.norm()                    
	v0 = v2 - v1
	v0.norm()
	
	v = w[1] - q
	a = v.norm()
	if a < r*1e-5 :		# w[1]  と q が一致しているとみなす
		vt = v0
	else :
		d = v.dot(v0)
		vt = v0 - d*v
		vt.norm()
		
	return vt
		
		
		
		
#  s_sweep(ob1 as shape, ob2 as shape, threshold1 as float, threshold2 as float, limitter1 as float = 10., limitter2 as float = 3.) as shape
#
#	掃引体を Shade 上に作成し、その shape を返す
#		ob1 :	掃引中心線形状
#		ob2 :	掃引断面形状
#		threshold1:	分割しきい値 ( 全曲率ベース )							##  2つの threshold  #############
#		threshold2:	分割しきい値 ( 面法線ベース )							##  2つの threshold  #############
#		limitter1 :	直線コーナーの scale liitter	
#		limitter2 :	曲線コーナーの scale liitter	

def s_sweep(ob1, ob2, threshold1, threshold2, limitter1 = 10., limitter2 = 3.) :	
	if (ob1 == None) or (ob2 == None) :
		return	
		
	##  引数として渡された object の check  ##################################################
	if (not isinstance(ob1, xshade.line)) or (not ob1.has_dad) :
		return
	if (not isinstance(ob1.dad, xshade.part)) or (ob1.dad.part_type != 1) :	
		return
	if (ob1.dad.number_of_sons < 2) or (ob1.dad.total_number_of_control_points < 4) :
		return
	if not isinstance(ob2, xshade.line) :														
		return
	##############################################################################
		
	#  check
	if (threshold1 < 0) or (threshold2 < 0) :								##  2つの threshold  #############
		return
	if limitter1 < 1 or limitter2 < 1 :
		return
	

	#  掃引中心線データ取得 nop, closed, bZ
	ob1.activate()
	nop = ob1.number_of_control_points	#  ポイント数
#	if nop < 2 :															##   check 済みにつき、不要  #############
#		return																##   check 済みにつき、不要  #############
	closed = ob1.closed						#  開閉情報
	bZs = labo.get_bezier(xshade, -1)		#  全区間の bezier data [matrix, matrix, … ]
	
	
	
	###   追加  #######################################################################
	
	#  掃引中心線データ取得  bP
	if not closed :
		n = nop - 1							#  U ブロック数
	else :
		n = nop
		
	ob0 = ob1.dad
	m = ob1.ordinal - ob0.ordinal - 1		#  V ブロック No.
	lastL = False

	if (not ob0.surface_closed) and (m == ob0.number_of_sons - 1) :
		m += -1
		lastL = True
		Mrv = matrix()
		Mrv.append([0, 0, 0, 1])
		Mrv.append([0, 0, 1, 0])
		Mrv.append([0, 1, 0, 0])
		Mrv.append([1, 0, 0, 0])
		
	ob0.activate()
	bPs = []					#  ob1 の bezier patch list
	for i in range(n) :
		bp = labo.get_bezier_patch(xshade, m, i)
		if lastL :
			bp = Mrv*bp			#  中心線が交差方向が開いた自由曲面の最後のラインの場合、V 方向の向きを逆転
		bPs.append(bp)
		
	#################################################################################
	
	
	
	# mx1, mx2, mx22, ps, pe
	mx1 = matrix(ob2.world_to_local_matrix)
	mx2 = matrix(ob2.local_to_world_matrix)
	if not closed :
		mx22 = matrix(ob1.local_to_world_matrix)
		ps = vec3(ob1.control_point(0).in_handle)*mx22			#  掃引体の交差方向 handle 作成時に使用
		pe = vec3(ob1.control_point(nop - 1).out_handle)*mx22	#  掃引体の交差方向 handle 作成時に使用
		
	#  r							
	bb = labo.vec3_bounding_box_size(bZs)	
	r = bb[3]/16												#  sweep_handle( ), get_vt( ) の引数に使用	
	
		
	#  掃引中心線 の handle 情報取得
	has_in_handle = []			#  掃引中心線の inhandle の有無のリスト
	has_out_handle = []			#  掃引中心線の outhandle の有無のリスト
	linked = []					#  掃引中心線の hanndle link リスト
	smooth_handle = []			#  中心線の smooth handle lリスト
	for i in range(nop) :
		has_in_handle.append(ob1.control_point(i).has_in_handle)
		has_out_handle.append(ob1.control_point(i).has_out_handle)
		linked.append(ob1.control_point(i).linked)
		if has_in_handle[i] and has_out_handle[i] :
			if i == 0 :
				if closed :
					smooth_handle.append(get_smooth_handle(bZs[nop - 1][2], bZs[0][0], bZs[0][1]))
				else :
					smooth_handle.append(False)	
			elif i == nop - 1 :
				if closed :
					smooth_handle.append(get_smooth_handle(bZs[nop - 2][2], bZs[nop - 1][0], bZs[nop - 1][1]))
				else :
					smooth_handle.append(False)		
			else :
				smooth_handle.append(get_smooth_handle(bZs[i - 1][2], bZs[i][0], bZs[i][1]))		
		else :
			smooth_handle.append(False)
			
		
	#  bZs, bPs を分割	
	#   屈曲点で分割
	bZs1 = []
	bPs1 = []																##  追加  #############							
	insert_list = []
	k = 0
	
	for bz in bZs :
		bp = bPs[k]															##  追加  #############
		k += 1
		t = bezier_inflection_point(bz)
		n = len(t)
		if n == 0 :				#  屈曲点なし
			bZs1.append(bz)
			bPs1.append(bp)													##  追加  #############
		else :
			t.insert(0, 0)
			t.append(1)
			n += 2
			for i in range(n - 1) :
				bZs1.append(labo.bezier_line_subdivision(bz, t[i], t[i + 1]))
				bPs1.append(labo.bezier_patch_subdivision(bp, 0, 1, t[i], t[i + 1]))		##  追加  #############
				insert_list.append(k)
				nop += 1
			insert_list.pop()
			nop += -1

	if len(insert_list) > 0 :
		insert_list.reverse()
		for k in insert_list :
			has_in_handle.insert(k, True)
			has_out_handle.insert(k, True)
			linked.insert(k, True)
			smooth_handle.insert(k, True)

	
	if threshold1 > 0 :	
		# threshold1 / threshold2 に従って 全曲率 あるいは 面法線基準で分割
		bZs2 = []
		bPs2 = []																			##  追加  #############
		insert_list = []
		k = 0
		
		for bz in bZs1 :
			bp = bPs1[k]																	##  追加  #############
			k += 1
			t = bezier_divide_by_total_curvature_or_torsion(bz, bp, threshold1, threshold2)	##  関数変更  #############
			if t == None :	#  面法線が求まらない場合、None が返され、途中終了する					##  追加  #############
				print '基準となる自由曲面で面法線が求まらないので、終了します。'						##  追加  #############
				return																		##  追加  #############
			n = len(t)
			if n == 0 : 
				bZs2.append(bz)
				bPs2.append(bp)																##  追加  #############
			else :
				for i in range(n - 1) :
					bZs2.append(labo.bezier_line_subdivision(bz, t[i], t[i + 1]))
					bPs2.append(labo.bezier_patch_subdivision(bp, 0, 1, t[i], t[i + 1]))	##  追加  #############
					insert_list.append(k)
					nop += 1
				insert_list.pop()
				nop += -1
					
		if len(insert_list) > 0 :
			insert_list.reverse()
			for k in insert_list :
				has_in_handle.insert(k, True)
				has_out_handle.insert(k, True)
				linked.insert(k, True)
				smooth_handle.insert(k, True)
				
		bZs = bZs2
		bPs = bPs2																			##  追加  #############						
											
	
	#  ob0 : 掃引体を格納する自由曲面
	ob2.activate()
	scene.create_surface_part(None)
	ob0 = scene.active_shape()				#   掃引体を格納する自由曲面
	ob0.surface_closed = closed
	ob2.place_child(1)

	
	#  掃引断面を配置
	p1 = bZs[0][0]
	Mt1 = matrix().translate(-p1[0], -p1[1], -p1[2])
	Q1 = quaternion()
	n = nop - 1
	v2_prev = None
	vn_prev = None																			##  追加  #############
		
	for i in range(n) :
		p2 = bZs[i][3]
		v2 = labo.bezier_line_tangent(bZs[i], 1)
		Mt2 = matrix().translate(p2[0], p2[1], p2[2])
		
		##  s_sweep_section_quaternion( ) に変更  ######################################################
		Q2 = s_sweep_section_quaternion(bZs[i], bPs[i])	
		
		#  Q2 が求まることは、既に実行した関数 bezier_divide_by_total_curvature_or_torsion( ) でチェック済み
#		if Q2 == None :									#  面法線が求まらない場合、途中終了
#			print '基準となる自由曲面で面法線が求まらないので、終了します。'
#			return
		######################################################################################
			
		Q2 = Q2*Q1
		Mr = Q2.matrix()
		
		if i == n - 1 and not closed :
			M = mx2*Mt1*Mr*Mt2*mx1
			
		else :
			v3 = labo.bezier_line_tangent(bZs[i + 1], 0)
			if v3.abs2() == 0 :							#  重点
				M = mx2*Mt1*Mr*Mt2*mx1
				if v2_prev == None :					#  最初の重点
					v2_prev = v2
					vn_prev = labo.bezier_surface_normal_2(bPs[i], 0., 1.)						##  追加  #############
				
			elif v2.abs2() == 0 :						#  最後の重点
				if v2_prev == None or v2_prev.abs2() == 0 :				#   始点が重点から始まっている
					M = mx2*Mt1*Mr*Mt2*mx1	
				else :
					vn = labo.bezier_surface_normal_2(bPs[i + 1], 0., 0.)								##  追加  #############
					Qc = labo.attitude_control_quaternion(v2_prev, vn_prev, v3, vn, 0)				##  変更  #############								
					Q2 = Qc*Q2							
					Mc = Qc.matrix()						
					M = mx2*Mt1*Mr*Mc*Mt2*mx1			
					
				v2_prev = None
				vn_prev = None																	##  追加  #############
				
			else :
				if v2.dot(v3) >= 0.99999994 :			#  掃引中心線の point が corner point  ではない
					M = mx2*Mt1*Mr*Mt2*mx1
				
				elif v2.dot(v3) <= -0.99999994 :		#  掃引中心線の point が 180 度 corner point
					M = mx2*Mt1*Mr*Mt2*mx1
					vn = labo.bezier_surface_normal_2(bPs[i + 1], 0., 0.)							##  追加  #############
					Qc = labo.attitude_control_quaternion(v2, -vn, v3, vn, 0)						##  変更  #############
					Q2 = Qc*Q2							#   次の point での掃引断面配置に必要

				else :									#  掃引中心線の point が corner point 
					w = v2 + v3
					w.norm()
					v1 = labo.bezier_line_tangent(bZs[i], 0)
					v4 = labo.bezier_line_tangent(bZs[i + 1], 1)

					if v1.dot(v2) + v3.dot(v4) > 1.98 :	#  直線コーナー
						limitter = limitter1
					else :								#  曲線コーナー
						limitter = limitter2
					Q3 = quaternion().slerp(v2, w)
					Mq = Q3.matrix()
					
					v = v2 - v3
					cosT = w.dot(v2)
					Q4 = quaternion().slerp(v, vec3(0, 1, 0))
					Mqs1 = Q4.matrix()
					Q5 = Q4.inverse()
					Mqs2 = Q5.matrix()
					Ms = matrix().scale(1, min(limitter, 1/cosT), 1)	
					M = mx2*Mt1*Mr*Mq*Mqs1*Ms*Mqs2*Mt2*mx1

					Q6 = quaternion().slerp(v2, v3)										##  quaternion().slerp( ) に変更  #############
					Q2 = Q6*Q2						#   次の point での掃引断面配置に必要
					
		ob2.copy_object(M)
		ob = ob2.bro
		ob.place_brother(i)
	
		Q1 = Q2

	
	#  掃引体の交差方向 handle をセット
	n = ob2.number_of_control_points
	ob0.switch()
	ob = ob0.son
	
	flatL = s_flatness_list(bZs, bPs, closed)					#  中心線の各 bezier patch 区間の平坦性 list				##  変更  #############
	
	for i in range(n) :
		ob = ob.bro												#  掃引体内部の交差方向線形状
		pL = get_poisition_list(ob, nop, closed, mx2)			# ob の 隣接する3つの anchor point 座標の lisit
		vt = None												#  交差方向 handle の方向 ( inhandle と outhandle が平行な場合 )	
		
		for j in range(nop) :
			w = pL[j][1]										#  交差方向線形状の point 座標
			
			if not has_in_handle[j] :
				inH = w											#  inhandle 座標
			else :
				if (j == 0) and not closed :
					inH = w + ps - bZs[0][0]					#  inhandle 座標
				else :
					if j == 0 :
						k = nop - 1
					else :
						k = j - 1
					vp = sweep_plane_normal(bZs[k][3], bZs[k][2], bZs[k][1], bZs[k][0])
					u = bZs[k][2] - bZs[k][3] 
					u.norm()
					p = bZs[k][2]
				
					if flatL[j] or flatL[j + 1] or not smooth_handle[j] :	#  掃引中心線の当該 point  前後の bezier 区間のいずれかが非ねじれ面、あるいは not smooth handle	
						inH = h_sweep_handle(w, u, r, p, vp)	#  inhandle 座標
					else :
						vt = get_vt(pL[j], bZs[k][3], r)
						inH = h_sweep_handle(w, -vt, r, p, vp)	#  inhandle 座標
					
	
			if not has_out_handle[j] :
				outH = w										#  outhandle 座標
			else :
				if (j == nop - 1) and not closed :
					outH = w + pe - bZs[nop - 2][3]				#  outhandle 座標
				else :
					vp = sweep_plane_normal(bZs[j][0], bZs[j][1], bZs[j][2], bZs[j][3])
					u = bZs[j][1] - bZs[j][0] 
					u.norm()
					p = bZs[j][1]
	
					if flatL[j] or flatL[j + 1] or not smooth_handle[j] :	#  掃引中心線の当該 point  前後の bezier 区間のいずれかが非ねじれ面、あるいは not smooth handle
						outH = h_sweep_handle(w, u, r, p, vp)	#  outhandle 座標
					else :
						if vt == None :
							vt = get_vt(pL[j], bZs[j][0], r)
						outH = h_sweep_handle(w, vt, r, p, vp)	#  outhandle 座標
						vt = None																			
						

			ob.control_point(j).in_handle = inH*mx1
			ob.control_point(j).out_handle = outH*mx1					
			ob.control_point(j).linked = linked[j]
	
		
	ob0.switch()
	ob0.activate()
	ob0.disclosed = False
	
	return ob0
	

	
	
#  Shade  上で 掃引中心線 , 掃引断面 の順で選択して実行

scene = xshade.scene()

#  dialog uuid は Online GUID Generator により生成		https://www.guidgenerator.com
dialog = xshade.create_dialog_with_uuid('91b274c6-179a-4a48-a390-cfa6f74a1701')

dialog.begin_tab_group('分割')
level_ = dialog.append_selection('level/90度/60度/45度/30度/22.5度/15度')
dialog.end_tab_group()

dialog.begin_tab_group('limitter')
limitter1_ = dialog.append_float('直線コーナー', '( defalt 10 )')
limitter2_ = dialog.append_float('曲線コーナー', '( defalt 3 )')
dialog.end_tab_group()

dialog.set_value(4, 10)
dialog.set_value(5, 3)

if dialog.ask('曲面に沿った掃引'):
	level = dialog.get_value(level_)
	from math import pi
	if level == 0 :	
		threshold1 = pi/2	
	elif level == 1 :		
		threshold1 = pi/3
	elif level == 2 :	
		threshold1 = pi/4
	elif level == 3 :	
		threshold1 = pi/6
	elif level == 4 :	
		threshold1 = pi/8
	elif level == 5 :	
		threshold1 = pi/12
	threshold2 = threshold1
	limitter1 = dialog.get_value(limitter1_)
	limitter2 = dialog.get_value(limitter2_)

	scene.exit_modify_mode()
	[ob1, ob2] = scene.active_shapes
	ob2.copy_object(None)
	ob3 = ob2.bro
	ob4 = s_sweep(ob1, ob3, threshold1, threshold2, limitter1, limitter2)
	if ob4 != None :
		ob4.activate()    
	else :
		scene.active_shapes = [ob1, ob2]

#8

30 - 8  サンプル形状


サンプル形状出力 script として、[ script 30 - 5 ], [ script 30 - 6 ] を記します。
      [ fig 30 - 10 ]
[ script 30 - 5 ]
from vec3 import *
from math import sin, cos, pi


scene = xshade.scene()
scene.exit_modify_mode()
scene.lathe_rotation = [1, 0, 0, 0]


#   面
scale = 5000
ds = 2*pi/24
dt = 2*pi/24

scene.create_surface_part('基準面')
ob1 = scene.active_shape()
ob1.surface_closed = True
#ob1.disclosed = False

scene.begin_creating()
for j in range(24) :
	p = []
	for i in range(24) :
		s = i*ds
		t = j*dt
		x = cos(s + t)/(2**0.5 + cos(t - s))
		z = -sin(s + t)/(2**0.5 + cos(t - s))
		y = sin(t - s)/(2**0.5 + cos(t - s))
		p.append(scale*vec3(x, y, z))
	
	scene.create_line(None, p, True)
scene.end_creating()

ob1.activate()
scene.smooth()
ob1.switch()
scene.smooth()
ob1.switch()

ob2 = ob1.son.bro
ob2.name = '中心線'
	
p0 = ob2.control_point(0).position


#  断面
r = 0.5
p = []
p.append(scale*vec3(-r, 0, 0))
p.append(scale*vec3(0, r/2, 0))
p.append(scale*vec3(0, -r/2, 0))

scene.begin_creating()
scene.create_line('断面', p, True)
scene.end_creating()

scene.move_object(None, None, [-45, 0, 0], p0)
ob3 = scene.active_shape()


scene.active_shapes = [ob2, ob3]

[ script 30 - 6 ]
from vec3 import *
from matrix import *
import labo


scene = xshade.scene()
scene.exit_modify_mode()
scene.lathe_rotation = [1, 0, 0, 0]


#   面
scale = 10000


scene.create_surface_part('基準面')
ob1 = scene.active_shape()
ob1.surface_closed = True
#ob1.disclosed = False

scene.begin_creating()
scene.begin_line(None, False)
scene.append_point(scale*vec3(-1, 1, 0.5), None, scale*vec3(-0.7, -0.3, 0.5), None, None)
scene.append_point(scale*vec3(1.5, -1.5, 0.5), scale*vec3(0, -1.5, 0.5), None, None, None)
scene.end_line()
scene.end_creating()

ob = scene.active_shape()
s = 1.
for t in [0.8, 0.6, 0.5, 0.4, 0.3, 0.15] :
	ob.insert_control_point(t/s)
	s = t
	
scene.move_object([0, 0, 0], None, [0, 30, 0], None)
for i in range(17) :
	scene.copy_object([0, 0, 0], None, [0, 20, 0], None)
ob1.activate()
ob1.switch()
scene.smooth()
ob1 = scene.active_shape()
ob1.switch()
ob2 = ob1.son.bro
ob2.name = '中心線'


#  断面
bp = labo.get_bezier_patch(xshade, 0, 0)
u = labo.bezier_surface_u_tangent(bp, 0, 0)
v = labo.bezier_surface_normal(bp, 0, 0)
pp = labo.bezier_surface_position(bp, 0, 0)

r = 0.5
p = []
p.append(scale*vec3(0, r, 0))
p.append(scale*vec3(-r/2,0, 0))
p.append(scale*vec3(r/2, 0, 0))

scene.begin_creating()
scene.create_line('断面', p, True)
scene.end_creating()
ob3 = scene.active_shape()
Q = labo.attitude_control_quaternion(vec3(0, 0, 1), vec3(0, 1, 0), u, v, 1)
Mq = Q.matrix()
Mt = matrix().translate(pp[0], pp[1], pp[2])
ob3.move_object(Mq*Mt)


scene.active_shapes = [ob2, ob3]