HDK
karma_procedurals/mandelbulb_intersect.vfl
/*%vcc -l intersect.otl -O vop % */
#include "math.h"
// Distance estimator adapted from Mikael Hvidtfeldt's blog:
// http://blog.hvidtfeldts.net/index.php/2011/09/distance-estimated-3d-fractals-v-the-mandelbulb-different-de-approximations/~]
float
DE(vector pos; float power; int iters; float bail; float phi_off; float theta_off)
{
vector z = pos;
float dr = 1.0;
float r = 0.0;
float theta, phi, r_pow;
float sin_theta, cos_theta, sin_phi, cos_phi;
for (int i = 0; i < iters; ++i)
{
r = length(z);
if (r > bail)
break;
// convert to polar coords
theta = acos(z.y/r) + theta_off;
phi = atan2(z.z, z.x) + phi_off;
// calcuate scalar derivative approximation
r_pow = pow(r, power-1.0);
dr = r_pow*power*dr + 1.0;
// scale and rotate point
r_pow *= r;
theta *= power;
phi *= power;
// convert back to cartesian coords
sin_theta = sin(theta);
cos_theta = cos(theta);
sin_phi = sin(phi);
cos_phi = cos(phi);
z = r_pow * set(sin_theta*cos_phi,
cos_theta,
sin_theta*sin_phi);
z += pos;
}
return 0.5*log(r)*r/dr;
}
vector
estimateNormal(vector pos; float power; int iters; float bail; float phi; float theta)
{
float eps = 2e-6f;
float n = DE(pos, power, iters, bail, phi, theta);
float dx = DE(pos + set(eps, 0, 0), power, iters, bail, phi, theta);
float dy = DE(pos + set(0, eps, 0), power, iters, bail, phi, theta);
float dz = DE(pos + set(0, 0, eps), power, iters, bail, phi, theta);
vector grad = set(dx, dy, dz) - n;
}
// intersect sphere
int
intersectSphere(vector ray_org; vector ray_dir; float radius; export float dist)
{
int hit = 0;
float b = 2*dot(ray_org, ray_dir);
float d = b*b - 4.0f*c;
// check discriminant for hit
if (d > 0)
{
d = sqrt(d);
float dist0 = 0.5f * (-b - d);
float dist1 = 0.5f * (-b + d);
if (dist0 > 0 && dist0 < dist1)
{
hit = 1;
dist = dist0;
}
else if (dist1 > 0)
{
hit = 1;
dist = dist1;
}
}
return hit;
}
cvex
// Inputs
vector ray_org = 0; // ray origin given by procedural
vector ray_dir = 0; // ray direction given by procedural
// Parameters
int iters = 6;
int max_steps = 100;
float eps = 2e-6f;
float power = 8;
float bailout = 15;
float size = 1;
float phi = 0.0f;
float theta = 0.0f;
// Intersection outputs
export float distance = 0; // distance to hit
export vector Ng = 0; // normal at hit
export vector uvw = 0; // uvw
export int prim_id = 0; // primitive id
export int hit = 0; // boolean determining whether there was a hit
// Attributes
export vector displayColor = { 0, 0, 0 };
)
{
float t = 0.0f;
// this scales everything down to the unit cube
vector pos = ray_org/size;
// intersect with bounding sphere
hit = intersectSphere(pos, ray_dir, 2.0f, t);
pos += t * ray_dir;
distance = t;
// steps
float steps = 0;
// bounds
float xsign = -sign(pos.x);
float ysign = -sign(pos.y);
float zsign = -sign(pos.z);
distance = 0;
// hit bounding sphere
if (hit)
{
hit = 1;
// ray march using the distance estimator
for (steps = 0; steps < max_steps; ++steps)
{
t = DE(pos, power, iters, bailout, phi, theta);
distance += t;
pos += t * ray_dir;
if (t < eps)
{
uvw.y = t;
break;
}
// check if we left the bounding box (sign accounts for original direction)
if (xsign*pos.x < -2 ||
ysign*pos.y < -2 ||
zsign*pos.z < -2)
{
hit = 0;
break;
}
}
}
if (hit)
{
Ng = estimateNormal(pos, power, iters, bailout, phi, theta);
uvw.x = fit(steps, 0, max_steps, 0, 1);