LoginSignup
TrendQuestionOfficial EventOfficial ColumnOpportunitiesOrganization
3
1

More than 3 years have passed since last update.

@sage-git

posted at

updated at

３次元スカラー場をレイマーチングして等高面描写

JavaScript,GLSL

レイマーチングでは距離関数から光線の進む距離を決めるのが普通ですが、光線の進む距離さえ何らかの手段で決められたら距離関数がなくとも同様のことができるのでは？と試してみました。

結果

三次元中のスカラー場について、目的とする値を取る地点を集めた等高面の描写に成功しました。

p_d_solid.gif

鏡面反射するd軌道とか初めて見た。

導入

３次元スカラー場と言っているのは、要は３次元中で定義された関数 $f(x, y, z)$ です。やりたいことはこの関数が $c$を取るような点$(x, y, z)$を集めた平面$\Omega$を描くことです。

\Omega = \{(x, y, z) \mid f(x, y, z) = c\}

例えば、３次元空間中の波動関数の等高面を描写することを考えています。実際結果の動画はp軌道とd軌道です。

アイデア

光線が $i$ ステップ目に 

\vec{r_i} = \left(\begin{matrix}x_i\\ y_i\\ z_i\end{matrix}\right)

にいた時、

f_i = f(\vec{r_i})\\
\vec{d_i} = \nabla f(\vec{r_i})

を計算します。

そして次に進むべき距離の決定ですが、普通のレイマーチングなら距離関数から素直に計算できます。というかそれが距離関数の定義です。しかし今回はそれができません。

一次元だと、このような絵を描きます。

test.png

図の$L$のように、現在位置での値と目的の値の差を今いるところの傾きで割った量は、進むべき距離の参考になるはずだというのがアイデアの根本です。

この図の$L$を更に具体的に決めていきます。とりあえず以下のように$\beta$を足します。

L = p \frac{ c - f_i }{\beta + \left|\vec{d_i}\right|}

ここで、$\beta$ は $\left|\vec{d_i}\right|$ が $0$ だったときに発散しないようにするためです。

なんか図だと $x$ での傾きを延長させたら $f(x) = c$ となる $x$ を遥か通りすぎている気がしますが、実際これは起こりうる可能性がありますので、 $p$ や $\beta$ で行き過ぎないように抑えます。

とりあえずこれで実装してみて、ダメだったら試行錯誤していくスタイル。

法線

$\Omega$ 上での $(x, y, z)$ の法線 $\vec{n}$ は実は難しいことはなくて、

\vec{n} = \nabla f

で計算できます。
つまり、$\nabla f$の計算を行う関数は、光線をすすめる距離を計算するだけでなく、法線ベクトルを得るのにも用います。

なので、今回定義しなければならないのは $f(\vec{r})$ と $ \nabla f(\vec{r})$ の具体的な表式と、
これらの実際に計算した値と$c$から進む距離を決定するアルゴリズムです。

暫定のパラメータ

p = 1.0\\
\beta = 0.5

他に、波動関数中でのパラメータとして $\alpha = 2.0$としています。

暫定の関数

p軌道

f(\vec{r}) = \frac{x}{|\vec{r}|}\exp\left(-\frac{\vec{r}^2}{\alpha^2}\right)\\
\nabla f(\vec{r}) = \left(\begin{matrix}
\vec{r}^2 - x^2\\ -xy\\ -xz\end{matrix}\right)
\exp\left(-\frac{\vec{r}^2}{\alpha^2}\right) - \frac{2\vec{r}f(\vec{r})}{\alpha^2}

d軌道

f(\vec{r}) = \frac{2z^2 - x^2 - y^2}{\vec{r}^2}\exp\left(-\frac{\vec{r}^2}{\alpha^2}\right)\\
\nabla f(\vec{r}) = \left(\begin{matrix}-2x\\ -2y\\ 4z\end{matrix}\right)\exp\left(-\frac{\vec{r}^2}{\alpha^2}\right) - 2\vec{r}f(\vec{r})\left(\frac{1}{\vec{
r}^2} + \frac{1}{\alpha^2} \right)

ソースコード

  • 2020/3/12 追記

Githubにまとめました。

