Shade3D 公式

16 Bezier Patch 分割 [ Shade Labo ]

16 - 1 区分曲面


関連記事:
  • 15 Bezier Line 分割

bezier patch [ **P** ] で構成される bezier surface **P**( s, t ) の区分曲面 **P”**( s, t )( s1 <= s <= s2, t1 <= t <= t2 )の bezier patch [ **P”** ] を求めます。
まず最初に区間 s1 <= s <= s2 に対する区分曲面 **P’**( s, t )を考えます。

bezier line の区分曲線を求める [ 式 15 - 10 ] がそのまま bezier patch にも当てはまります。


          [ 式 15 - 10 ]


次に区間 t1 <= t <= t2 での区分曲面を考えると、次の手続きにより同様に [ 式 15 - 10 ] を適用できます。

  • [ P ] の転置 matrix をとる
  • それに対して parameter t を用いて [ 式 15 - 10 ] を適用
  • 得られた matrix の転置 matrix をとる

よって、parameter s1 <= s <= s2, t1 <= t <= t2 における区分曲面の bezier patch [ **P”** ] は次にように得られます。
           [ 式 16 - 1 ]
以上から、区分曲面を求める関数は次のようになり、labo に登録してあります。
[ script 16 - 1 ]
#  bezier_patch_subdivision( bP as vec3 matrix, s1 as float, s2 as float, t1 as float, t2 as float ) as vec3 matrix
#		bezier patch bP の parameter 区間 s1 ~ s2, t1 ~ t2 における区分曲面の  bezier patch を返す

def bezier_patch_subdivision(bP, s1, s2, t1, t2) :
 	if bP == None :
 	 	return None

	bP1 = bezier_line_subdivision(bP, s1, s2)
	bP1.transpose()
	bP2 = bezier_line_subdivision(bP1, t1, t2)
	bP2.transpose()		
	return bP2

次の script では選択された自由曲面の 区間 [ 0, 0 ] 面の分割面をを対角線状に3つ並べます。
[ script 16 - 2 ]
import labo
	
scene = xshade.scene()

#  bezier patch 取得
bP = labo.get_bezier_patch(xshade, 0, 0)

#  区分面
bP1 = labo.bezier_patch_subdivision(bP, 0, 1./3, 0, 0.25)
bP2 = labo.bezier_patch_subdivision(bP, 1./3, 2./3, 0.25, 0.7)
bP3 = labo.bezier_patch_subdivision(bP, 2./3, 1, 0.7, 1)

#  区分面出力
scene.create_part('区分曲面')
labo.make_bezier_surface(xshade, bP1)
labo.make_bezier_surface(xshade, bP2)
labo.make_bezier_surface(xshade, bP3)
scene.select_parent(1)

16 - 2 分割後の Bezier Patch


関連記事:
  • 12 Bezier Patch

[ script 16 - 3 ] では、サンプル自由曲面を分割した自由曲面を作成し、次の2つの bezier patchを出力します。
  • script 内で 分割計算を行った時に得られる bezier patch
  • Shade 上に出力された分割面から得られる bezier patch

[ script 16 - 3 ]
import labo
from matrix import *
	
scene = xshade.scene()
scene.create_part()

#  sample  bezier surface
bP = []
bP.append([[0, 10000, 0], [0, 10000, 7000], [0, -3000, 10000], [0, -10000, 10000]])
bP.append([[0, 10000, 0], None, None, [5000, -4000, 10000]])
bP.append([[0, 10000, 0], None, None, [10000, -6000, 5000]])
bP.append([[0, 10000, 0], [5000, 10000, 0], [10000, 8000, 0], [10000, 2000, 0]])
bP = matrix(bP)
labo.make_bezier_surface(xshade, bP, '<<sample')

#  区分面
labo.set_bezier_patch(bP)
bP1 = labo.bezier_patch_subdivision(bP, 0, 0.8, 0, 0.8)
labo.make_bezier_surface(xshade, bP1, '<<区分面')
bP2 = labo.get_bezier_patch(xshade, 0, 0)		#  Shade に出力された自由曲面から patch を取得

