Shade3D 公式

11 精度 [ Shade Labo ]


#1

11 - 1 有効桁数


有効桁は次のようになります。

浮動小数点 4 byte single の有効桁は、2進数で24桁、10進数では約7桁
浮動小数点 8 byte double の有効桁は、2進数で53桁、10進数では約15桁

python の float は 8byte double
Shade から script に data として渡されるのは 4byte single

従って python 内で計算をしている限りは double の有効桁は保たれますが、一旦形状として Shade に出力すると、そこから得られるデータの精度は single に落ちてしまいます。


#2

11 - 2 平行判定


2016 / 01 / 23      **[ script 11 - 3 ] にコード追加**
精度の一例として、平行判定を考えてみます。

ベクトル v1, v2 の外積を用いた平行判定は | v1 x v2 | = 0 ですが、実際には丸め誤差があるので、 | v1 x v2 | < ε としなければなりません。


この判定しきい値 ε は次の基準で求めるのが妥当でしょう。

  1. モデルのサイズに依存しない相対誤差ベースで判定を行う

  2. Shade で平行であるつもりでモデリングしたものに対し、正しい判定結果を返すこと



**1)相対誤差ベースでの判定**
今、判定対象となる二つのベクトル **v1**, **v2** の長さに極端な差がないとことを前提とします。( 一方が極端に短い場合、後述する **0** ベクトル判定という別の話が出てきます)

相対誤差をどのように定義するかは色々な考えがありますが、ここでは次のような判定を考えます。

          [ 式 11 - 1 ]


幾何的な解釈は、

  • v1, v2 で作られる扁平平行四辺形の面積 S1 が、一辺が | v1 |, | v2 | の長方形の面積 S0 の ε 倍よりも小さければ平行と見なす

あるいは、

  • v1, v2 のなす角をθとしたとき、sin(θ) <= εであれば平行と見なす


              [ fig 11 - 1 ]



2)平行であるつもりでモデリングしたものに対し、正しい判定結果を返す


これは少々厄介な問題を含みます。

Shade の X 軸上に直線状に並んだ 3点 p1, p2, p3 の平行判定を考えます。

ベクトル V は2点の座標値の差で与えられますから、座標値の有効桁数( Shade 上のデータは 7 桁 )によって保持できる精度が異なってしまいます。


p1 = ( 1.000000, 0, 0 )
p2 = ( 1.100000, 0, 0 )
V = p2 - p1 = ( 0.100000, 0, 0 )
| V | = 0.100000          V の有効桁は 6 桁


p1 = ( 100.0000, 0, 0 )
p2 = ( 100.1000, 0, 0 )
V = p2 - p1 = ( 0.1000, 0, 0 )
| V | = 0.1000          V の有効桁は 4 桁


一見すると両者とも同じ値に見えますが、その保持する有効桁数は異なります。

これは Shade 上でモデルのある部分の平行判定を行うとき、同一の判定しきい値 ε を用いても、モデルが原点付近にあるか離れた所にあるかによって、異なる結果を得てしまう可能性を示唆します。

これは点列 p1, p2, p3 に回転やスケーリングを施すと顕わになります。


次の script では同じモデルを原点付近と原点から離れた所に置いて、回転やスケーリングを施した後、座標値を読み取って| v1 x v2 | / | v1 | | v2 | なる比を求めます。


[ script 11 - 1 ]

import labo
from vec3 import *


def test(scale) :
	p1 = vec3(0, 0, 0)
	p2 = vec3(1, 0, 0)
	p3 = vec3(2, 0, 0)
	pL = [p1, p2, p3]				#  一直線上に並ぶ3点

	#  Shade 上に出力し、移動, 回転 , スケール処理	
	labo.make_simple_line(xshade, pL, False)
	scene.move_object([0, 0, 0], None, None, [0.8*scale, 1.4*scale, -1*scale])
	scene.move_object([0, 0, 0], None, [0, -115, 0], None)
	scene.move_object([0, 0, 0], None, [40, 0, 0], None)
	scene.move_object([0, 0, 0], None, [0, 0, 48], None)
	scene.move_object([0, 0, 0], [1.5, 0.7, 0.8], None, None)

	ob = scene.active_shape()
	p1 = vec3(ob.control_point(0).position)
	p2 = vec3(ob.control_point(1).position)
	p3 =vec3(ob.control_point(2).position)

	v1 = p1 - p2
	v2 = p3 - p2
	v = v2 - v1
	ss = v1.abs2()*v2.abs2()		
	v = v1*v2
	s = v.abs2()

	return (s/ss)**0.5