非日本語圏兄貴がそれなりに見られますので、README.mdはクソザコイングリッシュです。

改めて、必要なファイルを並べます。

https://qiita.com/sage-git/items/662ec5aae1995c51f035
にもあるように、
https://github.com/gam0022/webgl-sandbox/tree/master/js
を用意します。

style.css(クリックして展開)
style.css
body {
    background-color: black;
    margin: 0;
    padding: 0;
}
a { color: skyblue }
canvas {
    display: block;
    position: absolute;
    top: 0;
    left: 0;
    right: 0;
    bottom: 0;
    margin: auto;
}
#info {
    color: white;
    font-size: 13px;
    position: absolute;
    bottom: 10px;
    width: 100%;
    text-align: center;
    z-index: 100;
}
#fragment-code {
    overflow: scroll;
    white-space: nowrap;
    font-family: 'Courier New', Courier, monospace;
}

以下のfragファイル群をfrag_filesフォルダに置きます。

headerとrendering２種類
header.frag
precision mediump float;
uniform vec2 resolution;
uniform mat4 viewMatrix;
uniform vec3 cameraPosition;
uniform mat4 cameraWorldMatrix;
uniform mat4 cameraProjectionMatrixInverse;
const float EPS = 0.01;
const float OFFSET = EPS * 5.0;
//const vec3 lightDir = vec3( -0.48666426339228763, 0.8111071056538127, -0.3244428422615251 );
const vec3 lightDir = vec3( 0.0, 1.0, 0.0 );

const int ITER = 80;
const float ITER_SHADOW = 128.0;

const float ISOVAL = 0.2;
const float step = 0.5;
float wavefunc( vec3 p );
vec3 wavefunc_d( vec3 p, float f );
render_wavefunc_solid.frag
float sceneDist( vec3 p, float f, vec3 df, vec3 ray ) {
    float nearest = length(p);
    float dlen = length(df);
    float angle = dot(df, ray);
    float dr = (ISOVAL - abs(f))/(0.5 + abs(angle));
    return dr;
}

vec3 sceneColor( vec3 p, float f, vec3 df ) {
    if( abs(abs(f) - ISOVAL) < EPS*2.0){
        return (f > 0.0)? vec3(1.0, 0.0, 0.0): vec3(0.0, 0.0, 1.0);
   }
   return vec3(0.0, 0.0, 0.0);
}

vec3 getNormal( vec3 p, float f, vec3 df ) {
    float len = length(df);
    return - sign(f)*df/len;
}

float getShadow( vec3 ro, vec3 rd, float f, vec3 df ) {
    float h = 0.0;
    float c = 0.0;
    float r = 1.0;
    float shadowCoef = 0.5;
    for ( float t = 0.0; t< ITER_SHADOW; t++ ) {
        h = sceneDist( ro + rd * c, f, df, rd );
        if ( h < EPS ) return shadowCoef;
        r = min( r, h * 16.0 / (c + EPS));
        c += h;
    }
    return 1.0 - shadowCoef + r * shadowCoef;
}

