Posted By

mikaelhc on 05/03/10


Tagged


Versions (?)

Distance Estimator for Menger Sponge


 / Published in: C
 

  1. /*
  2.   This is a simple modification of Tom Beddard's MandelbulbQuick.pbk.
  3.   It is based on Knighty's ideas from the Kaleidoscopic IFS fractal:
  4.   http://www.fractalforums.com/3d-fractal-generation/kaleidoscopic-%28escape-time-ifs%29/
  5.  
  6.   Mikael Hvidtfeldt Christensen (Syntopia), May 3rd 2010.
  7. */
  8.  
  9. /**
  10.  * MandelbulbQuick.pbk
  11.  * Last update: 14 December 2009
  12.  *
  13.  * Changelog:
  14.  * 1.0 - Initial release
  15.  * 1.0.1 - Fixed a missing asymmetry thanks to Chris King (http://www.dhushara.com)
  16.  * - Refinements in the colouring
  17.  * 1.0.2 - Added radiolaria option for a funky hair-like effect
  18.  * - Incorporated the scalar derivative method as described here:
  19.  * - http://www.fractalforums.com/mandelbulb-implementation/realtime-renderingoptimisations/
  20.  * 1.0.3 - Created a quick version of the script as using a boolean flag to determine
  21.  * which distance estimation method created long compilation times.
  22.  * 1.0.4 - Fixed issue with older graphic cards and the specular highlights
  23.  *
  24.  *
  25.  * Copyright (c) 2009 Tom Beddard
  26.  * http://www.subblue.com
  27.  *
  28.  * For more Flash and PixelBender based generative graphics experiments see:
  29.  * http://www.subblue.com/blog
  30.  *
  31.  * Licensed under the MIT License:
  32.  * http://www.opensource.org/licenses/mit-license.php
  33.  *
  34.  *
  35.  * Credits and references
  36.  * ======================
  37.  * For the story behind the 3D Mandelbrot see the following page:
  38.  * http://www.skytopia.com/project/fractal/mandelbulb.html
  39.  *
  40.  * The original forum disussion with many implementation details can be found here:
  41.  * http://www.fractalforums.com/3d-fractal-generation/true-3d-mandlebrot-type-fractal/
  42.  *
  43.  * This implementation references the 4D Quaternion GPU Raytracer by Keenan Crane:
  44.  * http://www.devmaster.net/forums/showthread.php?t=4448
  45.  *
  46.  * and the NVIDIA CUDA/OptiX implementation by cbuchner1:
  47.  * http://forums.nvidia.com/index.php?showtopic=150985
  48.  *
  49.  */
  50.  
  51. #define PI 3.141592653
  52. #define MIN_EPSILON 3e-7
  53.  
  54. <languageVersion : 1.0;>
  55.  
  56. kernel MandelbulbQuickDerivative
  57. < namespace : "com.subblue.filters";
  58. vendor : "Tom Beddard";
  59. version : 1;
  60. description : "Mandelbulb Fractal Ray Tracer - the quick version (Kaleidoscopic IFS Mod by Syntopia)";
  61. >
  62.  
  63.  
  64. {
  65. parameter int antialiasing
  66. <
  67. minValue:1;
  68. maxValue:3;
  69. defaultValue:1;
  70. description:"Super sampling quality. Number of samples squared per pixel.";
  71. >;
  72.  
  73. parameter bool phong
  74. <
  75. defaultValue:true;
  76. description: "Enable phong shading.";
  77. >;
  78.  
  79.  
  80. parameter float shadows
  81. <
  82. minValue:0.0;
  83. maxValue:1.0;
  84. defaultValue:0.0;
  85. description: "Enable ray traced shadows.";
  86. >;
  87.  
  88. parameter float ambientOcclusion
  89. <
  90. minValue:0.0;
  91. maxValue:1.0;
  92. defaultValue:0.8;
  93. description: "Enable fake ambient occlusion factor based on the orbit trap.";
  94. >;
  95.  
  96. parameter float ambientOcclusionEmphasis
  97. <
  98. minValue:0.0;
  99. maxValue:1.0;
  100. defaultValue:0.58;
  101. description: "Emphasise the structure edges based on the number of steps it takes to reach a point in the fractal.";
  102. >;
  103.  
  104. parameter float bounding
  105. <
  106. minValue:1.0;
  107. maxValue:16.0;
  108. defaultValue:12.5;
  109. description: "Sets the bounding sphere radius to help accelerate the raytracing.";
  110. >;
  111.  
  112. parameter float bailout
  113. <
  114. minValue:0.5;
  115. maxValue:12.0;
  116. defaultValue:4.0;
  117. description: "Sets the bailout value for the fractal calculation. Lower values give smoother less detailed results.";
  118. >;
  119.  
  120. parameter float power
  121. <
  122. minValue:2.0;
  123. maxValue:32.0;
  124. defaultValue:8.0;
  125. description: "The power of the fractal.";
  126. >;
  127.  
  128.  
  129. parameter float3 cameraPosition
  130. <
  131. minValue:float3(-4, -4, -4);
  132. maxValue:float3(4, 4, 4);
  133. defaultValue:float3(0, -2.88, 0);
  134. description: "Camera position.";
  135. >;
  136.  
  137. parameter float3 cameraPositionFine
  138. <
  139. minValue:float3(-0.1, -0.1, -0.1);
  140. maxValue:float3(0.1, 0.1, 0.1);
  141. defaultValue:float3(0, 0.0, 0.0);
  142. description: "Fine tune position.";
  143. >;
  144.  
  145. parameter float3 cameraRotation
  146. <
  147. minValue:float3(-180, -180, -180);
  148. maxValue:float3(180, 180, 180);
  149. defaultValue:float3(-68.4, 7.2, -90);
  150. description: "Pointing angle in each axis of the camera.";
  151. >;
  152.  
  153. parameter float cameraZoom
  154. <
  155. minValue:0.0;
  156. maxValue:10.0;
  157. defaultValue:0.0;
  158. description: "Zoom the camera view.";
  159. >;
  160.  
  161. parameter float3 light
  162. <
  163. minValue:float3(-50, -50, -50);
  164. maxValue:float3(50, 50, 50);
  165. defaultValue:float3(38, -42, 38);
  166. description: "Position of point light.";
  167. >;
  168.  
  169. parameter float3 colorBackground
  170. <
  171. minValue:float3(0, 0, 0);
  172. maxValue:float3(1, 1, 1);
  173. defaultValue:float3(0.0, 0.0, 0.0);
  174. description: "Background colour.";
  175. aeUIControl: "aeColor";
  176. >;
  177.  
  178. parameter float colorBackgroundTransparency
  179. <
  180. minValue:float(0.0);
  181. maxValue:float(1.0);
  182. defaultValue:float(1.0);
  183. description: "Background transparency.";
  184. >;
  185.  
  186. parameter float3 colorDiffuse
  187. <
  188. minValue:float3(0, 0, 0);
  189. maxValue:float3(1, 1, 1);
  190. defaultValue:float3(0.0, 0.85, 0.99);
  191. description: "Diffuse colour.";
  192. aeUIControl: "aeColor";
  193. >;
  194.  
  195. parameter float3 colorAmbient
  196. <
  197. minValue:float3(0, 0, 0);
  198. maxValue:float3(1, 1, 1);
  199. defaultValue:float3(0.67, 0.85, 1.0);
  200. description: "Ambient light colour.";
  201. aeUIControl: "aeColor";
  202. >;
  203.  
  204. parameter float colorAmbientIntensity
  205. <
  206. minValue:float(0);
  207. maxValue:float(1);
  208. defaultValue:float(0.4);
  209. description: "Ambient light intensity.";
  210. >;
  211.  
  212. parameter float3 colorLight
  213. <
  214. minValue:float3(0, 0, 0);
  215. maxValue:float3(1, 1, 1);
  216. defaultValue:float3(0.48, 0.59, 0.66);
  217. description: "Light colour.";
  218. aeUIControl: "aeColor";
  219. >;
  220.  
  221. parameter float colorSpread
  222. <
  223. minValue:0.0;
  224. maxValue:1.0;
  225. defaultValue:0.2;
  226. description: "Vary the colour based on the normal direction.";
  227. >;
  228.  
  229. parameter float rimLight
  230. <
  231. minValue:0.0;
  232. maxValue:1.0;
  233. defaultValue:0.0;
  234. description: "Rim light factor.";
  235. >;
  236.  
  237. parameter float specularity
  238. <
  239. minValue:0.0;
  240. maxValue:1.0;
  241. defaultValue:0.66;
  242. description: "Phone specularity";
  243. >;
  244.  
  245. parameter float specularExponent
  246. <
  247. minValue:0.1;
  248. maxValue:50.0;
  249. defaultValue:15.0;
  250. description: "Phong shininess";
  251. >;
  252.  
  253. parameter float3 rotation
  254. <
  255. minValue:float3(-180, -180, -180);
  256. maxValue:float3(180, 180, 180);
  257. defaultValue:float3(0, 36, 39.6);
  258. description: "Rotate the Mandelbulb in each axis.";
  259. >;
  260.  
  261. parameter int maxIterations
  262. <
  263. minValue:1;
  264. maxValue:20;
  265. defaultValue:6;
  266. description: "More iterations reveal more detail in the fractal surface but takes longer to calculate.";
  267. >;
  268.  
  269. parameter int stepLimit
  270. <
  271. minValue:10;
  272. maxValue:2000;
  273. defaultValue:310;
  274. description: "The maximum number of steps a ray should take.";
  275. >;
  276.  
  277. parameter float epsilonScale
  278. <
  279. minValue:0.1;
  280. maxValue:1.0;
  281. defaultValue:1.0;
  282. description: "Scale the epsilon step distance. Smaller values are slower but will generate smoother results for thin areas.";
  283. >;
  284.  
  285. parameter int2 size
  286. <
  287. minValue:int2(100, 100);
  288. maxValue:int2(2048, 2048);
  289. defaultValue:int2(640, 480);
  290. description: "The output size in pixels.";
  291. >;
  292.  
  293. region generated()
  294. {
  295. return region(float4(0, 0, size.x, size.y));
  296. }
  297.  
  298. dependent float bailout_sr, aspectRatio, sampleStep, sampleContribution, pixel_scale, eps_start;
  299. dependent float3 eye;
  300. dependent float3x3 viewRotation, objRotation;
  301. output pixel4 dst;
  302.  
  303.  
  304. // The fractal calculation
  305. //
  306. // Calculate the closest distance to the fractal boundary and use this
  307. // distance as the size of the step to take in the ray marching.
  308. //
  309. float DE(float3 z0, inout float min_dist)
  310. {
  311. float r = length(z0);
  312. float scale = 3.1;
  313. float CX = 1.0;
  314. float CY = 1.0;
  315. float CZ = 1.0;
  316. float t = 0.0;
  317. int i = 0;
  318. for (i=0;i<maxIterations && r<60.0;i++){
  319.  
  320. // Rotation around z-axis
  321. float angle = -0.1;
  322. float x2 = cos(angle)*z0.x + sin(angle)*z0.y;
  323. float y2 = -sin(angle)*z0.x + cos(angle)*z0.y;
  324. z0.x = x2; z0.y = y2;
  325.  
  326.  
  327. z0.x=abs( z0.x); z0.y=abs( z0.y); z0.z=abs( z0.z);
  328. if( z0.x- z0.y<0.0){t= z0.y; z0.y= z0.x; z0.x=t;}
  329. if( z0.x- z0.z<0.0){t= z0.z; z0.z= z0.x; z0.x=t;}
  330. if( z0.y- z0.z<0.0){t= z0.z; z0.z= z0.y; z0.y=t;}
  331.  
  332. // Reverse rotation around z-axis
  333. angle = -angle;
  334. x2 = cos(angle)*z0.x + sin(angle)*z0.y;
  335. y2 = -sin(angle)*z0.x + cos(angle)*z0.y;
  336. z0.x = x2; z0.y = y2;
  337.  
  338. z0.x=scale* z0.x-CX*(scale-1.0);
  339. z0.y=scale* z0.y-CY*(scale-1.0);
  340. z0.z=scale* z0.z;
  341. if( z0.z>0.5*CZ*(scale-1.0)) z0.z-=CZ*(scale-1.0);
  342.  
  343.  
  344. // Extra folding step
  345. if(z0.x+z0.y<1.0){t=-z0.y;z0.y=-z0.x;z0.x=t;}
  346. if(z0.x+z0.z<0.0){t=-z0.z;z0.z=-z0.x;z0.x=t;}
  347. if(z0.y+z0.z<0.0){t=-z0.z;z0.z=-z0.y;z0.y=t;}
  348.  
  349.  
  350. r=length(z0);
  351. }
  352. return (sqrt(r)-2.0)*pow(scale,float(-i));//the estimated distance
  353. }
  354.  
  355.  
  356. // Intersect bounding sphere
  357. //
  358. // If we intersect then set the tmin and tmax values to set the start and
  359. // end distances the ray should traverse.
  360. bool intersectBoundingSphere(float3 origin,
  361. float3 direction,
  362. out float tmin,
  363. out float tmax)
  364. {
  365. bool hit = false;
  366.  
  367. float b = dot(origin, direction);
  368. float c = dot(origin, origin) - bounding;
  369. float disc = b*b - c; // discriminant
  370. tmin = tmax = 0.0;
  371.  
  372. if (disc > 0.0) {
  373. // Real root of disc, so intersection
  374. float sdisc = sqrt(disc);
  375. float t0 = -b - sdisc; // closest intersection distance
  376. float t1 = -b + sdisc; // furthest intersection distance
  377.  
  378. if (t0 >= 0.0) {
  379. // Ray intersects front of sphere
  380. float min_dist;
  381. float3 z = origin + t0 * direction;
  382. tmin = DE(z, min_dist);
  383. tmax = t0 + t1;
  384. } else if (t0 < 0.0) {
  385. // Ray starts inside sphere
  386. float min_dist;
  387. float3 z = origin;
  388. tmin = DE(z, min_dist);
  389. tmax = t1;
  390. }
  391. hit = true;
  392. }
  393.  
  394. return hit;
  395. }
  396.  
  397.  
  398. // Calculate the gradient in each dimension from the intersection point
  399. float3 estimate_normal(float3 z, float e)
  400. {
  401. float min_dst; // Not actually used in this particular case
  402. float3 z1 = z + float3(e, 0, 0);
  403. float3 z2 = z - float3(e, 0, 0);
  404. float3 z3 = z + float3(0, e, 0);
  405. float3 z4 = z - float3(0, e, 0);
  406. float3 z5 = z + float3(0, 0, e);
  407. float3 z6 = z - float3(0, 0, e);
  408.  
  409. float dx = DE(z1, min_dst) - DE(z2, min_dst);
  410. float dy = DE(z3, min_dst) - DE(z4, min_dst);
  411. float dz = DE(z5, min_dst) - DE(z6, min_dst);
  412.  
  413. return normalize(float3(dx, dy, dz) / (2.0*e));
  414. }
  415.  
  416.  
  417. // Computes the direct illumination for point pt with normal N due to
  418. // a point light at light and a viewer at eye.
  419. float3 Phong(float3 pt, float3 N, out float specular)
  420. {
  421. float3 diffuse = float3(0, 0, 0); // Diffuse contribution
  422. float3 color = float3(0, 0, 0);
  423. specular = 0.0;
  424.  
  425. float3 L = normalize(light * objRotation - pt); // find the vector to the light
  426. float NdotL = dot(N, L); // find the cosine of the angle between light and normal
  427.  
  428. if (NdotL > 0.0) {
  429. // Diffuse shading
  430. diffuse = colorDiffuse + abs(N) * colorSpread;
  431. diffuse *= colorLight * NdotL;
  432.  
  433. // Phong highlight
  434. float3 E = normalize(eye - pt); // find the vector to the eye
  435. float3 R = L - 2.0 * NdotL * N; // find the reflected vector
  436. float RdE = dot(R,E);
  437.  
  438. if (RdE <= 0.0) {
  439. specular = specularity * pow(abs(RdE), specularExponent);
  440. }
  441. } else {
  442. diffuse = colorDiffuse * abs(NdotL) * rimLight;
  443. }
  444.  
  445. return (colorAmbient * colorAmbientIntensity) + diffuse;
  446. }
  447.  
  448.  
  449. // Define the ray direction from the pixel coordinates
  450. float3 rayDirection(float2 p)
  451. {
  452. float3 direction = float3( 2.0 * aspectRatio * p.x / float(size.x) - aspectRatio,
  453. -2.0 * p.y / float(size.y) + 1.0,
  454. -2.0 * exp(cameraZoom));
  455. return normalize(direction * viewRotation * objRotation);
  456. }
  457.  
  458.  
  459. // Calculate the output colour for each input pixel
  460. float4 renderPixel(float2 pixel)
  461. {
  462. float tmin, tmax;
  463. float3 ray_direction = rayDirection(pixel);
  464. float4 color;
  465. color.rgb = colorBackground.rgb;
  466. color.a = colorBackgroundTransparency;
  467.  
  468. if (intersectBoundingSphere(eye, ray_direction, tmin, tmax)) {
  469. float3 ray = eye + tmin * ray_direction;
  470.  
  471. float dist, ao;
  472. float min_dist = 4.0;
  473. float ray_length = tmin;
  474. float eps = MIN_EPSILON;
  475.  
  476. // number of raymarching steps scales inversely with factor
  477. int max_steps = int(float(stepLimit) / epsilonScale);
  478. int i;
  479. float f;
  480.  
  481. for (i = 0; i < max_steps; ++i) {
  482. dist = DE(ray, min_dist);
  483.  
  484. // March ray forward
  485. f = epsilonScale * dist;
  486. ray += f * ray_direction;
  487. ray_length += f * dist;
  488.  
  489. // Are we within the intersection threshold or completely missed the fractal
  490. if (dist < eps || ray_length > tmax) break;
  491.  
  492. // Set the intersection threshold as a function of the ray length from the camera
  493. eps = max(MIN_EPSILON, pixel_scale * ray_length);
  494. }
  495.  
  496.  
  497. // Found intersection?
  498. if (dist < eps) {
  499. ao = 1.0 - clamp(1.0 - min_dist * min_dist, 0.0, 1.0) * ambientOcclusion;
  500.  
  501. if (phong) {
  502. float3 normal = estimate_normal(ray, eps/2.0);
  503. float specular = 0.0;
  504. color.rgb = Phong(ray, normal, specular);
  505.  
  506. if (shadows > 0.0) {
  507. // The shadow ray will start at the intersection point and go
  508. // towards the point light. We initially move the ray origin
  509. // a little bit along this direction so that we don't mistakenly
  510. // find an intersection with the same point again.
  511. float3 light_direction = normalize((light - ray) * objRotation);
  512. ray += normal * eps * 2.0;
  513.  
  514. float min_dist2;
  515. dist = 4.0;
  516.  
  517. for (int j = 0; j < max_steps; ++j) {
  518. dist = DE(ray, min_dist2);
  519.  
  520. // March ray forward
  521. f = epsilonScale * dist;
  522. ray += f * light_direction;
  523.  
  524. // Are we within the intersection threshold or completely missed the fractal
  525. if (dist < eps || dot(ray, ray) > bounding * bounding) break;
  526. }
  527.  
  528. // Again, if our estimate of the distance to the set is small, we say
  529. // that there was a hit and so the source point must be in shadow.
  530. if (dist < eps) {
  531. color.rgb *= 1.0 - shadows;
  532. } else {
  533. // Only add specular component when there is no shadow
  534. color.rgb += specular;
  535. }
  536. } else {
  537. color.rgb += specular;
  538. }
  539. } else {
  540. // Just use the base colour
  541. color.rgb = colorDiffuse;
  542. }
  543.  
  544. ao *= 1.0 - (float(i) / float(max_steps)) * ambientOcclusionEmphasis * 2.0;
  545. color.rgb *= ao;
  546. color.a = 1.0;
  547. }
  548. }
  549.  
  550. return clamp(color, 0.0, 1.0);
  551. }
  552.  
  553.  
  554. // Common values used by all pixels
  555. void evaluateDependents()
  556. {
  557. aspectRatio = float(size.x) / float(size.y);
  558.  
  559. // Camera orientation
  560. float c1 = cos(radians(-cameraRotation.x));
  561. float s1 = sin(radians(-cameraRotation.x));
  562. float3x3 viewRotationY = float3x3( c1, 0, s1,
  563. 0, 1, 0,
  564. -s1, 0, c1);
  565.  
  566. float c2 = cos(radians(-cameraRotation.y));
  567. float s2 = sin(radians(-cameraRotation.y));
  568. float3x3 viewRotationZ = float3x3( c2, -s2, 0,
  569. s2, c2, 0,
  570. 0, 0, 1);
  571.  
  572. float c3 = cos(radians(-cameraRotation.z));
  573. float s3 = sin(radians(-cameraRotation.z));
  574. float3x3 viewRotationX = float3x3( 1, 0, 0,
  575. 0, c3, -s3,
  576. 0, s3, c3);
  577.  
  578. viewRotation = viewRotationX * viewRotationY * viewRotationZ;
  579.  
  580. // Object rotation
  581. c1 = cos(radians(-rotation.x));
  582. s1 = sin(radians(-rotation.x));
  583. float3x3 objRotationY = float3x3( c1, 0, s1,
  584. 0, 1, 0,
  585. -s1, 0, c1);
  586.  
  587. c2 = cos(radians(-rotation.y));
  588. s2 = sin(radians(-rotation.y));
  589. float3x3 objRotationZ = float3x3( c2, -s2, 0,
  590. s2, c2, 0,
  591. 0, 0, 1);
  592.  
  593. c3 = cos(radians(-rotation.z));
  594. s3 = sin(radians(-rotation.z));
  595. float3x3 objRotationX = float3x3( 1, 0, 0,
  596. 0, c3, -s3,
  597. 0, s3, c3);
  598.  
  599. objRotation = objRotationX * objRotationY * objRotationZ;
  600.  
  601. eye = (cameraPosition + cameraPositionFine);
  602. if (eye == float3(0, 0, 0)) eye = float3(0, 0.0001, 0);
  603.  
  604. eye *= objRotation;
  605.  
  606. // Super sampling
  607. sampleStep = 1.0 / float(antialiasing);
  608. sampleContribution = 1.0 / pow(float(antialiasing), 2.0);
  609. pixel_scale = 1.0 / max(float(size.x), float(size.y));
  610. }
  611.  
  612.  
  613. // The main loop
  614. void evaluatePixel()
  615. {
  616. float4 color = float4(0, 0, 0, 0);
  617.  
  618. if (antialiasing > 1) {
  619. // Average antialiasing^2 points per pixel
  620. for (float i = 0.0; i < 1.0; i += sampleStep)
  621. for (float j = 0.0; j < 1.0; j += sampleStep)
  622. color += sampleContribution * renderPixel(float2(outCoord().x + i, outCoord().y + j));
  623. } else {
  624. color = renderPixel(outCoord());
  625. }
  626.  
  627. // Return the final color which is still the background color if we didn't hit anything.
  628. dst = color;
  629. }
  630. }

Report this snippet  

You need to login to post a comment.