scene = xshade.scene()

print '原点付近 : ', test(1)
print 'scale 100 離れた位置 : ', test(100)
print 'scale 1000 離れた位置 : ', test(1000)

結果は,

     原点付近 :                         4.034192451e-07
     scale 100 離れた位置 :     1.59501635829e-05
     scale 1000 離れた位置 :   0.000259411825431


このように平行判定するベクトルのサイズと、そのベクトルが定義される元座標の位置によって、必要な ε のサイズが決まることになりますが、ケースによって都度 ε を決めるのも煩わしく、予め代表モデルケースを想定して定数としておくのが妥当でしょう。

たとえば、ε = 5e-4 とすれば、原点から 10m 離れた位置で、1cm サイズのベクトルに対して、丸め誤差があっても正しい判定が期待できます。




以上から一つの平行判定の script として、次が考えられます。


[ script 11 - 2 ]

#  is_parallel(v1 as vec3, v2 as vec3,  epsilon as double = 5e-4) as bool
#
#	2 つのベクトル v1, v2 が その向きを問わず平行であるか否かを返す

def is_parallel(v1, v2, epsilon = 5e-4) :
	ss =  v1.abs2()*v2.abs2()
	v = v1*v2
	s = v.abs2()
	return s/ss <= epsilon**2		#  面積の二乗で判定


さらに、ベクトル **v1**, **v2** の長さ L1, L2 の比にも注意を払います。

いま L1, L2 が同じ L と仮定し、2 つのベクトルで作られる扁平平行四辺形の面積を S とすると、平行判定のしきい面積は、

          [ 式 11 - 2 ]

となりますから、L1 << L2 の場合、

          [ 式 11 - 3 ]

が平行判定の限界条件、つまり、短い方のベクトルは長い方のベクトルに対して実質長さのない 0 ベクトルとみなされる 条件とすることができます。


従って 0 ベクトル判定を含めた、より一般的な平行判定は次のようになり、、module labo に登録してあります。


[ script 11 - 3 ]

#  is_parallel(v1 as vec3, v2 as vec3,  epsilon as double = 5e-4) as bool
#
#	2 つのベクトル v1, v2 が その向きを問わず平行であるか否かを返す
#

def is_parallel(v1, v2, epsilon = 5e-4) :
	L1 = v1.abs2()
	L2 = v2.abs2()
	ratio = min(L1, L2)/max(L1, L2)
	
	if ratio <= epsilon :
		return False					#   一方が 0 ベクトルと見なされる
	else :
		v = v1*v2
		s = v.abs2()
		return s/(L1*L2) <= epsilon**2	#  面積の二乗で判定


#  2016 / 01 / 23	追加
#
#	opposite = False ならば、
#		2 つのベクトル v1, v2 が 同じ方向で平行であるか否かを返す
#	opposite = True ならば、
#		2 つのベクトル v1, v2 が逆向きでで平行であるか否かを返す
 
def is_parallel_2(v1, v2, opposite = False, epsilon = 1e-4) :
	if is_parallel(v1, v2, epsilon) :
	 	return (v1.dot(v2) < 0) == opposite		#  方向チェック
	else :
	 	return False


なお、 ε = 5e-4 は、Shade から得た座標値を元に判定する場合のものであり、 script 内で double ベースで組み立てたモデルに対しての判定では、single の有効桁数が確保できる ε =1e-7 として問題ないでしょう。

#3

11 - 2 linked


平行判定の一例として、線形状の inhandle, outhandle の平行状態を返す property xshade.active_shape().control_point(k).linked の判定があります。

scriptで線形状を描くとき、あるアンカーポイントの inhandle とouthandle とが直線状に並んでいると Shade が判断すると、それらのハンドルは自動的にリンクされる仕様になっており、property linked はこの link 状態を返します。


判定は次のようになっています。

     inhandle と outhandle の先端を結ぶ線分 L を基準長とし、
     anchor point と線分 L との距離 r が、

          [ 式 11 - 4 ]

     であれば平行とみなし、linked = True とする


           [ fig 11 - 2 ]


| v1 | = | v2 | の条件下で 11 - 1 で述べた ε = 5e-4 と比較すると、ほぼ同じレベルの判定となっています。


          [ 式 11 - 5 ]


[ script 11 - 4 ]