vec3 getRayColor( vec3 origin, vec3 ray, out vec3 pos, out vec3 normal, out bool hit ) {
    // marching loop
    float dist;
    float depth = 0.0;
    pos = origin;
    for ( int i = 0; i < ITER; i++ ){
        float f = wavefunc(pos);
        vec3 df = wavefunc_d(pos, f);
        dist = sceneDist( pos, f, df, ray );
        depth += dist;
        pos = origin + depth * ray;
        if ( abs(dist) < EPS ) break;
    }
    // hit check and calc color
    vec3 color;
    const float matness = 1.0;

    if ( abs(dist) < EPS ) {
        float f = wavefunc(pos);
        vec3 df = wavefunc_d(pos, f);
        normal = getNormal( pos, f, df );
        float diffuse = clamp( dot( lightDir, normal ), 0.1, 1.0 );
        float specular = clamp( dot( reflect( lightDir, normal ), ray ), 0.0, 1.0 );
        //float shadow = getShadow( pos + normal * OFFSET, lightDir, f, df );
        //color = ( sceneColor( pos, f, df ) + vec3( 1.0 - matness ) * specular + vec3( matness) * diffuse ) * max( 0.5, shadow );
        color = ( sceneColor( pos, f, df ) + vec3( 1.0 - matness ) * specular + vec3( matness) * diffuse )*0.5;
        hit = true;
    } else {
        color = vec3( 0.4 );
    }
    return color;// - pow( clamp( 0.05 * depth, 0.0, 0.6 ), 2.0 ) * 0.1;
}
void main(void) {
    // screen position
    vec2 screenPos = ( gl_FragCoord.xy * 2.0 - resolution ) / min( resolution.x, resolution.y );
    // convert ray direction from screen coordinate to world coordinate
    vec3 ray = (cameraWorldMatrix * cameraProjectionMatrixInverse * vec4( screenPos.xy, 1.0, 1.0 )).xyz;
    ray = normalize( ray );
    // camera position
    vec3 cPos = cameraPosition;
    // cast ray
    vec3 color = vec3( 0.0 );
    vec3 pos, normal;
    bool hit;
    float alpha = 1.0;
    for ( int i = 0; i < 3; i++ ) {
        color += alpha * getRayColor( cPos, ray, pos, normal, hit );
        alpha *= 0.3;
        ray = normalize( reflect( ray, normal ) );
        cPos = pos + normal * OFFSET;
        if ( !hit ) break;
    }
    gl_FragColor = vec4( color, 1.0 );
}
render_wavefunc_fog.frag
const float FAREST = 30.0;
float sceneDist( vec3 p, float f, vec3 df, vec3 ray, out vec3 fog_color ) {
    float nearest = length(p);
    float dlen = length(df);
    fog_color = (f > 0.0)? vec3(1.0, 0.0, 0.0): vec3(0.0, 0.0, 1.0);
    float dr = clamp(step/(1.0 + dlen), EPS*2.0, step*2.0);
    fog_color = dr*fog_color*(f*f)*2.0;
    return dr;
}
vec3 sceneColor( vec3 p, float f, vec3 df ) {
    return vec3(0.3, 0.8, 0.9);
}
vec3 getNormal( vec3 p, float f, vec3 df ) {
    return p/length(p);
}

float getShadow( vec3 ro, vec3 rd, float f, vec3 df, out vec3 fog ) {
    float h = 0.0;
    float c = 0.0;
    float r = 1.0;
    float shadowCoef = 0.5;
    fog = vec3(0.0);
    for ( float t = 0.0; t< ITER_SHADOW; t++ ) {
        vec3 fog_i = vec3(0.0);
        h = sceneDist( ro + rd * c, f, df, rd, fog_i );
        fog += fog_i;
        if ( h < EPS ) return shadowCoef;
        r = min( r, h * 16.0 / (c + EPS));
        c += h;
    }
    return 1.0 - shadowCoef + r * shadowCoef;
}
vec3 getRayColor( vec3 origin, vec3 ray, out vec3 pos, out vec3 normal, out bool hit ) {
    // marching loop
    float dist = 1.0;
    float depth = 0.0;
    pos = origin;
    vec3 fog = vec3(0.0);
    float f = wavefunc(pos);
    for ( int i = 0; i < ITER; i++ ){
        f = wavefunc(pos);
        vec3 df = wavefunc_d(pos, f);
        vec3 fog_i = vec3(0.0);
        dist = sceneDist( pos, f, df, ray, fog_i );
        fog += fog_i;
        depth += dist;
        pos = origin + depth * ray;
        if ( abs(dist) < EPS || abs(depth) > FAREST) break;
    }
    // hit check and calc color
    vec3 color;
    const float matness = 1.0;

    if ( abs(dist) < EPS ) {
        float f = wavefunc(pos);
        vec3 df = wavefunc_d(pos, f);
        normal = getNormal( pos, f, df );
        float diffuse = clamp( dot( lightDir, normal ), 0.1, 1.0 );
        float specular = clamp( dot( reflect( lightDir, normal ), ray ), 0.0, 1.0 );
        //float shadow = getShadow( pos + normal * OFFSET, lightDir, f, df );
        //color = ( sceneColor( pos, f, df ) + vec3( 1.0 - matness ) * specular + vec3( matness) * diffuse ) * max( 0.5, shadow );
        color = ( sceneColor( pos, f, df ) + vec3( 1.0 - matness ) * specular + vec3( matness) * diffuse )*0.5;
        hit = true;
    } else {
        color = vec3( 0.4 );
    }
    return clamp(fog + color, 0.0, 1.0);// - pow( clamp( 0.05 * depth, 0.0, 0.6 ), 2.0 ) * 0.1;
}
void main(void) {
    // screen position
    vec2 screenPos = ( gl_FragCoord.xy * 2.0 - resolution ) / min( resolution.x, resolution.y );
    // convert ray direction from screen coordinate to world coordinate
    vec3 ray = (cameraWorldMatrix * cameraProjectionMatrixInverse * vec4( screenPos.xy, 1.0, 1.0 )).xyz;
    ray = normalize( ray );
    // camera position
    vec3 cPos = cameraPosition;
    // cast ray
    vec3 color = vec3( 0.0 );
    vec3 pos, normal;
    bool hit;
    float alpha = 1.0;
    for ( int i = 0; i < 3; i++ ) {
        color += alpha * getRayColor( cPos, ray, pos, normal, hit );
        alpha *= 0.3;
        ray = normalize( reflect( ray, normal ) );
        cPos = pos + normal * OFFSET;
        if ( !hit ) break;
    }
    gl_FragColor = vec4( color, 1.0 );
}

