投稿者:しばっち
投稿日:2018年 3月10日(土)21時15分55秒
|
|
|
ここではParaViewのボリュームレンダリング機能を使用します。
格子状にデータを書き込んでいきます。
3次元陰関数 f(x,y,z)=0 を定義し、その内側なら255(※この値に意味はありません)
外側なら0を書き込んでいます(2値データ)
プログラムでは3次元陰関数に2次曲面を定義しています。
また、3次元空間での探索なので実行に時間がかかります。
生成されるファイルサイズも大きくなります。
表示方法の「Surface」から「Volume」にしてください。
ボリュームレンダリング表示されます。
更に3行目電卓アイコンの右隣、半球みたいなアイコン「Contour」のボタンを押してください。
又は「Filters」メニューの「Alphabetical」から「Contour」を選択してください。
meshデータが自動生成されて「Surface」で表示されます。
重なって表示されるときは、Pipeline Browser内の目のようなアイコンをクリックして非表示にしてください。
マシンパワーによって各格子のサイズを調整してください。
RANDOMIZE
LET XSIZE=100 !'各格子の大きさ
LET YSIZE=100
LET ZSIZE=100
LET A=INT(RND*20)-10
LET B=INT(RND*20)-10
LET C=INT(RND*20)-10
LET D=INT(RND*20)-10
LET E=INT(RND*20)-10
LET F=INT(RND*20)-10
LET G=INT(RND*20)-10
LET XS=-5 !'探索範囲
LET XE=5
LET YS=-5
LET YE=5
LET ZS=-5
LET ZE=5
FILE GETSAVENAME F$,"vtkファイル|*.vtk"
IF F$="" THEN STOP
IF POS(UCASE$(F$),".VTK")=0 THEN LET F$=F$&".vtk"
OPEN #1:NAME F$
ERASE #1
PRINT #1:"# vtk DataFile Version 2.0"
PRINT #1:"inequality"
PRINT #1:"ASCII"
PRINT #1:"DATASET STRUCTURED_POINTS"
PRINT #1:"DIMENSIONS";XSIZE;YSIZE;ZSIZE
PRINT #1:"ORIGIN 0 0 0"
PRINT #1:"SPACING";(XE-XS)/XSIZE;(YE-YS)/YSIZE;(ZE-ZS)/ZSIZE
PRINT #1:"POINT_DATA";XSIZE*YSIZE*ZSIZE
PRINT #1:"SCALARS scalar float"
PRINT #1:"LOOKUP_TABLE default"
FOR ZZ=0 TO ZSIZE-1
FOR YY=0 TO YSIZE-1
FOR XX=0 TO XSIZE-1
LET X=XS+XX*(XE-XS)/XSIZE
LET Y=YS+YY*(YE-YS)/YSIZE
LET Z=ZS+ZZ*(ZE-ZS)/ZSIZE
IF X*X+Y*Y+Z*Z+A*X*Y+B*X*Z+C*Y*Z+D*X+E*Y+F*Z+G<0 THEN LET DA=255 ELSE LET DA=0
PRINT #1:DA
NEXT XX
NEXT YY
NEXT ZZ
CLOSE #1
END
|
|
|
投稿者:しばっち
投稿日:2018年 3月21日(水)12時58分17秒
|
|
|
> No.4530[元記事へ]
GLSL言語(シェーディング言語)のレイ・マーチング法(レイ・トレーシング法の一種)で使用される距離関数のいくつかを
十進BASICに移植してみました。
ここではParaViewのボリュームレンダリング機能を使用し、距離関数によって3D形状を定義しています。
ParaViewでの表示を「volume」にしてください。
距離関数で角が丸い立方体を定義しています。
LET XSIZE=100 !'各格子の大きさ
LET YSIZE=100
LET ZSIZE=100
DIM P(3),B(3)
CALL VEC3(B,.7,.7,.7)
FILE GETSAVENAME F$,"vtkファイル|*.vtk"
IF F$="" THEN STOP
IF POS(UCASE$(F$),".VTK")=0 THEN LET F$=F$&".vtk"
OPEN #1:NAME F$
ERASE #1
PRINT #1:"# vtk DataFile Version 2.0"
PRINT #1:"round box"
PRINT #1:"ASCII"
PRINT #1:"DATASET STRUCTURED_POINTS"
PRINT #1:"DIMENSIONS";XSIZE;YSIZE;ZSIZE
PRINT #1:"ORIGIN 0 0 0"
PRINT #1:"SPACING 1 1 1"
PRINT #1:"POINT_DATA";XSIZE*YSIZE*ZSIZE
PRINT #1:"SCALARS scalar float"
PRINT #1:"LOOKUP_TABLE default"
FOR ZZ=0 TO ZSIZE-1
FOR YY=0 TO YSIZE-1
FOR XX=0 TO XSIZE-1
LET X=-1.5+3/XSIZE*XX !'-1.5 ~ 1.5 探索範囲
LET Y=-1.5+3/YSIZE*YY !'-1.5 ~ 1.5
LET Z=-1.5+3/ZSIZE*ZZ !'-1.5 ~ 1.5
CALL VEC3(P,X,Y,Z)
IF ROUNDBOX(P,B,.4)<.1 THEN LET DA=255 ELSE LET DA=0 !'距離関数による判定
PRINT #1:DA
NEXT XX
NEXT YY
NEXT ZZ
CLOSE #1
END
EXTERNAL FUNCTION ROUNDBOX(P(),B(),R) !'float udRoundBox( vec3 p, vec3 b, float r )
DIM V(3),A(3)
CALL VABS(V,P)
MAT V=V-B
CALL NMAX(A,V,0)
LET ROUNDBOX=LENGTH(A)-R !' { return length(max(abs(p)-b,0.0))-r;}
END FUNCTION
EXTERNAL SUB VABS(A(),B())
FOR I=1 TO 3
LET A(I)=ABS(B(I))
NEXT I
END SUB
EXTERNAL FUNCTION LENGTH(A())
LET LENGTH=SQR(A(1)^2+A(2)^2+A(3)^2)
END FUNCTION
EXTERNAL SUB NMAX(A(),B(),N)
FOR I=1 TO 3
LET A(I)=MAX(B(I),N)
NEXT I
END SUB
EXTERNAL SUB VEC3(A(),X,Y,Z)
LET A(1)=X
LET A(2)=Y
LET A(3)=Z
END SUB
|
|
|
投稿者:しばっち
投稿日:2018年 3月21日(水)12時59分2秒
|
|
|
> No.4534[元記事へ]
ParaViewのボリュームレンダリング機能を使用しています。
表示を「volume」にしてください。
楕円体を定義しています。
マシンパワーによって各格子のサイズを調整してください。
LET XSIZE=100 !'各格子の大きさ
LET YSIZE=100
LET ZSIZE=100
DIM P(3),R(3)
CALL VEC3(R,.2,.3,.4)
FILE GETSAVENAME F$,"vtkファイル|*.vtk"
IF F$="" THEN STOP
IF POS(UCASE$(F$),".VTK")=0 THEN LET F$=F$&".vtk"
OPEN #1:NAME F$
ERASE #1
PRINT #1:"# vtk DataFile Version 2.0"
PRINT #1:"ellipsoid"
PRINT #1:"ASCII"
PRINT #1:"DATASET STRUCTURED_POINTS"
PRINT #1:"DIMENSIONS";XSIZE;YSIZE;ZSIZE
PRINT #1:"ORIGIN 0 0 0"
PRINT #1:"SPACING 1 1 1"
PRINT #1:"POINT_DATA";XSIZE*YSIZE*ZSIZE
PRINT #1:"SCALARS scalar float"
PRINT #1:"LOOKUP_TABLE default"
FOR ZZ=0 TO ZSIZE-1
FOR YY=0 TO YSIZE-1
FOR XX=0 TO XSIZE-1
LET X=-1+2/XSIZE*XX !'-1 ~ 1 探索範囲
LET Y=-1+2/YSIZE*YY !'-1 ~ 1
LET Z=-1+2/ZSIZE*ZZ !'-1 ~ 1
CALL VEC3(P,X,Y,Z)
IF ELLIPSOID(P,R)<.1 THEN LET DA=255 ELSE LET DA=0
PRINT #1:DA
NEXT XX
NEXT YY
NEXT ZZ
CLOSE #1
END
EXTERNAL FUNCTION ELLIPSOID(P(),R()) !'float sdEllipsoid( in vec3 p, in vec3 r ){
LET ELLIPSOID=(LENGTH2(P(1)/R(1),P(2)/R(2),P(3)/R(3))-1)*MIN(MIN(R(1),R(2)),R(3)) !'return (length( p/r ) - 1.0) * min(min(r.x,r.y),r.z);}
END FUNCTION
EXTERNAL FUNCTION LENGTH2(A,B,C)
LET LENGTH2=SQR(A^2+B^2+C^2)
END FUNCTION
EXTERNAL SUB VEC3(A(),X,Y,Z)
LET A(1)=X
LET A(2)=Y
LET A(3)=Z
END SUB
|
|
|
投稿者:しばっち
投稿日:2018年 3月21日(水)13時00分18秒
|
|
|
> No.4534[元記事へ]
ParaViewのボリュームレンダリング機能を使用しています。
表示を「volume」にしてください。
トーラス(ドーナツ型)を定義しています。
マシンパワーによって各格子のサイズを調整してください。
LET XSIZE=100
LET YSIZE=100
LET ZSIZE=100
DIM P(3),T(3)
CALL VEC2(T,.7,.1)
FILE GETSAVENAME F$,"vtkファイル|*.vtk"
IF F$="" THEN STOP
IF POS(UCASE$(F$),".VTK")=0 THEN LET F$=F$&".vtk"
OPEN #1:NAME F$
ERASE #1
PRINT #1:"# vtk DataFile Version 2.0"
PRINT #1:"torus"
PRINT #1:"ASCII"
PRINT #1:"DATASET STRUCTURED_POINTS"
PRINT #1:"DIMENSIONS";XSIZE;YSIZE;ZSIZE
PRINT #1:"ORIGIN 0 0 0"
PRINT #1:"SPACING 1 1 1"
PRINT #1:"POINT_DATA";XSIZE*YSIZE*ZSIZE
PRINT #1:"SCALARS scalar float"
PRINT #1:"LOOKUP_TABLE default"
FOR ZZ=0 TO ZSIZE-1
FOR YY=0 TO YSIZE-1
FOR XX=0 TO XSIZE-1
LET X=-1+2/XSIZE*XX
LET Y=-1+2/YSIZE*YY
LET Z=-1+2/ZSIZE*ZZ
CALL VEC3(P,X,Y,Z)
IF TORUS(P,T)<.1 THEN LET DA=255 ELSE LET DA=0
PRINT #1:DA
NEXT XX
NEXT YY
NEXT ZZ
CLOSE #1
END
EXTERNAL FUNCTION TORUS(P(),T()) !'float sdTorus( vec3 p, vec2 t ){
DIM Q(3)
CALL VEC2(Q,LENGTH2(P(1),P(3),0)-T(1),P(2)) !' vec2 q = vec2(length(p.xz)-t.x,p.y);
LET TORUS=LENGTH(Q)-T(2) !' return length(q)-t.y;}
END FUNCTION
EXTERNAL FUNCTION LENGTH2(A,B,C)
LET LENGTH2=SQR(A^2+B^2+C^2)
END FUNCTION
EXTERNAL FUNCTION LENGTH(A())
LET LENGTH=SQR(A(1)^2+A(2)^2+A(3)^2)
END FUNCTION
EXTERNAL SUB VEC2(A(),X,Y)
LET A(1)=X
LET A(2)=Y
END SUB
EXTERNAL SUB VEC3(A(),X,Y,Z)
LET A(1)=X
LET A(2)=Y
LET A(3)=Z
END SUB
|
|
|
投稿者:しばっち
投稿日:2018年 3月21日(水)13時01分7秒
|
|
|
> No.4534[元記事へ]
ParaViewのボリュームレンダリング機能を使用しています。
表示を「volume」にしてください。
6角柱を定義しています。
マシンパワーによって各格子のサイズを調整してください。
LET XSIZE=100
LET YSIZE=100
LET ZSIZE=100
DIM H(3),P(3)
CALL VEC2(H,1,.5)
FILE GETSAVENAME F$,"vtkファイル|*.vtk"
IF F$="" THEN STOP
IF POS(UCASE$(F$),".VTK")=0 THEN LET F$=F$&".vtk"
OPEN #1:NAME F$
ERASE #1
PRINT #1:"# vtk DataFile Version 2.0"
PRINT #1:"Hex Prism"
PRINT #1:"ASCII"
PRINT #1:"DATASET STRUCTURED_POINTS"
PRINT #1:"DIMENSIONS";XSIZE;YSIZE;ZSIZE
PRINT #1:"ORIGIN 0 0 0"
PRINT #1:"SPACING 1 1 1"
PRINT #1:"POINT_DATA";XSIZE*YSIZE*ZSIZE
PRINT #1:"SCALARS scalar float"
PRINT #1:"LOOKUP_TABLE default"
FOR ZZ=0 TO ZSIZE-1
FOR YY=0 TO YSIZE-1
FOR XX=0 TO XSIZE-1
LET X=-2+4/XSIZE*XX
LET Y=-2+4/YSIZE*YY
LET Z=-2+4/ZSIZE*ZZ
CALL VEC3(P,X,Y,Z)
IF HEXPRISM(P,H)<.1 THEN LET DA=255 ELSE LET DA=0
PRINT #1:DA
NEXT XX
NEXT YY
NEXT ZZ
CLOSE #1
END
EXTERNAL FUNCTION HEXPRISM(P(),H()) !'float sdHexPrism( vec3 p, vec2 h ){
DIM Q(3)
CALL VABS(Q,P) !' vec3 q = abs(p);
LET HEXPRISM=MAX(Q(3)-H(2),MAX((Q(1)*SQR(3)/2+Q(2)*.5),Q(2))-H(1)) !' return max(q.z-h.y,max((q.x*0.866025+q.y*0.5),q.y)-h.x);}
END FUNCTION
EXTERNAL SUB VABS(A(),B())
FOR I=1 TO 3
LET A(I)=ABS(B(I))
NEXT I
END SUB
EXTERNAL SUB VEC3(A(),X,Y,Z)
LET A(1)=X
LET A(2)=Y
LET A(3)=Z
END SUB
EXTERNAL SUB VEC2(A(),X,Y)
LET A(1)=X
LET A(2)=Y
END SUB
|
|
|
投稿者:しばっち
投稿日:2018年 3月21日(水)13時01分45秒
|
|
|
> No.4534[元記事へ]
ParaViewのボリュームレンダリング機能を使用しています。
表示を「volume」にしてください。
円柱を定義しています。
マシンパワーによって各格子のサイズを調整してください。
LET XSIZE=100
LET YSIZE=100
LET ZSIZE=100
DIM P(3),H(3)
CALL VEC2(H,.8,.5)
FILE GETSAVENAME F$,"vtkファイル|*.vtk"
IF F$="" THEN STOP
IF POS(UCASE$(F$),".VTK")=0 THEN LET F$=F$&".vtk"
OPEN #1:NAME F$
ERASE #1
PRINT #1:"# vtk DataFile Version 2.0"
PRINT #1:"cylinder"
PRINT #1:"ASCII"
PRINT #1:"DATASET STRUCTURED_POINTS"
PRINT #1:"DIMENSIONS";XSIZE;YSIZE;ZSIZE
PRINT #1:"ORIGIN 0 0 0"
PRINT #1:"SPACING 1 1 1"
PRINT #1:"POINT_DATA";XSIZE*YSIZE*ZSIZE
PRINT #1:"SCALARS scalar float"
PRINT #1:"LOOKUP_TABLE default"
FOR ZZ=0 TO ZSIZE-1
FOR YY=0 TO YSIZE-1
FOR XX=0 TO XSIZE-1
LET X=-1+2/XSIZE*XX
LET Y=-1+2/YSIZE*YY
LET Z=-1+2/ZSIZE*ZZ
CALL VEC3(P,X,Y,Z)
IF CAPPEDCYLINDER(P,H)<.1 THEN LET DA=255 ELSE LET DA=0
PRINT #1:DA
NEXT XX
NEXT YY
NEXT ZZ
CLOSE #1
END
EXTERNAL FUNCTION CAPPEDCYLINDER(P(),H()) !'float sdCappedCylinder( vec3 p, vec2 h ){
DIM V(3),S(3),DD(3),D(3)
LET L=LENGTH2(P(1),P(3),0)
CALL VEC2(V,L,P(2))
CALL VABS(S,V)
MAT D=S-H !'vec2 d = abs(vec2(length(p.xz),p.y)) - h;
CALL NMAX(DD,D,0)
LET CAPPEDCYLINDER=MIN(MAX(D(1),D(2)) ,0)+LENGTH(DD) !' return min(max(d.x,d.y),0.0) + length(max(d,0.0));}
END FUNCTION
EXTERNAL FUNCTION LENGTH(A())
LET LENGTH=SQR(A(1)^2+A(2)^2+A(3)^2)
END FUNCTION
EXTERNAL FUNCTION LENGTH2(A,B,C)
LET LENGTH2=SQR(A^2+B^2+C^2)
END FUNCTION
EXTERNAL SUB NMAX(A(),B(),N)
FOR I=1 TO 3
LET A(I)=MAX(B(I),N)
NEXT I
END SUB
EXTERNAL SUB VABS(A(),B())
FOR I=1 TO 3
LET A(I)=ABS(B(I))
NEXT I
END SUB
EXTERNAL SUB VEC3(A(),X,Y,Z)
LET A(1)=X
LET A(2)=Y
LET A(3)=Z
END SUB
EXTERNAL SUB VEC2(A(),X,Y)
LET A(1)=X
LET A(2)=Y
END SUB
|
|
|
投稿者:しばっち
投稿日:2018年 3月21日(水)13時02分43秒
|
|
|
> No.4534[元記事へ]
ParaViewのボリュームレンダリング機能を使用しています。
表示を「volume」にしてください。
円錐形を定義しています。
マシンパワーによって各格子のサイズを調整してください。
LET XSIZE=100
LET YSIZE=100
LET ZSIZE=100
DIM P(3),C(3)
CALL VEC3(C,.5,.7,.5)
FILE GETSAVENAME F$,"vtkファイル|*.vtk"
IF F$="" THEN STOP
IF POS(UCASE$(F$),".VTK")=0 THEN LET F$=F$&".vtk"
OPEN #1:NAME F$
ERASE #1
PRINT #1:"# vtk DataFile Version 2.0"
PRINT #1:"cone"
PRINT #1:"ASCII"
PRINT #1:"DATASET STRUCTURED_POINTS"
PRINT #1:"DIMENSIONS";XSIZE;YSIZE;ZSIZE
PRINT #1:"ORIGIN 0 0 0"
PRINT #1:"SPACING 1 1 1"
PRINT #1:"POINT_DATA";XSIZE*YSIZE*ZSIZE
PRINT #1:"SCALARS scalar float"
PRINT #1:"LOOKUP_TABLE default"
FOR ZZ=0 TO ZSIZE-1
FOR YY=0 TO YSIZE-1
FOR XX=0 TO XSIZE-1
LET X=-1+2/XSIZE*XX
LET Y=-1+2/YSIZE*YY
LET Z=-1+2/ZSIZE*ZZ
CALL VEC3(P,X,Y,Z)
IF CAPPEDCONE(P,C)<.1 THEN LET DA=255 ELSE LET DA=0
PRINT #1:DA
NEXT XX
NEXT YY
NEXT ZZ
CLOSE #1
END
EXTERNAL FUNCTION CAPPEDCONE(P(),C()) !'float sdCappedCone( in vec3 p, in vec3 c ){
DIM Q(3),V(3),W(3),VV(3),QV(3),D(3),K(3)
CALL VEC2(Q,LENGTH2(P(1),P(3),0),P(2)) !'vec2 q = vec2( length(p.xz), p.y );
CALL VEC2(V,C(3)*C(2)/C(1),-C(3)) !'vec2 v = vec2( c.z*c.y/c.x, -c.z );
MAT W=V-Q !'vec2 w = v - q;
CALL VEC2(VV,DOT(V,V),V(1)^2) !'vec2 vv = vec2( dot(v,v), v.x*v.x );
CALL VEC2(QV,DOT(V,W),V(1)*W(1)) !'vec2 qv = vec2( dot(v,w), v.x*w.x );
CALL NMAX(K,QV,0) !'vec2 d = max(qv,0.0)*qv/vv;
LET D(1)=K(1)*QV(1)/VV(1)
LET D(2)=K(2)*QV(2)/VV(2)
WHEN EXCEPTION IN
LET CAPPEDCONE=SQR(DOT(W,W)-MAX(D(1),D(2)))*SGN(MAX(Q(2)*V(1)-Q(1)*V(2),W(2))) !'return sqrt( dot(w,w) - max(d.x,d.y) ) * sign(max(q.y*v.x-q.x*v.y,w.y));}
USE
LET CAPPEDCONE=100000
END WHEN
END FUNCTION
EXTERNAL FUNCTION LENGTH2(A,B,C)
LET LENGTH2=SQR(A^2+B^2+C^2)
END FUNCTION
EXTERNAL SUB VEC2(A(),X,Y)
LET A(1)=X
LET A(2)=Y
END SUB
EXTERNAL SUB VEC3(A(),X,Y,Z)
LET A(1)=X
LET A(2)=Y
LET A(3)=Z
END SUB
EXTERNAL SUB NMAX(A(),B(),N)
FOR I=1 TO 3
LET A(I)=MAX(B(I),N)
NEXT I
END SUB
|
|
|
投稿者:しばっち
投稿日:2018年 3月21日(水)13時03分23秒
|
|
|
> No.4534[元記事へ]
ParaViewのボリュームレンダリング機能を使用しています。
表示を「volume」にしてください。
カプセル型を定義しています。
マシンパワーによって各格子のサイズを調整してください。
LET XSIZE=100
LET YSIZE=100
LET ZSIZE=100
DIM A(3),B(3),P(3)
CALL VEC3(A,.75,0,1)
CALL VEC3(B,-.75,0,1)
FILE GETSAVENAME F$,"vtkファイル|*.vtk"
IF F$="" THEN STOP
IF POS(UCASE$(F$),".VTK")=0 THEN LET F$=F$&".vtk"
OPEN #1:NAME F$
ERASE #1
PRINT #1:"# vtk DataFile Version 2.0"
PRINT #1:"capsule"
PRINT #1:"ASCII"
PRINT #1:"DATASET STRUCTURED_POINTS"
PRINT #1:"DIMENSIONS";XSIZE;YSIZE;ZSIZE
PRINT #1:"ORIGIN 0 0 0"
PRINT #1:"SPACING 1 1 1"
PRINT #1:"POINT_DATA";XSIZE*YSIZE*ZSIZE
PRINT #1:"SCALARS scalar float"
PRINT #1:"LOOKUP_TABLE default"
FOR ZZ=0 TO ZSIZE-1
FOR YY=0 TO YSIZE-1
FOR XX=0 TO XSIZE-1
LET X=-2.5+5/XSIZE*XX
LET Y=-2.5+5/YSIZE*YY
LET Z=-2.5+5/ZSIZE*ZZ
CALL VEC3(P,X,Y,Z)
IF CAPSULE(P,A,B,1)<.1 THEN LET DA=255 ELSE LET DA=0
PRINT #1:DA
NEXT XX
NEXT YY
NEXT ZZ
CLOSE #1
END
EXTERNAL FUNCTION CAPSULE(P(),A(),B(),R) !'float sdCapsule( vec3 p, vec3 a, vec3 b, float r ){
DIM PA(3),BA(3),D(3)
MAT PA=P-A !'vec3 pa = p - a, ba = b - a;
MAT BA=B-A
LET H=CLAMP(DOT(PA,BA)/DOT(BA,BA),0,1) !'float h = clamp( dot(pa,ba)/dot(ba,ba), 0.0, 1.0 );
MAT D=H*BA
MAT PA=PA-D
LET CAPSULE=LENGTH(PA)-R !'return length( pa - ba*h ) - r;}
END FUNCTION
EXTERNAL FUNCTION LENGTH(A())
LET LENGTH=SQR(A(1)^2+A(2)^2+A(3)^2)
END FUNCTION
EXTERNAL SUB VABS(A(),B())
FOR I=1 TO 3
LET A(I)=ABS(B(I))
NEXT I
END SUB
EXTERNAL SUB VEC3(A(),X,Y,Z)
LET A(1)=X
LET A(2)=Y
LET A(3)=Z
END SUB
EXTERNAL FUNCTION CLAMP(X,A,B) !' A<=X<=B
LET CLAMP=MIN(B,MAX(X,A))
END FUNCTION
|
|
|
戻る