#  bezier patch を出力
labo.make_bezier_patch(xshade, bP1, 'script 内で計算')
labo.make_bezier_patch(xshade, bP2, 'Shade から読み取り')

scene.select_parent(1)

結果は以下のようになり、2つの bezier patch は異なってしまいます。
                    [ fig 16 - 1 ]
16 個の control point 座標の内、周囲の12個は当然一致しますが、中央4個の座標が異なります。

[ script 16 - 3 ] では Shade 上に出力した自由曲面から bzier patch を読み込むという遠回しな方法を用いていますが、要は、
  • [ 式 16 - 1 ] で分割した bezier patch はオリジナル面を正しくトレースする

  • Shade の自由曲面では与えられた bezier patch の中央4座標に対して独自のルールに基づいた再構成を行う( 12 Bezier Patch 参照 )

  • よって、分割の境界線は元の面をトレースするが、内部の形状は同一とは限らない

  • 形状にある程度の対称性が保たれていれば同一面を維持できることもある

ということです。


分割におけるこの問題は、自由曲面への線形状の追加の場合でも同様に発生します。( 線形状の追加とは分割された面をつなげた状態にしたものです )
この問題はさらに深い問題へとつながります。
  • 分割された面のデータを Shade 上から取得してをさらに分割すると、もはや分割境界線ですらオリジナル面をトレースできない

  • Shade 上で複数回の分割を行う場合、同じ分割でも 分割の順序によって形状の形が変わる


これに対し、script 内でオリジナルの bezier patch を保持し( 中央4個のデータを [ 式 12 - 1 ] での補正をかけない状態のまま保持 )、script 内部で分割し出力すれば、分割を繰り返したとしても、分割の順序にも関係なく、**少なくとも分割境界線はオリジナル面をトレースし続けることが可能です**。

16 - 3 Bezier Line 分割での lateral handle


関連記事:
  • 13 Bezier Surface 基礎

bezier line を分割する場合の交差ハンドルの扱いは、bezier patch の分割と同じロジックで行われます。

bezier line を parameter t で分割したときの交差方向の handle 座標は、bezier line を bezier patch の第一列に見立てて、parameter t における曲面上の V 方向線形状 Pv の control point 配列を求める場合の方法がそのまま適用されます。


          [ fig 16 - 2 ]

          [ 式 16 - 2 ]


つまり、bezier line を parameter t で分割したときの交差方向の handle 座標 Pv1 は control point 配列 [ P10, P11, P12, P13 ] なる bezier line の parameter t における座標になります。


[ script 16 - 5 ] はその検証 script です。

検証結果

          [ fig 16 - 3 ]


[ script 16 - 5 ]

from vec3 import *

scene = xshade.scene()
scene.create_part('test')
ob0 = scene.active_shape()
	
#  サンプル線形状
u = vec3(0, 1, 0)
v = vec3(1, 0, 0)
p1 = 10000*u
p2 = -p1
h1 = 2500*v
h2 = 10000*v

scene.begin_creating()
scene.begin_line(None, False)
scene.append_point(p1, None, p1/2, None, p1 + h1)
scene.append_point(p2, p2/2, None, None, p2 + h2)
scene.end_line()
scene.end_creating()
ob1 = scene.active_shape()
ob1.copy_object(None)
ob2 = ob1.bro
ob2.insert_control_point(0.5)
ob2.activate()
scene.create_surface_part()
ob3 = scene.active_shape()
ob2.place_child(1)
ob3.switch()

ob3.activate()
scene.begin_creating()
scene.begin_line(None, False)
scene.append_point(p1 + h1, None, p1/2 + h1, None, None)
scene.append_point(p2 + h2, p2/2 + h2, None, None, None)
scene.end_line()
scene.end_creating()

