はじめに
intersection法でトーラスを描いたことがある。目線の位置からレイをとばし、トーラスと交叉する場合にその点における法線とかを計算して、色を付ける。過去にそういうのを書いたんだが、
OpenProcessingの仕様変更で壊れてしまったので、復活させた。
方針
通常の3D描画と同じように、eye,center,topを用意する。ただし射影が必要なポリゴン描画と違い、ダイレクトに色を出せるので、道のりはだいぶ異なる。それでもビュー座標系の算出までは一緒である。詳しい解説は省略。なお透視投影です。
コード全文
uvを原点中心、上下が-1～1,左右が-1～1になっているようなuv空間で2D描画をする独自のレンダー機構でやってますがshadertoyと同じような感じなのでなんとなくわかるでしょう。
// これのサルベージ
// https://openprocessing.org/@u161812/882125
// もう動かないので。
// 四次方程式が解きたい人
// https://qiita.com/TETSURO1999/items/e0d8ab5055ab56a0180b
// 自分の記事
// https://qiita.com/inaba_darkfox/items/22a330ba4b3ae8aac7b7
const {Render2D} = fisce.foxApplications;
function setup() {
createCanvas(400,400,WEBGL);
smooth();
const gl = drawingContext;
// cam: [0,2,4] [0,0,0] [0,1,0] Math.PI/3 1 0.04 100;
// ray marching/intersectionの場合、near/farに意味は無いです。fovとaspectから画角などが決まります。
// やることは単純で対象のトーラスを...
// 対象のトーラスは？
// 原点中心、yUpで、長半径1で、短半径は0.3ですね。半径というのはそのまんまの意味。
// eyeから-frontに1だけ離れたところにキャンバスを置く。
// それの上端はtan(fov/2)*upで出る
// それの右端はtan(fov/2)*aspect*sideで出る
// ゆえに
// dir = -front + tan(fov/2)*uv.x*up + tan(fov/2)*aspect*uv.y*side となるわけ。
const r2d = Render2D.create(gl);
const src = `
@fs
<global>
# 4次方程式を解くっぽい。
// 4x^3 - 3qx = rの解を求める。
float getSub(float q, float r){
// rが0とみなせるならふつうに0.0やらqの平方根でやる。
if(abs(r) < 1e-10){
if(q < 0.0){ return 0.0; }
return 0.5 * sqrt(3.0 * q);
}
// 以降、|r|>0とする。符号をとっておく。
float sign_r = sign(r);
// qが0とみなせるなら3乗根を取ればよい
if(abs(q) < 1e-10){
return sign_r * pow(0.25 * abs(r), 1.0 / 3.0);
}
// 以降は|r|>0かつ|q|>0とする。
float d = r * r - q * q * q; // 判別式
// この辺りの処理ではxの係数を変数変換で1にしている
float cf = sqrt(abs(q));
float h = abs(r) / pow(cf, 3.0);
// まず判別式が負の場合は三角関数解が得られる
if(d < 0.0){
float alpha = acos(sign_r * h) / 3.0;
return cf * cos(alpha);
}
// 判別式が0以上の場合は実数解はひとつなのでそれを
// 双曲線関数を用いて取得する
// qが正の時はコサインハイポ
if(q > 0.0){
float x = log(h + sqrt(h * h - 1.0)) / 3.0;
// rが負の時はマイナスを付ける
return sign_r * cf * cosh(x);
}
// qが負の時はサインハイポ
float x = log(h + sqrt(h * h + 1.0)) / 3.0;
return sign_r * cf * sinh(x);
}
// 4次方程式の実数解を取得するパート。
// x^4 + 4(k3)x^3 + 4(k2)x^2 + 8(k1)x + 4(k0) = 0の解を求める
vec4 solve4(float k3, float k2, float k1, float k0){
// あらかじめ-1.0で埋めておいて解が見つかったら置き換える感じ。
// どうせ正の数しか使わないでしょ
vec4 ans = vec4(-1.0);
// 変数変換してx^3の係数をなくす処理
// 具体的には解の和が0になるようにグラフの平行移動をしている
float c2 = (2.0 * k2 - 3.0 * k3 * k3) / 3.0;
float c1 = 2.0 * (k3 * k3 * k3 - k2 * k3 + k1);
float c0 = (-3.0 * pow(k3, 4.0) + 4.0 * k3 * k3 * k2 - 8.0 * k1 * k3 + 4.0 * k0) / 3.0;
// c1が0とみなせるなら退化として計算
if(abs(c1) < 1e-10){
if(3.0 * c2 * c2 < c0){ return ans; }
float beta = sqrt(9.0 * c2 * c2 - 3.0 * c0);
float alpha = -3.0 * c2;
if(alpha + beta >= 0.0){
ans.x = sqrt(alpha + beta) - k3;
ans.y = -sqrt(alpha + beta) - k3;
}
if(alpha - beta >= 0.0){
ans.z = sqrt(alpha - beta) - k3;
ans.w = -sqrt(alpha - beta) - k3;
}
return ans;
}
// 以下では|c1| > 0.0とする。続きは、帰ってから（えー）
float q = c0 + c2 * c2;
float r = c1 * c1 + c2 * c2 * c2 - 3.0 * c0 * c2;
// q, rに対して4x^3 - qx = rの解のうち最大の実数値を取る
float w = getSub(q, r);
// するとv = w - c2が正の数になる
float v = w - c2;
// vと係数から解を出す
float j = sqrt(v);
float h = c1 / j;
// 判別式
float d1 = -v - 3.0 * c2 - h;
float d2 = -v - 3.0 * c2 + h;
if(d1 >= 0.0){
d1 = sqrt(d1);
ans.x = j + d1 - k3;
ans.y = j - d1 - k3;
}
if(d2 >= 0.0){
d2 = sqrt(d2);
ans.z = -j + d2 - k3;
ans.w = -j - d2 - k3;
}
return ans;
}
// トーラス。
// cは中心、nは法線ベクトル、aは軸半径、bは胴体半径。
float getTorus(vec3 ori, vec3 dir, vec3 c, vec3 n, float a, float b){
// まずk3, k2, k1, k0を出す。
float q0 = dot(ori - c, ori - c);
float q1 = dot(ori - c, dir);
float q2 = pow(dot(ori - c, n), 2.0);
float q3 = dot(ori - c, n) * dot(dir, n);
float q4 = pow(dot(dir, n), 2.0);
float k = 0.5 * (q0 - a * a - b * b);
float k3 = q1;
float k2 = k + q1 * q1 + a * a * q4;
float k1 = q1 * k + a * a * q3;
float k0 = k * k - a * a * b * b + a * a * q2;
vec4 ans = solve4(k3, k2, k1, k0);
if(ans.x < 0.0 && ans.y < 0.0 && ans.z < 0.0 && ans.w < 0.0){ return -1.0; }
float t = 1e20;
if(ans.x >= 0.0){ t = min(t, ans.x); }
if(ans.y >= 0.0){ t = min(t, ans.y); }
if(ans.z >= 0.0){ t = min(t, ans.z); }
if(ans.w >= 0.0){ t = min(t, ans.w); }
return t;
}
vec3 getNormalOfTorus(vec3 p, vec3 c, vec3 n, float a, float b){
vec3 q = c + normalize(p - c - dot(p - c, n) * n) * a;
return p-q;
}
<main>
vec3 eye = vec3(0.0, 2.0, 3.0);
vec3 center = vec3(0.0);
vec3 top = vec3(0.0, 1.0, 0.0);
vec3 front = normalize(eye-center);
vec3 side = normalize(cross(top, front));
vec3 up = cross(front, side);
// uvは-1～1の-1～1っすね。それで...
float fov = 1.0;
float aspect = 1.0;
vec3 dir = normalize(-front + tan(fov*0.5)*uv.y*up + tan(fov*0.5)*aspect*uv.x*side);
vec3 torusCenter = vec3(0.0);
vec3 torusNormal = vec3(0.0, 1.0, 0.0);
float a = 1.0;
float b = 0.3;
float t = getTorus(eye, dir, torusCenter, torusNormal, a, b);
if(t < 0.0){
color = vec4(0.0, 0.0, 0.0, 1.0);
}else{
vec3 p = eye + t*dir;
vec3 normal = getNormalOfTorus(p, torusCenter, torusNormal, a, b);
if(length(normal) > b+0.05){
color = vec4(vec3(0.0), 1.0);
}else{
normal = normalize(normal);
float blt0 = max(0.2, dot(normal, normalize(vec3(0.0, 1.0, 0.0))));
vec3 col0 = blt0*vec3(0.5, 1.0, 1.5);
color = vec4(col0, 1.0);}
}
`;
r2d.build({
shader:{align:'center_yUp'},
layout:src
});
draw = () => {
gl.clearColor(0,0,0,1);
gl.clear(gl.COLOR_BUFFER_BIT);
r2d.useProgram();
r2d.render();
r2d.flush();
}
}
結果：
最後の方、表面上の点であることを仮定して法線を取るところでちょっとインチキしていますが、まあいいでしょう。できました。当時の自分が見ていた景色がちょっとだけ追体験出来ました。これで満足するとしましょう。まあこういうのは通常レイマーチングやポリゴン描画でやるんですよね...
ああー...今見返すとなんか黒い点が見えてますね。たぶんこのトーラスは太陽かなんかで、黒点ですね。青く見えますがきっと暑いんだと思います。
ぽこポケやりたい。
おわりに
ここまでお読みいただいてありがとうございました。
自分の記事がリンクしてもらえてうれしかったので、ちょっと過去を懐かしんでみたくなったのでした。
追記
Waterfoxだから動かなかったようで、Chrome系だと動きました。
だからなに？
ブラウザによって挙動が不安定になるようなコードなど、書きたくありません。今回書いたコードはどちらでも問題なく動きます。それが一番です。私が今愛用してるのはWaterfoxなので、それで動かないならそれは駄目ということです。