曲面の３次元CG（４）陰関数曲面の描き方

前回までがほぼ基本的な曲面の描き方です。今日は陰関数表現された曲面の描き方についてです。半径５の球は

というパラメータ表示できました。スクリプトでは

SusrTest=({
" 5*cos(u)*cos(v)",
"5*sin(u)*cos(v)",
"5*sin(v)"
});
uvmesh(20,20);
etrgmesh(SusrTest);
end;

となります。ところが良く知られている様に球は次のような式で表す事も出来ます。

こういう表現を陰関数表現と言います。所が陰関数表現による曲面の表示はとても難しい事が知られています。簡単に言えば方程式の解（x,y,z）の集合だからです。さて面倒な陰関数表現による曲面の表示は次のように簡単に行う事が出来ます。

expr = "x^2 + y^2 + z^2 - 5";
implicitsurface( getstr(expr), [-3,-3,-3],[3,3,3],[0.5, 0.5, 0.5], 3);
end;

exprに数式を記述しています。[-3,-3,-3],[3,3,3]は曲面が収まるボックスの頂点座標です。
格子間隔[0.5, 0.5, 0.5]が必要です。この間隔で作成した３次元の区画で方程式の値value

を求めて行きます。このvalueに応じて各格子の頂点に0～255の値が設定されます。ちょうどグレースケールの画像のようになっています。最後の3は方程式の解である境界を0～255の範囲で指定します。
さらにもう一つの描き方としては以下のようにimplicitsurface2を使う方法があります。

expr = "x^2 + y^2 + z^2 - 5";
implicitsurface2( getstr(expr), [-3,-3,-3],[3,3,3],[0.0, 0.0, 0.0]);
end;

exprに数式を記述し、[-3,-3,-3],[3,3,3]は曲面が収まるボックスの頂点座標なのは同じです。
それと初期点[0.5, 0.5, 0.5]が必要です。初期点はこの近傍に方程式の解が在る事を知らせるためのヒントのようなものです。

expr= = "(1/(x^2/9+4*y^2+4*z^2)^4+1/(y^2/9+4*x^2+4*z^2)^4+1/(z^2/9+4*y^2+4*x^2)^4+
1/((4*x/3-4)^2+16*y^2/9+16*z^2/9)^4+1/((4*x/3+4)^2+16*y^2/9+16*z^2/9)^4+
1/((4*y/3-4)^2+16*x^2/9+16*z^2/9)^4+1/((4*y/3+4)^2+16*x^2/9+16*z^2/9)^4)^(-1/4)-1";

implicitsurface( getstr(expr), [-6,-6,-6],[6,6,6],[0.15,0.15,0.15], 1, 0.001);
end;

#%
sub func(n, xyz)
{
x = xyz[0];
y = xyz[1];
z = xyz[2];

hx =(x-1)*x*x*(x+1);
f = pow((hx+y*y),2)+z*z;
func.ret = f;
return (f);
}


put 1;
implicitsurface2("u[0]=x;u[1]=y;u[2]=z;t_=$func(3,u);t_-0.02",[-0.7,-0.73,-0.3],[0.7,0.73,0.3],[1,1,1],0.025);
implicitsurface2("u[0]=x;u[1]=y;u[2]=z;t_=$func(3,u);t_-0.02",[-1.5,-0.73,-0.3],[0.7,0.73,0.3],[1,1,1],0.025);
implicitsurface2("u[0]=x;u[1]=y;u[2]=z;t_=$func(3,u[]);t_-0.02",[0.6,-0.73,-0.3],[1.5,0.73,0.3],[1.5,0,0],0.025);
end;

今日はこんなところで。