#  is_linked(inh as vec3, p as vec3,  outh as vec3) as bool
#
#	点 p からのびる inhandle と outhandle が一直線上に並んでいるか判定する

def is_linked(inh, p, outh) :
	u = inh - outh
	v = (inh - p)*(outh - p)
		
	ratio = v.abs()/u.abs2()
	return (ratio <= 0.00025) or (ratio >= 1000)

次の script では linked の返す値と、is_linked( ) で返される判定結果とを比較します。

[ script 11 - 5 ]

import labo
from vec3 import *


#  is_linked(inh as vec3, p as vec3,  outh as vec3) as bool
#
#	点 p からのびる inhandle と outhandle が一直線上に並んでいるか判定する

def is_linked(inh, p, outh) :
	u = inh - outh
	v = (inh - p)*(outh - p)
		
	ratio = v.abs()/u.abs2()
	print 'ratio = ', ratio
	return (ratio <= 0.00025) or (ratio >= 1000)
	
	
#  線形状を作成
def make_line(pL) :
	scene.begin_creating()
	scene.begin_line(None, False)
	for p in pL :
		scene.append_point(p[1], p[0], p[2], None, None)
	scene.end_line()
	scene.end_creating()

	
	
		
scene = xshade.scene()

pL = []
p1 = vec3(-1, 0, 0)
outh1 = vec3(-0.9, 0, 0)
pL.append([None, p1, outh1])

inh = vec3(-0.5, 0, 0)
p = vec3(0, 0, 0)
outh = vec3(0.5, 0, 0)
pL.append([inh, p, outh])

inh2 = vec3(0.9, 0, 0)
p2 = vec3(1, 0, 0)
pL.append([inh2, p2, None])


print ''
p[1] = 0.00025
make_line(pL)
linked = is_linked(inh, p, outh)
print 'linked = ', scene.active_shape().control_point(1).linked	== 1
print 'is_linked = ', linked

print ''
p[1] = 0.0002501
make_line(pL)
linked = is_linked(inh, p, outh)
print 'linked = ', scene.active_shape().control_point(1).linked	== 1	
print 'is_linked = ', linked

[ script 11 - 6 ]
import labo
from vec3 import *


#  is_linked(inh as vec3, p as vec3,  outh as vec3) as bool
#
#	点 p からのびる inhandle と outhandle が一直線上に並んでいるか判定する

def is_linked(inh, p, outh) :
	u = inh - outh
	v = (inh - p)*(outh - p)
		
	ratio = v.abs()/u.abs2()
	print 'ratio = ', ratio
	return (ratio <= 0.00025) or (ratio >= 1000)
	
	
#  線形状を作成
def make_line(pL) :
	scene.begin_creating()
	scene.begin_line(None, False)
	for p in pL :
		scene.append_point(p[1], p[0], p[2], None, None)
	scene.end_line()
	scene.end_creating()

	
	
		
scene = xshade.scene()

pL = []
p1 = vec3(1, 1, 0)
outh1 = vec3(0.9, 1, 0)
pL.append([None, p1, outh1])

inh = vec3(1, 0, 0)
p = vec3(0, 0, 0)
outh = vec3(1, 0, 0)
pL.append([inh, p, outh])

inh2 = vec3(0.9, -1, 0)
p2 = vec3(1, -1, 0)
pL.append([inh2, p2, None])

print ''
outh[1] = 0.0005
inh[1] = -0.0005
make_line(pL)
linked = is_linked(inh, p, outh)
print 'linked = ', scene.active_shape().control_point(1).linked	== 1	
print 'is_linked = ', linked

print ''
outh[1] = 0.00050005
inh[1] = -0.00050005
make_line(pL)
linked = is_linked(inh, p, outh)
print 'linked = ', scene.active_shape().control_point(1).linked	== 1	
print 'is_linked = ', linked


なお、この property linked は control point と handle が作成された時点でセットされ、その後 Shade 上で 移動, 回転, スケールを行っても更新されず、 マニュアルや script で point / handle 座標が変更されたときに更新される仕様になっています。

これにより、本来 linked = True のものが 移動, 回転, スケール処理による丸め誤差の増大で linked = False と判定されないよううまく回避されます。


次の script で、これを確認します。

  • 形状を作成して linked, is_linked をレポート
  • 形状を 移動, 回転, スケーリングして、linked, is_linked をレポート
  • 読み取った handle 座標値をセットし直して、linked, is_linked をレポート

[ script 11 - 7 ]
import labo
from vec3 import *