ob4 = scene.active_shape()
ob4.copy_object(None)
ob5 = ob4.bro
ob5.insert_control_point(0.5)

ob0.activate()
scene.enter_modify_mode()
scene.select_all()


**16.2** で触れたように、bezier patch の分割では分割の順序によって形状が変わってしまいますが、bezier line 分割時の交差方向 handle についても同じ現象が生じ、分割順の違いによってハンドル長さが異なってしまいます。
[ script 16 - 6 ] はその検証 script で、次の操作を行います。
  • 交差ハンドルを持つ線形状を作成し、その copy を作る
  • 2つの線形状に対し、1/3 と 2/3 の2点にアンカーポイントを追加
  • 2つの線形状ではポイント追加の順序を変えておく
  • 追加されたポイントにおける交差方向のハンドルを比較する

**検証結果**

          [ fig 16 - 4 ]


[ script 16 - 6 ]

from vec3 import *


scene = xshade.scene()
scene.create_part('test')
ob0 = scene.active_shape()
	
#  サンプル線形状
u = vec3(0, 1, 0)
v = vec3(1, 0, 0)
p1 = 10000*u
p2 = -p1
h1 = 2500*v
h2 = 10000*v

scene.begin_creating()
scene.begin_line(None, False)
scene.append_point(p1, None, p1/2, None, p1 + h1)
scene.append_point(p2, p2/2, None, None, p2 + h2)
scene.end_line()
scene.end_creating()
ob1 = scene.active_shape()
ob1.copy_object(None)
ob2 = ob1.bro

#  線形状を分割
ob1.insert_control_point(1./3)
ob1.insert_control_point(1.5)
ob2.insert_control_point(2./3)
ob2.insert_control_point(0.5)

scene.create_surface_part()
ob3 = scene.active_shape()
ob1.place_child(1)
ob2.activate()
scene.create_surface_part()
ob4 = scene.active_shape()
ob2.place_child(1)

ob3.switch()
ob4.switch()
ob0.activate()
scene.enter_modify_mode()
scene.select_all()

コーナーでの自由形式の表面テッセレーション動作 (Shade16, Shade 17.1 以降)

関連記事:

  • 12 - 7 Bezier Patch 内挿法( Shade 17 以降 )

Shade 17.1以降、4つの内部制御点の計算方法は、反対側のハンドルが特定の制御点の最終位置に寄与する差動システムで機能します。
From Shade 17.1, the way the four internal control points are calculated works in a differential system where the opposite handle contributes to the final position of a particular control point.

Shadeの変形システムは、対向するハンドルが自由曲面のテッセレーション動作に寄与するという点で、Shade 17.1の方法と同様に機能します。
The deformation system in Shade works similarly to the method of Shade 17.1, in that opposing handles contribute to the tessellation behavior of the free-form surface.

Shade 17.1の方法を使用した自由曲面の場合、コーナーでのテッセレーション動作は通常のベジェサーフェスと同じです。
For a free-form surface using the method of Shade 17.1, the tessellation behavior at the corner is the same as a normal Bezier surface:

shade17_corner

ただし、従来の方法を使用したフリーフォームサーフェスでは、コーナーを拡大すると、Shade 17.1メソッドのように見えるものが表示されます。
However, with a freeform surface using the legacy method, we see something that looks like the Shade 17.1 method when zooming in on a corner:

shade16_corner

よく調べてみると、ハンドル長が等しくない2D自由曲面のテッセレーション動作は次のとおりです。
Investigating closer, it turns out that the tessellation behavior for a 2D freeform surface with unequal handle lengths is the following:

一般的な場合、ハンドルの荷重みは、指定されたハンドルペア(U1、U2など)のハンドル間の角度、およびそれらの長さの比率によって変調されます。前者は3次多項式で制御され、後者は前述の多項式の垂直スケーリングを制御します。
多項式の形状により、Shadeはその領域を[0、asin(0.99))に制限します。
For the general case, the handle weighting is modulated by the angle between the handles in a given handle pair (e.g. U1, U2), and the ratio of their lengths. The former is controlled by a third degree polynomial, while the latter controls the vertical scaling of the aforementioned polynomial.
Due to the shape of the polynomial, Shade limits its domain to [0, asin(0.99)).