波動関数
wave_px.frag
const float ALPHA = 2.0;
float wavefunc( vec3 p ){
    float r = length(p);
    return exp( - r*r/ALPHA/ALPHA)*p.x/r;
}
vec3 wavefunc_d( vec3 p, float f ){
    float r = length(p);
    if ( r == 0.0 ){ return vec3(1.0, 0.0, 0.0); }
    float pexp = exp( - r*r/ALPHA/ALPHA);
    vec3 ret = vec3(r*r - p.x*p.x, - p.x*p.y, -p.x*p.z)/(r*r*r)*pexp - f*2.0*r*p/(ALPHA*ALPHA*r);
    return ret;
}
wave_dzr.frag
const float ALPHA = 2.0;
float wavefunc( vec3 p ){
    float r = length(p);
    return exp( - r*r/ALPHA/ALPHA)*(2.0*p.z*p.z - p.x*p.x - p.y*p.y)/r/r;
}
vec3 wavefunc_d( vec3 p, float f ){
    float r = length(p);
    if ( r == 0.0 ){ return vec3(1.0, 0.0, 0.0); }
    float pexp = exp( - r*r/ALPHA/ALPHA);
    vec3 ret = vec3(- 2.0*p.x, -2.0*p.y, 4.0*p.z)*pexp/r/r - 2.0*f*p/r/r - 2.0*f*p/ALPHA/ALPHA;
    return ret;
}

最後に本体のHTMLファイルです。

index.html
<!DOCTYPE html>
<html lang="en">
    <head>
        <title>three.js webgl - raymarching - reflect</title>
        <meta charset="utf-8">
        <meta name="viewport" content="width=device-width, user-scalable=no, minimum-scale=1.0, maximum-scale=1.0">
        <link rel="stylesheet" href="style.css" type="text/css">
    </head>
    <body>
<script id="vertex_shader" type="x-shader/x-vertex">
    attribute vec3 position;
    void main(void) {
        gl_Position = vec4(position, 1.0);
    }
</script>

<script src="js/three.min.js"></script>
<script src="js/controls/FlyControls.js"></script>
<script src="js/libs/stats.min.js"></script>
<script src="js/libs/dat.gui.min.js"></script>

<script>
var camera, scene, controls, renderer;
var geometry, material, mesh;
var mouse = new THREE.Vector2( 0.5, 0.5 );
var canvas;
var view_center;
var stats;
var rotation_function = function(x, y){};
var clock = new THREE.Clock();
var config = {
    saveImage: function() {
        renderer.render( scene, camera );
        window.open( canvas.toDataURL() );
    },
    resolution: '512',
    cameraReset: camera_init,
    renderPattern: "p wave solid",
    atom_on: "ON",
    load_code: function() {
        fragment_code = document.getElementById("fragment-code");
        loadResources([], fragment_code);
    },
};

