前回の最後で パックマンもどき(その3) では、なぜだか円柱の表示が消えてる事が確認されました。
そんなバカな。 でも、cross.c をv58の状態に戻すと円柱は表示されます。 ということは原因は v59.patch の変更
$ cat v59.patch
diff -urN rt_v58/cross.c rt_v59/cross.c
--- rt_v58/cross.c 2018-06-29 09:20:42.000000000 +0900
+++ rt_v59/cross.c 2018-07-06 03:39:31.000000000 +0900
@@ -334,19 +334,27 @@
if(kind_idx != KIND_BALL){
double *p = ret->l_nv_g.p;
double *v = ret->l_nv_g.v;
- double z[3];
+ double x[3], z[3];
struct vecs vs;
if(is_plane(kind_idx)){
v_set(z, v_x1);
}else if(kind_idx == KIND_PIPE_SIDE){
- v_add(p, v_z1, z);
- }else{ /* KIND_CONE_SIDE */
v_set(z, v_z1);
+ }else{ /* KIND_CONE_SIDE */
+ v_sub(v_z1, ret->p, z);
}
+ vecs_by_fix_ref(&vs, 'z', z, 'y', nv);
+ v_set( x, vecs_get_by_name(&vs, 'x') );
+
+ v_add(x, ret->p, x);
+ v_add(z, ret->p, z);
+ lstx_tr_p(l2g, L2G, x, x);
lstx_tr_p(l2g, L2G, z, z);
+ v_sub(x, p, x);
v_sub(z, p, z);
- vecs_by_fix_ref(&vs, 'z', z, 'y', v);
+ vecs_by_fix_ref(&vs, 'x', x, 'z', z);
+ v_set( v, vecs_get_by_name(&vs, 'y') );
}
if(rev){
最後のv_set()をコメントアウトしてvの値を更新しなければ、表示は出ます。 修正したつもりの法線ベクトルの処理に、何かまちがいが...
それでも、円柱の場合のKIND_PIPE_SIDEの処理を確認しても、問題なさげです。 なぜに円柱が消える?
デバッグアウト出して追い込んでいくと、見えてきました。
平面の法線ベクトルが、逆向きになってる。 というか、逆向きになる事がある。
円柱は、無限の長さの筒と、2つの円の平面の「AND」で構成されます。 上下の円の法線が逆を向いていれば、その時点で中身は空になり、確かに円柱は消えます。
法線がなぜ逆向きになる事があるか? 先の cross.c cross_get_nv() の処理
if(kind_idx != KIND_BALL){
double *p = ret->l_nv_g.p;
double *v = ret->l_nv_g.v;
double x[3], z[3];
struct vecs vs;
if(is_plane(kind_idx)){
v_set(z, v_x1);
}else if(kind_idx == KIND_PIPE_SIDE){
v_set(z, v_z1);
}else{ /* KIND_CONE_SIDE */
v_sub(v_z1, ret->p, z);
}
// ローカル座標系で算出した法線を 'y'
// それに明らかに直行する 'z' を用意しておいて
vecs_by_fix_ref(&vs, 'z', z, 'y', nv);
v_set( x, vecs_get_by_name(&vs, 'x') );
// 'y', 'z' に直行する 'x' を算出
v_add(x, ret->p, x);
v_add(z, ret->p, z);
// ローカル座標系での交点 ret->p から
// 'x', 'z' の方向に位置する2点を用意
lstx_tr_p(l2g, L2G, x, x);
lstx_tr_p(l2g, L2G, z, z);
// 2点をグローバル座標系に変換して
v_sub(x, p, x);
v_sub(z, p, z);
// グローバル座標系での交点 p を引き算して
// 2点を、2つのグローバル座標系のベクトル 'x', 'z' に
vecs_by_fix_ref(&vs, 'x', x, 'z', z);
v_set( v, vecs_get_by_name(&vs, 'y') );
// 2つのグローバル座標系のベクトル 'x', 'z' に直行する
// グローバル座標系のベクトル 'y' を算出して、v を更新
}
な訳ですが、ローカル座標系とグローバル座標系。 アフィン変換で変換してますが、回転だけじゃなく拡大縮小もあるので、例えば、X方向に-1倍とかもあります。
という事は、ローカル座標系とグローバル座標系で、右手系と左手系の違いが、、、あり得ます。
という事はという事は、平面に限らず円柱も円すいも法線ベクトルが逆向きになってる場合もあったかも。
いや〜、法線が逆向きでも、何とか表示が出てると判りにくいもんですね。 修正します。
$ mv rt_v60 rt_v61 $ cat v61.patch | ( cd rt_v61 ; patch -p1 ) $ cd rt_v61 $ make clean $ make
それでは、またしてもリテイク
$ ./pac_dat.py $ ( cd pac ; wget http://kondoh2.html.xdomain.jp/rt/pac.log ) $ ./pac_dat.py $ ./cg.py eyep=[0,0,20],[50,50,10],20 sec=30 yaml=pac.yaml data_name=all name=out_v61/pac_tst3_ fps=0.2 : wh : 307200/307200(100.0%) : fin 50.53s frm : 6/6(100.0%) : fin 5m 13.18s estimated 13.05 hour at 640*480 30fps $ ls -lt out_v61/ | head -rw-r--r-- 1 kondoh staff 144 7 13 20:22 pac_tst3_.mp4 -rw-r--r-- 1 kondoh staff 136861 7 13 20:22 pac_tst3_00006.jpg -rw-r--r-- 1 kondoh staff 123617 7 13 20:21 pac_tst3_00005.jpg -rw-r--r-- 1 kondoh staff 133281 7 13 20:20 pac_tst3_00004.jpg -rw-r--r-- 1 kondoh staff 161416 7 13 20:19 pac_tst3_00003.jpg -rw-r--r-- 1 kondoh staff 88597 7 13 20:19 pac_tst3_00002.jpg :
法線ベクトルを修正して円柱が表示されるようになりました。
が、が、ヘアライン!
ヘアライン(その2) で修正したはずなのに、またしても!
まずは、例のカットアンドトライで切り出した再現データで試してみます。
$ cat hair_line.yaml
all:
- kind: or
args:
- kind: pipe
rtd: { base: 0.1, diff: 0.3, reflect: 0.5, reflact: 0.2, density: 2 }
m2g:
- ax.rot_y(90)
- ax.slide_y(2)
# EOF
$ ./cg.py eyep=[0,-50,20],0,1 yaml=hair_line.yaml data_name=all name=out_v61/h n=1 show_sec=-1 div=2
確かに、このデータではヘアラインの現象は再現しません。
とりあえず再現するpac.yamlをhh.yamlとしてコピーして、 手動でinclude行を展開しながら、 削っていって追い込んでみます。
$ cat hh.yaml
blk_h:
- kind: pipe
rtd: { base: 0.1, diff: 0.3, reflect: 0.5, reflact: 0.2, density: 2 }
m2g:
- ax.slide_z(1)
- 'ax.zoom([ 0.2 , 0.2 , 0.5 ])'
- ax.rot_y(90)
blk_hs:
- export: blk_h
kind: export
m2g:
- ax.zoom_x(2)
- ax.slide([12,8,0])
blk:
- kind: or
args:
- kind: exports
exports:
- blk_hs
def_col: [200,50,50]
pac:
- kind: exports
exports: [ blk ]
m2g:
- 'ax.slide([-21.0/2,-23.0/2,0])'
- ax.zoom_y(-1)
all:
- kind: exports
exports: [ pac ]
# EOF
$ ./cg.py eyep=[0,0,20],[50,50,10],20 sec=30 yaml=hh.yaml data_name=all name=out_v61/hh n=1 show_sec=-1 init_sec=5
:
さらに追い込んで
$ cat hh.haml all: - kind: or args: - kind: pipe m2g: - ax.slide_z(1) - 'ax.zoom([ 0.2 , 0.2 , 0.5 ])' - ax.rot_y(90) - ax.zoom_x(2) - ax.slide([12,8,0]) - 'ax.slide([-21.0/2,-23.0/2,0])' - ax.zoom_y(-1) # EOF $ ./cg.py eyep=[0,0,20],[50,50,10],20 sec=30 yaml=hh.yaml data_name=all name=out_v61/hh2_ n=1 show_sec=-1 init_sec=5 :
ここから
$ cat hh.yaml all: #- kind: or # args: - kind: pipe m2g: - ax.slide_z(1) - 'ax.zoom([ 0.2 , 0.2 , 0.5 ])' - ax.rot_y(90) - ax.zoom_x(2) - ax.slide([12,8,0]) - 'ax.slide([-21.0/2,-23.0/2,0])' - ax.zoom_y(-1) # EOF $ ./cg.py eyep=[0,0,20],[50,50,10],20 sec=30 yaml=hh.yaml data_name=all name=out_v61/hh3_ n=1 show_sec=-1 init_sec=5 :
とすると、再現しません。
$ cat dat.py : if kind == 'pipe': f = for_op() d0 = d_copy(d) d0['kind'] = 'pipe_side' d1 = d_copy(d) d1['kind'] = 'circle' d1['l2m'] = [ 'ax.slide_z(1)' ] + d1.get( 'l2m', [] ) d2 = d_copy(d) d2['kind'] = 'circle' d2['l2m'] = [ 'ax.slide_z(1)', 'ax.rot_y(180)' ] + d2.get( 'l2m', [] ) return ret_lst( [d0,d1,d2], f ) :
なので、pipe は pipe_side と2つの circle の AND に展開されます。
$ cat hh.haml
all:
- kind: or
args:
- kind: and
args:
- kind: pipe_side
- kind: circle
l2m:
- ax.slide_z(1)
- kind: circle
l2m:
- ax.slide_z(1)
- ax.rot_y(180)
m2g:
- ax.slide_z(1)
- 'ax.zoom([ 0.2 , 0.2 , 0.5 ])'
- ax.rot_y(90)
- ax.zoom_x(2)
- ax.slide([12,8,0])
- 'ax.slide([-21.0/2,-23.0/2,0])'
- ax.zoom_y(-1)
# EOF
$ ./cg.py eyep=[0,0,20],[50,50,10],20 sec=30 yaml=hh.yaml data_name=all name=out_v61/hh4_ n=1 show_sec=-1 init_sec=5
これでも再現します。
$ cat hh.yaml
all:
- kind: and
args:
- kind: pipe_side
- kind: circle
l2m:
- ax.slide_z(1)
- kind: circle
l2m:
- ax.slide_z(1)
- ax.rot_y(180)
m2g:
- ax.slide_z(1)
- 'ax.zoom([ 0.2 , 0.2 , 0.5 ])'
- ax.rot_y(90)
- ax.zoom_x(2)
- ax.slide([12,8,0])
- 'ax.slide([-21.0/2,-23.0/2,0])'
- ax.zoom_y(-1)
# EOF
そして、pipeを手動で展開した後なら、これでも再現します。
デバッグアウトで位置をば rt.c set_cols()
$ cat rt.c
:
for(i=0; i<n; i++){
lp[0] = (idx + i) % w;
lp[1] = (idx + i) / w;
lstx_tr_p(wh2g, L2G, lp, p);
line_new_p2(&l_g, eyep, p);
get_col(&l_g, sec, 1.0, -1, col);
if(!is_no_col(col)){
col_max_255(col, col);
#if 1
dbg_enable = 1;
DBG("%d %d: ", (int)lp[0], (int)lp[1]);
DBG("%d %d %d\n", col[0], col[1], col[2]);
#endif
}
col += 3;
}
:
$ cat rt.py
:
names = [ 'srv.rt', 'srv.rt-15', 'srv.rt-16', 'srv.rt-17', 'srv.rt-18' ][:2]
#names += [ 'srv.rt-2', 'srv.rt-3', 'srv.rt-4', 'srv.rt-5', 'srv.rt-6' ][:2]
names = [ 'srv.rt' ] * 2 + [ 'srv.rt-15' ] * 2
names = [ 'srv.rt' ]
:
# デバッグアウトが入り乱れないように、スレッドを1つにして...
$ ./skill.sh
$ make
$ ./cg.py eyep=[0,0,20],[50,50,10],20 sec=30 yaml=hh.yaml data_name=all name=out_v61/hh5_ n=1 init_sec=5
:
638 36: 85 85 85
639 36: 85 85 85
320 37: 74 74 74
457 37: 216 216 216
458 37: 216 216 216
459 37: 73 73 73
460 37: 73 73 73
461 37: 73 73 73
:
638 37: 72 72 72
639 37: 72 72 72
320 38: 63 63 63
457 38: 216 216 216
458 38: 216 216 216
459 38: 62 62 62
460 38: 62 62 62
461 38: 62 62 62
:
638 38: 62 62 62
639 38: 62 62 62
320 39: 63 63 63
457 39: 216 216 216
458 39: 216 216 216
459 39: 63 63 63
460 39: 63 63 63
461 39: 63 63 63
:
ぱっと見てもヘアラインの位置はxが中央の辺りなので x=320 で確定ですね。
$ cat rt.c
:
lp[2] = -d_;
for(i=0; i<n; i++){
lp[0] = (idx + i) % w;
lp[1] = (idx + i) / w;
lstx_tr_p(wh2g, L2G, lp, p);
line_new_p2(&l_g, eyep, p);
#if 1
dbg_enable = ( (int)lp[0] == 320 && (int)lp[1] == 37 );
#endif
get_col(&l_g, sec, 1.0, -1, col);
if(!is_no_col(col)){
col_max_255(col, col);
}
col += 3;
}
:
$ cat cross.c
:
cross_op_ins(struct line *l_g, int prev_idx, struct data *data, struct area *a, struct cross_ret *r)
{
struct data_one *d = &data->lst[r->idx];
int n = r->ts_n;
area_init(a, 0);
#if 1
DBG("kind=%d n=%d\n", d->kind, n);
#endif
if(n == 0){
:
$ ./skill.sh
$ make
$ ./cg.py eyep=[0,0,20],[50,50,10],20 sec=30 yaml=hh.yaml data_name=all name=out_v61/hh5_ n=1 init_sec=5
:
kind=4 n=2
kind=1 n=0
kind=1 n=1
:
ふむ。
$ cat cross.c
:
area_init(a, 0);
#if 1
DBG("kind=%d n=%d\n", d->kind, n);
#endif
if(n == 0){
if(is_plane(d->kind) && r->l.p[2] < 0){
area_not(a);
}
}else if(n == 1){
r->t = r->ts[0];
line_on_line_p( &r->l, r->t, r->p );
cross_get_nv( d->kind, &d->l2g, l_g, 0, r );
#if 1
DBG("and_nv_eyev=%f =(%d) >(%d) t=%f\n", r->ang_nv_eyev, (r->ang_nv_eyev == 0), (r->ang_nv_eyev > 0), r->t );
#endif
if(r->ang_nv_eyev == 0){
r->ts_n = 0;
cross_op_ins(l_g, prev_idx, data, a, r);
return;
}
r->flags = 0;
if( prev_idx == r->idx ){
r->flags |= CROSS_FLAG_PREV;
}
area_ins( a, r->t, r, sizeof(*r) );
if(r->ang_nv_eyev > 0){
area_not(a);
}
:
$ ./skill.sh
$ make
$ ./cg.py eyep=[0,0,20],[50,50,10],20 sec=30 yaml=hh.yaml data_name=all name=out_v61/hh5_ n=1 init_sec=5
:
kind=4 n=2
kind=1 n=0
kind=1 n=1
and_nv_eyev=nan =(0) >(0) t=-6755399441055697.000000
:
かなり視点の後方側(tが負)の位置に交点があって 視線のベクトルと法線ベクトルの内積が無限大?
$ cat cross.c
:
cross_get_nv(int kind_idx, struct lstx *l2g, struct line *l_g, int rev, struct cross_ret *ret)
{
:
#if 1
DBG("nv=[%f %f %f] eyev=[%f %f %f] dot=%f\n",
ret->nv[0],ret->nv[1],ret->nv[2],
ret->eyev[0],ret->eyev[1],ret->eyev[2],
ret->ang_nv_eyev);
#endif
}
:
$ ./skill.sh
$ make
$ ./cg.py eyep=[0,0,20],[50,50,10],20 sec=30 yaml=hh.yaml data_name=all name=out_v61/hh5_ n=1 init_sec=5
:
kind=4 n=2
kind=1 n=0
kind=1 n=1
nv=[nan nan nan] eyev=[0.000000 0.837759 -0.546039] dot=nan
and_nv_eyev=nan =(0) >(0) t=-6755399441055697.000000
:
法線ベクトルの各成分が無限大
ということは、その直前のここ
: double *v = ret->l_nv_g.v; : v_set(v, y); : unit(ret->l_nv_g.v, ret->nv);
vが(0,0,0)で単位ベクトルにしようとして、無限大。 ということは、やっぱり法線ベクトル処理か。
if(kind_idx != KIND_BALL){
double *p = ret->l_nv_g.p;
double *v = ret->l_nv_g.v;
double x[3], z[3];
struct vecs vs;
if(is_plane(kind_idx)){
v_set(z, v_x1);
}else if(kind_idx == KIND_PIPE_SIDE){
v_set(z, v_z1);
}else{ /* KIND_CONE_SIDE */
v_sub(v_z1, ret->p, z);
}
// ローカル座標系で算出した法線を 'y'
// それに明らかに直行する 'z' を用意しておいて
vecs_by_fix_ref(&vs, 'z', z, 'y', nv);
v_set( x, vecs_get_by_name(&vs, 'x') );
// 'y', 'z' に直行する 'x' を算出
v_add(x, ret->p, x);
v_add(z, ret->p, z);
// ローカル座標系での交点 ret->p から
// 'x', 'z' の方向に位置する2点を用意
lstx_tr_p(l2g, L2G, x, x);
lstx_tr_p(l2g, L2G, z, z);
// 2点をグローバル座標系に変換して
v_sub(x, p, x);
v_sub(z, p, z);
// グローバル座標系での交点 p を引き算して
// 2点を、2つのグローバル座標系のベクトル 'x', 'z' に
vecs_by_fix_ref(&vs, 'x', x, 'z', z);
v_set( v, vecs_get_by_name(&vs, 'y') );
// 2つのグローバル座標系のベクトル 'x', 'z' に直行する
// グローバル座標系のベクトル 'y' を算出して、v を更新
}
とりあえずデバッグアウト追加
$ cat cross.c
:
if(kind_idx != KIND_BALL){
double *p = ret->l_nv_g.p;
double *v = ret->l_nv_g.v;
double x[3], y[3], z[3];
struct vecs vs;
if(is_plane(kind_idx)){
v_set(z, v_x1);
}else if(kind_idx == KIND_PIPE_SIDE){
v_set(z, v_z1);
}else{ /* KIND_CONE_SIDE */
v_sub(v_z1, ret->p, z);
}
vecs_by_fix_ref(&vs, 'z', z, 'y', nv);
v_set( x, vecs_get_by_name(&vs, 'x') );
#if 1
DBG("kind=%d z=%f %f %f nv=%f %f %f x=%f %f %f\n",
kind_idx, z[0],z[1],z[2],nv[0],nv[1],nv[2],x[0],x[1],x[2]);
#endif
v_add(x, ret->p, x);
v_add(z, ret->p, z);
lstx_tr_p(l2g, L2G, x, x);
lstx_tr_p(l2g, L2G, z, z);
v_sub(x, p, x);
v_sub(z, p, z);
#if 1
DBG("x=%f %f %f z=%f %f %f\n",
x[0],x[1],x[2],z[0],z[1],z[2]);
#endif
vecs_by_fix_ref(&vs, 'x', x, 'z', z);
v_set( y, vecs_get_by_name(&vs, 'y') );
if(dot(y, v) < 0){
v_neg(y, y);
}
v_set(v, y);
}
:
$ ./skill.sh
$ make
$ ./cg.py eyep=[0,0,20],[50,50,10],20 sec=30 yaml=hh.yaml data_name=all name=out_v61/hh5_ n=1 init_sec=5
:
kind=1 n=1
kind=1 z=1.000000 0.000000 0.000000 nv=0.000000 0.000000 1.000000 x=0.000000 1.000000 0.000000
x=0.000000 0.000000 0.000000 z=0.000000 0.000000 0.000000
nv=[nan nan nan] eyev=[0.000000 0.837759 -0.546039] dot=nan
and_nv_eyev=nan =(0) >(0) t=-6755399441055697.000000
:
xのベクトルをローカル座標系からグローバル座標系に変換すると(0,0,0)
: v_add(x, ret->p, x); v_add(z, ret->p, z); lstx_tr_p(l2g, L2G, x, x); lstx_tr_p(l2g, L2G, z, z); v_sub(x, p, x); v_sub(z, p, z); :
$ cat cross.c
:
if(kind_idx != KIND_BALL){
double *p = ret->l_nv_g.p;
double *v = ret->l_nv_g.v;
double x[3], y[3], z[3];
struct vecs vs;
if(is_plane(kind_idx)){
v_set(z, v_x1);
}else if(kind_idx == KIND_PIPE_SIDE){
v_set(z, v_z1);
}else{ /* KIND_CONE_SIDE */
v_sub(v_z1, ret->p, z);
}
vecs_by_fix_ref(&vs, 'z', z, 'y', nv);
v_set( x, vecs_get_by_name(&vs, 'x') );
#if 1
DBG("kind=%d z=%f %f %f nv=%f %f %f x=%f %f %f\n",
kind_idx, z[0],z[1],z[2],nv[0],nv[1],nv[2],x[0],x[1],x[2]);
#endif
v_add(x, ret->p, x);
#if 1
DBG("ret->p=%f %f %f x=%f %f %f\n",
ret->p[0],ret->p[1],ret->p[2],x[0],x[1],x[2]);
#endif
v_add(z, ret->p, z);
lstx_tr_p(l2g, L2G, x, x);
#if 1
DBG("tr x=%f %f %f\n", x[0],x[1],x[2]);
#endif
lstx_tr_p(l2g, L2G, z, z);
v_sub(x, p, x);
#if 1
DBG("p=%f %f %f x=%f %f %f\n",
p[0],p[1],p[2],x[0],x[1],x[2]);
#endif
v_sub(z, p, z);
#if 1
DBG("x=%f %f %f z=%f %f %f\n",
x[0],x[1],x[2],z[0],z[1],z[2]);
#endif
vecs_by_fix_ref(&vs, 'x', x, 'z', z);
v_set( y, vecs_get_by_name(&vs, 'y') );
if(dot(y, v) < 0){
v_neg(y, y);
}
v_set(v, y);
}
:
$ ./skill.sh
$ make
$ ./cg.py eyep=[0,0,20],[50,50,10],20 sec=30 yaml=hh.yaml data_name=all name=out_v61/hh5_ n=1 init_sec=5
:
kind=1 n=1
kind=1 z=1.000000 0.000000 0.000000 nv=0.000000 0.000000 1.000000 x=0.000000 1.000000 0.000000
ret->p=18443572058546440.000000 28296999662978408.000000 0.000000 x=18443572058546440.000000 28296999662978408.000000 0.000000
tr x=-1.210423 -5659399932595678.000000 3688714411709288.000000
p=-1.210423 -5659399932595678.000000 3688714411709288.000000 x=0.000000 0.000000 0.000000
x=0.000000 0.000000 0.000000 z=0.000000 0.000000 0.000000
nv=[nan nan nan] eyev=[0.000000 0.837759 -0.546039] dot=nan
and_nv_eyev=nan =(0) >(0) t=-6755399441055697.000000
:
x=(0,1,0) に ret->p=( 大きな値a, 大きな値b, 0 ) を足すと x=( 大きな値a, 大きな値b, 0 ) で、y成分が1増えてない。
pythonの4バイトのfloatの精度の問題!?
見えてきました。 交点と、交点から(0,1,0)進めた位置が交点と同じに扱われて、、、 その同じ2点をグローバル座標系の位置に変換しても、同じ位置。 2つの同じ位置を引き算してベクトルを出すと 0 ベクトル。
成分に1足しても変化しないほど離れた位置に交点がある場合、 それはもう「交差してなくて並行な扱い」としていいのではなかろうか。
v62 へつづく。