#  is_linked(inh as vec3, p as vec3,  outh as vec3) as bool
#
#	点 p からのびる inhandle と outhandle が一直線上に並んでいるか判定する

def is_linked(inh, p, outh) :
	u = inh - outh
	v = (inh - p)*(outh - p)
		
	ratio = v.abs()/u.abs2()
	print 'ratio = ', ratio
	return (ratio <= 0.00025) or (ratio >= 1000)
	
	
#  線形状を作成
def make_line(pL) :
	scene.begin_creating()
	scene.begin_line(None, False)
	for p in pL :
		scene.append_point(p[1], p[0], p[2], None, None)
	scene.end_line()
	scene.end_creating()

	
	
		
scene = xshade.scene()

pL = []
p1 = vec3(-1, 0, 0)
outh1 = vec3(-0.9, 0, 0)
pL.append([None, p1, outh1])

inh = vec3(-0.5, 0, 0)
p = vec3(0, 0, 0)
outh = vec3(0.5, 0, 0)
pL.append([inh, p, outh])

inh2 = vec3(0.9, 0, 0)
p2 = vec3(1, 0, 0)
pL.append([inh2, p2, None])

p[1] = 0.00024999
scale = 100



# Shade 上に出力して report
print ''
print 'Shade 上に出力'
make_line(pL)
ob = scene.active_shape()

inh_ = vec3(ob.control_point(1).in_handle)
p_ = vec3(ob.control_point(1).position)
outh_ = vec3(ob.control_point(1).out_handle)
linked_2 = is_linked(inh_, p_, outh_)
print 'linked = ', ob.control_point(1).linked == 1
print 'is_linked = ', linked_2


# Shade 上で 移動, 回転, スケーリング して report
print ''
print 'Shade 上で 移動, 回転, スケーリング'
scene.move_object([0, 0, 0], None, None, [0.5*scale, 1.2*scale, -0.3*scale])
scene.move_object([0, 0, 0], None, [0, -115, 0], None)
scene.move_object([0, 0, 0], None, [40, 0, 0], None)
scene.move_object([0, 0, 0], None, [0, 0, 48], None)
scene.move_object([0, 0, 0], [1.5, 0.7, 0.8], None, None)

inh_ = vec3(ob.control_point(1).in_handle)
p_ = vec3(ob.control_point(1).position)
outh_ = vec3(ob.control_point(1).out_handle)
linked_2 = is_linked(inh_, p_, outh_)
print 'linked = ', ob.control_point(1).linked == 1
print 'is_linked = ', linked_2


# Shade 上で handel 座標値をリセットして report
print ''
print 'Shade 上で handel 座標値をリセット'

ob.control_point(1).in_handle = inh_
ob.control_point(1).out_handle = outh_

inh_ = vec3(ob.control_point(1).in_handle)
p_ = vec3(ob.control_point(1).position)
outh_ = vec3(ob.control_point(1).out_handle)
linked_2 = is_linked(inh_, p_, outh_)
print 'linked = ', ob.control_point(1).linked == 1
print 'is_linked = ', linked_2

#4

11 - 3 要精度管理項目


**11 - 2** で述べたように、linked では 移動, 回転, スケール処理による丸め誤差増大の影響を防ぐようになっているとはいえ、Shade 上での処理を script 内で計算によって行い、その結果を Shade 上に出力すれば、同じ形状であっても linked 状態が異なってしまう可能性を排除できません。

このように A 点から B 点へのベクトルのように座標値の差で与えられる要素に対しては、常に 有効桁数の問題 がつきまといます。

凝った処理を行おうとする場合、オブジェクトのサイズや原点からの位置によって精度に差が生じてしまうことが無視できない問題であるかか否か、注意が必要です。

平行判定は最も身近に頻繁に出て来て、誤判定の影響が大になる可能性もある 要精度管理項目 だと思います。


#5

11 - 4 桁落ちを防ぐ


**関連記事**
  • 09 bezier line 接線ベクトル

これはかなり特殊なケースで、実用上必要とされるケースは極めてレアと思われますが、ハンドルが出てないポイント極近傍での接線ベクトルを計算する場合、[ 式 09 - 1 ] を用いるよりも、[ 式 09 - 4 ] の方が、良い結果が得られます。

[ 式 09 - 1 ] では、p0 ~ p3 に掛かる全ての係数に、桁落ちの原因となる t に関する二乗の計算が含まれてしまうのに対し、[ 式 09 - 2 ] では桁落ちの計算を一カ所に集約して影響が最小限に抑えられています。


          [ 式 09 - 1 ]

          [ 式 09 - 4 ]