[image 4]

labo171.zip (13。8 キロバイト)
The following function is registered in the labo module Ver. 1.7.1:

def set_blending_weights(bP) :
if bP == None :
	return
	
bb = vec3_bounding_box_size(bP)			#  bP の bounding box size
epsilon = bb[3]/100000						#  0 handle 判定しきい値

#  overlap :	patch の一点収束 ( 極 ) 存在状態
pL = []
pL.append(bP[0])										#   第一行
pL.append([bP[0][0], bP[1][0], bP[2][0], bP[3][0]])		#   第一列
pL.append([bP[0][3], bP[1][3], bP[2][3], bP[3][3]])		#   第四列
pL.append(bP[3])										#   第四行

overlap = []
for p in pL :
	result = True
	for i in range(1,4) :
		v = p[i] - p[i - 1]
		if v.abs() > epsilon :
			result = False						
			break
	if result :
		overlap.append(True)				#  一点収束している
	else :
		overlap.append(False)				#  一点収束していない
		

#  4 つの bezier patch 荷重をセット
for [[i, j], [ii, jj], k1, k2] in [  [[0, 0],[1, 1], 0, 1]  ,  [[0, 3],[1, -1], 0, 2]  ,  [[3, 0], [-1, 1], 3, 1]  ,  [[3, 3], [-1, -1], 3, 2]  ] :
	
	p00 = bP[i][j]
	
	p01 = bP[i][j + jj]
	p10 = bP[i + ii][j]
	
	p30 = bP[i + 3*ii][j]
	p31 = bP[i + 3*ii][j + jj]
	
	p03 = bP[i][j + 3*jj]
	p13 = bP[i + ii][j + 3*jj]

	u1 = p01 - p00
	u2 = p31 - p30
	v1 = p10 - p00
	v2 = p13 - p03
	msin = math.asin(0.99)
	mcos = math.acos(0.99)
	lim = msin+(2.*mcos)
	u1p1 = vec3(u1)
	u2p1 = vec3(u2)
	unorm1 = u1p1.norm()
	unorm2 = u2p1.norm()
	uy = (u1p1 - u2p1)
	ux = (u1p1 + u2p1)
	uynorm = uy.norm()
	uxnorm = ux.norm()
	utheta = math.atan2((uynorm),(uxnorm))*2.
	v1p1 = vec3(v1)
	v2p1 = vec3(v2)
	vnorm1 = v1p1.norm()
	vnorm2 = v2p1.norm()
	vy = (v1p1 - v2p1)
	vx = (v1p1 + v2p1)
	vynorm = vy.norm()
	vxnorm = vx.norm()
	vtheta = math.atan2((vynorm),(vxnorm))*2.
	ustrength = (3.+((-0.116729182126221*(utheta**3)) - (0.0406410648893143*(utheta**2)) - (0.275506612356452*(utheta)))) # 経験的に導出された関数
	ustrength2 = (3.+((-0.116729182126221*((msin-(utheta-lim))**3)) - (0.0406410648893143*((msin-(utheta-lim))**2)) - (0.275506612356452*(msin-(utheta-lim)))))
	strength3 = (3.+((-0.116729182126221*(msin**3)) - (0.0406410648893143*(msin**2)) - (0.275506612356452*(msin))))
	vstrength = 3.+((-0.116729182126221*(vtheta**3)) - (0.0406410648893143*(vtheta**2)) - (0.275506612356452*(vtheta)))
	vstrength2 = (3.+((-0.116729182126221*((msin-(vtheta-lim))**3)) - (0.0406410648893143*((msin-(vtheta-lim))**2)) - (0.275506612356452*(msin-(vtheta-lim)))))

	
	#  u, v
	if overlap[k1] and overlap[k2] :					#  行, 列 両方向が一点収束 ( 極 )
		u = vec3(0, 0, 0)
		v = vec3(0, 0, 0)
	
	elif overlap[k1] :									#  行方向が一点収束 ( 極 )
		p20 = bP[i + 2*ii][j]
		v3 = p20 - p00
		r = v3.abs()
		if r > epsilon :
			u = u2*v1.abs()/r
		else :
			u = vec3(0, 0, 0)
		v = v1
		
	elif overlap[k2] :									#  列方向が一点収束 ( 極 )
		p02 = bP[i][j + 2*jj]
		u3 = p02 - p00
		r = u3.abs()
		if r > epsilon :
			v = v2*u1.abs()/r
		else :
			v = vec3(0, 0, 0)
		u = u1
		
	else :											#  行 / 列 方向が一点収束 ( 極 ) ではない
		if u1.abs() > epsilon :						#  行方向 handle あり
			if u2.abs() > epsilon :					#  反対側の行方向 handle あり
				if u1.length() > u2.length() > 0 :
					umod1 = ((3-((0.816456082525*(u2.length()/u1.length()))+2.183543917475312))/0.816456082525) # 経験的に導出された関数
					if utheta < msin :
						u = (((((3.-ustrength)*umod1)+3)*u1) + (u1p1*u2.length()))/4
					elif utheta > lim :
						u = (((((3.-ustrength2)*umod1)+3)*u1) + (u1p1*u2.length()))/4
					elif msin <= utheta <= lim :
						u = (((((3.-strength3)*umod1)+3)*u1) + (u1p1*u2.length()))/4
				else :
					if u2.length() > u1.length() > 0 :
						umod2 = ((3-((0.816456082525*(u1.length()/u2.length()))+2.183543917475312))/0.816456082525)
						if utheta < msin :
							u = (((((3.-ustrength)*umod2)+3)*u1) + (u1p1*u2.length()))/4
						elif utheta > lim :
							u = (((((3.-ustrength2)*umod2)+3)*u1) + (u1p1*u2.length()))/4
						elif msin <= utheta <= lim :
							u = (((((3.-strength3)*umod2)+3)*u1) + (u1p1*u2.length()))/4
			else :
				u = u1
		else :									#  反対側の行方向 handle なし
			u = vec3(0, 0, 0)
		

		if v1.abs() > epsilon :						#  列方向 handle あり
			if v2.abs() > epsilon :					#  反対側の列方向 handle あり
				if v1.length() > v2.length() > 0 :
					vmod1 = ((3-((0.816456082525*(v2.length()/v1.length()))+2.183543917475312))/0.816456082525)
					if vtheta < msin :
						v = (((((3.-vstrength)*vmod1)+3)*v1) + (v1p1*v2.length()))/4
					elif vtheta > lim :
						v = (((((3.-vstrength2)*vmod1)+3)*v1) + (v1p1*v2.length()))/4
					elif msin <= vtheta <= lim :
						v = (((((3.-strength3)*vmod1)+3)*v1) + (v1p1*v2.length()))/4
				else :
					if v2.length() > v1.length() > 0 :
						vmod2 = ((3-((0.816456082525*(v1.length()/v2.length()))+2.183543917475312))/0.816456082525)
						if vtheta < msin :
							v = (((((3.-vstrength)*vmod2)+3)*v1) + (v1p1*v2.length()))/4
						elif vtheta > lim :
							v = (((((3.-vstrength2)*vmod2)+3)*v1) + (v1p1*v2.length()))/4
						elif msin <= utheta <= lim :
							v = (((((3.-vtrength3)*vmod2)+3)*v1) + (v1p1*v2.length()))/4
			else :									#  反対側の列方向 handle なし
				v = v1
		else :
			v = vec3(0, 0, 0)
	
	bP[i + ii][j + jj] = p00 + u + v`整形済みテキスト`