|
> 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
|
|