[ script 11 - 8 ] では t = 1e-8, 1e-9, 1e-10, 1e-11, 1e-12, 1e-13, 1e-14, 1e-15, 0 での両ケースでの計算結果を比較します。

結果は次のとおりで、[ 式 09 - 1 ] を用いた最初のケースでは、t = 1e-11 ~ 1e-15 で、桁落ちによる累積誤差の影響が顕著に表れています。


t = 1e-8 : v = v [ 0.49805824455, 0.809344648865, -0.311286402512 ]
t = 1e-9 : v = v [ 0.498058230562, 0.809344656093, -0.3112864061 ]
t = 1e-10 : v = v [ 0.498058258399, 0.809344641143, -0.311286400431 ]
t = 1e-11 : v = v [ 0.498059961411, 0.809343728184, -0.311286049301 ]
t = 1e-12 : v = v [ 0.498082668407, 0.809331554936, -0.311281367283 ]
t = 1e-13 : v = v [ 0.497940732835, 0.809407634653, -0.311310628713 ]
t = 1e-14 : v = v [ 0.498224563954, 0.809255467846, -0.311252103018 ]
t = 1e-15 : v = v [ 0.46301474982, 0.827271414546, -0.318181313287 ]
t = 0 : v = v [ 0.498058245092, 0.809344648274, -0.311286403182 ]


t = 1e-8 : v = v [ 0.498058247428, 0.809344647323, -0.311286401918 ]
t = 1e-9 : v = v [ 0.498058245325, 0.809344648179, -0.311286403056 ]
t = 1e-10 : v = v [ 0.498058245115, 0.809344648265, -0.31128640317 ]
t = 1e-11 : v = v [ 0.498058245094, 0.809344648273, -0.311286403181 ]
t = 1e-12 : v = v [ 0.498058245092, 0.809344648274, -0.311286403182 ]
t = 1e-13 : v = v [ 0.498058245092, 0.809344648274, -0.311286403182 ]
t = 1e-14 : v = v [ 0.498058245092, 0.809344648274, -0.311286403182 ]
t = 1e-15 : v = v [ 0.498058245092, 0.809344648274, -0.311286403182 ]
t = 0 : v = v [ 0.498058245092, 0.809344648274, -0.311286403182 ]




[ script 11 - 8 ]

import labo
from vec3 import *


def bezier_line_tangent(bz, t, b) :
	if not b :
		v = (-3*(1 - t)**2)*bz[0] + 3*(1 - t)*(1 - 3*t)*bz[1] + 3*t*(2 - 3*t)*bz[2] + (3*t**2)*bz[3]
	else :
		v = 3*t**2*(-bz[0] + 3*bz[1] - 3*bz[2] + bz[3]) + 6*t*(bz[0] - 2*bz[1] + bz[2]) + 3*(-bz[0] + bz[1])
	v.norm()
	if v.abs2() < 0.5 :
		if t < 0.5 :				#  outhandle が出ていなくて t = 0
			v = bz[2] - bz[0]
		else :					#  inhandle がでていなくて t = 1
			v = bz[3] - bz[1]
		v.norm()
		
	return v
	


bz = [vec3(-8000, 0, 0), vec3(-8000, 0, 0), vec3(0, 13000, -5000), vec3(6000, 7500, 0)]
labo.make_bezier_line(xshade, bz)

for b in [False, True] :
	print ''
	v = bezier_line_tangent(bz, 1e-8, b)
	print 't = 1e-8   : v = ', v
	v = bezier_line_tangent(bz, 1e-9, b)
	print 't = 1e-9   : v = ', v
	v = bezier_line_tangent(bz, 1e-10, b)
	print 't = 1e-10   : v = ', v
	v = bezier_line_tangent(bz, 1e-11, b)
	print 't = 1e-11   : v = ', v
	v = bezier_line_tangent(bz, 1e-12, b)
	print 't = 1e-12   : v = ', v
	v = bezier_line_tangent(bz, 1e-13, b)
	print 't = 1e-13   : v = ', v
	v = bezier_line_tangent(bz, 1e-14, b)
	print 't = 1e-14   : v = ', v
	v = bezier_line_tangent(bz, 1e-15, b)
	print 't = 1e-15   : v = ', v
	v = bezier_line_tangent(bz, 0, b)
	print 't = 0   : v = ', v