var pwave_set = ['header.frag', 'wave_px.frag', 'render_wavefunc_solid.frag'];
var dwave_set = ['header.frag', 'wave_dzr.frag', 'render_wavefunc_fog.frag'];
var dwave_solid_set = ['header.frag', 'wave_dzr.frag', 'render_wavefunc_solid.frag'];

var fragment_set ={
    "p wave solid": ['header.frag', 'wave_px.frag', 'render_wavefunc_solid.frag'],
    "d wave fog": ['header.frag', 'wave_dzr.frag', 'render_wavefunc_fog.frag'],
    "d wave solid": ['header.frag', 'wave_dzr.frag', 'render_wavefunc_solid.frag'],
}

async function loadResources(urls, resource_store_textarea) {
    if(urls.length == 0){
        material_init(resource_store_textarea.value);
        render();
        return;
    }
    var promises = urls.map(function(url){
        return new Promise( function (resolve, reject){
            var request = new XMLHttpRequest();
            function error_callback() {
                console.log("Error while loading '" + url + "'.");
                reject(new Error(request.statusText));
            }
            request.open('GET', "frag_files/" + url, true);
            request.setRequestHeader('Pragma', 'no-cache');
            request.setRequestHeader('Cache-Control', 'no-cache');
            request.setRequestHeader('If-Modified-Since', 'Thu, 01 Jun 1970 00:00:00');
            request.responseType = "text";
            request.onload = function () {
                if (request.readyState == 4 || request.status == 200) {
                    console.log("Successfully '" + url + "' loaded.");
                    resolve(request.responseText);
                } else { 
                    error_callback();
                }
            };
            request.onerror = error_callback;
            request.send(null);
        });

    });

    result = (await Promise.all(promises)).join("\n");
    resource_store_textarea.value = result;
    compile_code(result);
}

function camera_init(){
    camera.position.set(0.0, 0.0, 5.0);
    camera.lookAt( new THREE.Vector3( 0.0, 0.0, 0.0 ));
    camera.up.set(0.0, 1.0, 0.0);
    view_center = new THREE.Vector3(0.0, 0.0, 0.0);
}

var rotation_horizon = function(angle){
    camera.position.sub(view_center);
    var r = camera.position.length();
    right_dir = camera.position.clone().cross(camera.up).normalize();
    camera.position.addScaledVector(right_dir, r*angle/2.0);
    camera.position.setLength(r);
    camera.position.add(view_center);
}
function rotation_vertical(angle){
    camera.position.sub(view_center);
    var r = camera.position.length();
    camera.position.addScaledVector(camera.up, r*angle/2.0);
    camera.position.setLength(r);
    camera.up = camera.position.clone().cross(camera.up).normalize();
    camera.up.cross(camera.position).normalize();
    camera.position.add(view_center);
}

function move_forward(dx){
    dr = camera.position.clone().sub(view_center).multiplyScalar(- dx);
    camera.position.add(dr);
    view_center.add(dr)
}

function camera_pan(angle){
    look_vec = view_center.clone().sub(camera.position);
    r = look_vec.length();
    right_dir = look_vec.clone().cross(camera.up);
    look_vec.addScaledVector(right_dir, angle/2.0).setLength(r);
    view_center.addVectors(camera.position, look_vec);
}

function camera_tilt(angle){
    look_vec = view_center.clone().sub(camera.position);
    r = look_vec.length();
    look_vec.addScaledVector(camera.up, angle*2.0).setLength(r);
    view_center.addVectors(camera.position, look_vec);
}

function material_init(fragment_code){
    scene = new THREE.Scene();
    geometry = new THREE.PlaneBufferGeometry( 2.0, 2.0 );
    material = new THREE.RawShaderMaterial( {
        uniforms: {
            resolution: { value: new THREE.Vector2( canvas.width, canvas.height ) },
            cameraWorldMatrix: { value: camera.matrixWorld },
            cameraProjectionMatrixInverse: { value: new THREE.Matrix4().getInverse( camera.projectionMatrix ) },
        },
        vertexShader: document.getElementById( 'vertex_shader' ).textContent,
        fragmentShader: fragment_code
    } );
    mesh = new THREE.Mesh( geometry, material );
    mesh.frustumCulled = false;
    scene.add( mesh );
}

function init() {
    renderer = new THREE.WebGLRenderer();
    renderer.setPixelRatio( window.devicePixelRatio );
    renderer.setSize( config.resolution, config.resolution );
    canvas = renderer.domElement;
    canvas.addEventListener( 'mousedown', onMouseDown, false );
    canvas.addEventListener( 'mouseup', onMouseUp, false );
    canvas.addEventListener( 'mouseout', function(e){onMouseUp(e);}, false );
    canvas.addEventListener( 'mousemove', onMouseMove, false );
    canvas.addEventListener( 'mousewheel', onMouseWheel, false);
    document.addEventListener( 'keydown', onKeyDown, false);
    window.addEventListener( 'resize', onWindowResize );
    document.body.appendChild( canvas );
    // Scene
    camera = new THREE.PerspectiveCamera( 60, canvas.width / canvas.height, 1, 2000 );
}

function init_controls(){
    // Controls
    controls = new THREE.FlyControls( camera, canvas );
    controls.autoForward = true;
    controls.dragToLook = false;
    controls.rollSpeed = Math.PI / 12;
    controls.movementSpeed = 0.5;
    // GUI
    var gui = new dat.GUI();
    gui.add( config, 'saveImage' ).name( 'Save Image' );
    gui.add( config, 'resolution', [ '256', '512', '800', 'full' ] ).name( 'Resolution' ).onChange( function( value ) {
        if ( value !== 'full' ) {
            canvas.width = value;
            canvas.height = value;
        }
        onWindowResize();
    } );
    gui.add( config, 'cameraReset').name("Reset Camera");
    gui.add( config, 'renderPattern', Object.keys(fragment_set)).name("Pattern").onChange( function( value ){
        fragment_files = fragment_set[value];
        fragment_code = document.getElementById("fragment-code");
        loadResources(fragment_files, fragment_code);
    } );
    gui.add( config, 'load_code').name("Load Code");
    stats = new Stats();
    document.body.appendChild( stats.domElement );
    fragment_code = document.createElement("textarea");
    fragment_code.id = "fragment-code";
    fragment_code.onChange = function(e){load_render_code();};
    document.body.appendChild(fragment_code);
}

function compile_code(){
    fragment_code = document.getElementById("fragment-code");
    material_init(fragment_code.value);
    render();
}

function render( timestamp ) {
    var delta = clock.getDelta();
    stats.begin();
    camera.lookAt(view_center);
    material.uniforms.resolution.value.set( canvas.width, canvas.height );
    material.uniforms.cameraProjectionMatrixInverse.value.getInverse( camera.projectionMatrix );
    renderer.render( scene, camera );
    stats.end();
    requestAnimationFrame( render );
}

function onMouseDown( event ) {
    rotation_function = function(x, y){
        rot_speed = 10.0;// controls.rollSpeed;
        rotation_horizon(x*rot_speed);
        rotation_vertical(y*rot_speed);
    }
}

function onMouseMove( event ) {
    px = mouse.x;
    py = mouse.y;
    mouse.x = event.clientX / canvas.width;
    mouse.y = event.clientY / canvas.height;
    rotation_function(mouse.x - px, mouse.y - py);
}

function onMouseUp( event ) {
    rotation_function = function(x, y){};
}

var delta = 0.05;
function onMouseWheel( event ){
    if(event.wheelDelta > 0){
        camera.position.multiplyScalar(1. + delta);
    }else{
        camera.position.multiplyScalar(1. - delta);
    }
}

function onKeyDown( event ){
    if(event.keyCode === 87){
        move_forward(0.01);
    }
    if(event.keyCode === 83){
        move_forward(-0.01);
    }
    if(event.keyCode === 37){
        camera_pan(0.05);
    }
    if(event.keyCode === 39){
        camera_pan(-0.05);
    }
    if(event.keyCode === 38){
        camera_tilt(0.05);
    }
    if(event.keyCode === 40){
        camera_tilt(-0.05);
    }
}

function onWindowResize( e ) {
    if ( config.resolution === 'full' ) {
        canvas.width = window.innerWidth;
        canvas.height = window.innerHeight;
    }
    renderer.setSize( canvas.width, canvas.height );
}

init();
init_controls();
camera_init();
loadResources(pwave_set, fragment_code);
        </script>
    </body>
</html>

申し訳ありませんが、コードの解説は力尽きたのでパスさせてください。

今後の工夫点

注目する点から遠くの点にいた時、スカラー場はほとんど変化がなくなると思われます。このとき、$\vec{d_i}$がほぼ$0$となります。すると、分母の補正$\beta$が小さすぎると次に進む距離$L$が大きくなりすぎ、等高面を描くべき地点を通りすぎてしまう可能性があります。しかし、$\beta$ を大きくする、つまり慎重に進むようにするとイテレーション回数が無駄に嵩み、最悪何も見えない状況になりえます。

そして$f(\vec{r})$の値のオーダーも重要です。単に定数倍した関数$g = 100 \times f$を考えてみると、$c$も１００倍することで描写する$\Omega$は変わらないのですが、$beta$に対して$\vec{d_i}$が十分小さいと、$L$は$f$から$g$に変えることでほぼ１００倍変わります。

つまり、$\beta$の他に、$f$のオーダーも考えたほうがいいかもしれません。
ですが適当な$\beta$でやってみたら案外うまくいったので、とりあえずそれでの結果を最初にお見せしました。

他にも、「原点から距離$L_0$より小さい範囲にしかない」というのがわかっているのなら、$1.5\times L_0$まではただの半径$L_0$の球だとして普通のレイマーチングの距離関数をつかうなども有効かもしれません。

実際十分近づくまではアタリをつけた球面として扱うというのは、こういった波動関数を描写するのに有効かと思われます。というのは、こういった波動関数は当たり前のようにgaussian（$e^{-x^2}$）で書かれています。これが数回あるだけでもWebGLとかで計算させるとか正気かとなりますが、実用的な波動関数ではこれが一回の$f$の計算で数十、数百回でてきます。なので、スカラー場の計算自体を減らすという点でも球殻モデルでサボる意義はあるんじゃないかなと。

ていうか、$f(x, y, z)$ を負荷の低い別の関数で近似するというのも必要かなと思ってきた。どうせこういうのは大まかな形が見れたらいいので、数値精度はそれほど気にしませんし。残念ながら良い近似関数を知らないという問題がありますが。

ちなみに、なんか複雑な物体を球殻で覆い、その球殻に達してから複雑な物体を評価するというのは、その道の人でもなされている工夫だったようです。
https://qiita.com/archeleeds/items/b3c8735c3e9ee4a9a5b7

失敗例

これらはいずれもパラメータの調整でなんとかなることもあります。しかし、多くのスカラー場に有効な汎用的なパラメータ、あるいはアルゴリズムがあるかと言われると、難しいかなと思います。

近くを横切る

目的の平面の近くを横切るときは、距離関数中の分数の分子に来る$f(\vec{r})-c$の値が０に近いので、距離関数の返す値が非常に小さくなり、ほとんど光線が進まないことがあります。その結果イテレーション回数上限に達してしまい、光がどの表面にも当たらず無限遠に行ったという判定になることがあります。まあこれは普通の距離関数でも起こりうることではあるのですが、今回は特にこの距離関数が実際の距離ではなくてスカラー場の値の差に依存するため、距離の評価が普通の表面と異なってしまい、問題が起こることがあります。

$c = 0.1$と小さめにしてみました。

edge_error.gif

境界が背景色になってしまっているのが解ると思います。

棘などの細かな構造

$c = 0.8$からはじめ、動画途中で$c = 0.7$に変えてみています。

spike_error.gif

ご覧のように、尖ったところは構造を通り越してしまうことがあります。

おまけ

光線が進むときに、距離を評価するのではなく、その場での波動関数の値に応じて色を重ねていくことを考えます。霧の中を進んでいく感じですね。

ソースコード中での render_wavefunc_fog.frag がこのおまけに相当します。

fog_d.gif

個人的には気に入っております。

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
What you can do with signing up
Sign upLogin
3
1