Ray Marching Distance Functions

[asimes]

I realize this may be an inappropriate place to ask this question as it is not directly fractal related, but in another topic about Ray Marching, ker2x pointed out a website by Iñigo Quilez that shows distance functions for primitives:

Distance Functions: http://www.iquilezles.org/www/articles/distfunctions/distfunctions.htm
Other topic: http://www.fractalforums.com/general-discussion/ray-marching/msg45288/

I've been trying to test out some of these in a general Ray Marching sketch (after having made one for Mandelbulb) but only the sphere makes sense to me. I've just been using the pythagorean theorem like this:

float distFromSphere = sqrt(sq(i-400)+sq(j-200)+sq(ray-200))-100;

Where 'i' and 'j' represent rays shot from the screen and "ray" is the depth of the ray. This makes a sphere at (400, 200, 200) with a radius of 100. If anyone could help me make a cube I would really appreciate it, I can't understand his explanation.

[fractower]

I feel you pain. I am a beginner to DE calculations, and I struggled with this page a while also. The notation is compact to the point of being opaque. On the other hand, a picture is worth 4 lines of code.

Code:

float udBox( vec3 p, vec3 b )
{
  return length(max(abs(p)-b,0.0));
}

b is a vector that describes a box. For example b=(1,1,2) describes a 2x2x4 box at the origin.

abs(p) in this case is an element abs operation. If p = (0, -5, 5) then abs(p) returns (0, 5, 5).

The subtraction is an element subtraction so abs(p) - r  gives (-1, 4, 3).

The max operation is zeros the negative elements giving (0, 4, 3).

The length returns 5 which is the distance from one of the edges of the box.

[eiffie]

It gets easier with time. For me now the simplest way to think of it is the DE for a plane...
float Plane(vec3 z){return z.y;}//if y is negative you've marched too far. so then a rectangular box is...
float Rect(vec3 z, vec3 r){return max(abs(z.x)-r.x, max(abs(z.y)-r.y, abs(z.z)-r.z));}//where r holds the half lengths

I write these with a fourth dimension z.w which starts at 1.0 and is used for scaling so...
float Cylinder(vec4 z, float r, float h){return max(length(z.xz)-r,abs(z.y)-h)/z.w;}//properly returns a scaled result

z*=2.0;
return Cylinder(z,0.25,1.0);//this is useful for difs

It is easy to string these together to make more complicated CSG shapes using min and max.
Then anything is possible...
<a href="http://www.youtube.com/v/Mbu3nRjOSS0&rel=1&fs=1&hd=1" target="_blank">http://www.youtube.com/v/Mbu3nRjOSS0&rel=1&fs=1&hd=1</a>

[Syntopia]

Then anything is possible...

Impressive videos, Eiffie! I'll need to keep an eye on your youtube-stream!

[marius]

Impressive videos, Eiffie! I'll need to keep an eye on your youtube-stream!

Eiffie is the king of fractal mash-up videos. Not to mention ingenious and prolific.
I keep thinking if only I had the spare time, but there's a lot more to it than that.  for eiffie

[asimes]

Thank you very much, I finally have a functional cube. I set it up like this:

Code:

float cubeDist = sqrt(sq(max(abs(i-400)-20, 0))+sq(max(abs(j-150)-20, 0))+sq(max(abs(ray-50)-20, 0)));

I have one more favor to ask, how do I display this correctly with perspective? I wrote this with Processing (Java) but I didn't get any feedback from their forum unfortunately. My sketch works correctly if I display it like this (flat):

Code:

pixVals[i+j*width] += steps;
pixDivs[i+j*width]++;

Here is an image of the flat result (with two spheres behind the cube):


This is my attempt at doing perspective using Processing's built in screenX and screenY. I'm not sure if the screenX / screenY are inaccurate or if I am doing this wrong. I instead display points like this:

Code:

int px = round(screenX(fi, fj, -ray));
int py = round(screenY(fi, fj, -ray));
if (px >= 0 && px < width && py >= 0 && py < height) {
  pixVals[px+py*width] += steps;
  pixDivs[px+py*width]++;
}

Here is an image of the perspective attempt:


If anyone has Processing installed and the patience to look through my code, this should be good to go. Thank you to anyone so patient, and thank you for the responses, I know that this is somewhat out of place for this forum:

Code:

float minThreshold = 0.0001;
int compact = 1;
int maxDepth = 800;

void setup() {
  size(800, 600, P3D);
  noLoop();
}

void draw() {
  loadPixels();
  float[] pixVals = new float[width*height];
  int[] pixDivs = new int[width*height];
  for (int i = 0; i < width*compact; i++) {
    for (int j = 0; j < height*compact; j++) {
      float ray = 0;
      int steps = 255;
      while (steps > 0 && ray < maxDepth) {
        float fi = (float)i/compact;
        float fj = (float)j/compact;
        float currDist1 = sqrt(sq(fi-400)+sq(fj-200)+sq(ray-200))-100;
        float currDist2 = sqrt(sq(fi-400)+sq(fj-300)+sq(ray-100))-100;
        float currDist3 = sqrt(sq(max(abs(fi-400)-20, 0))+sq(max(abs(fj-150)-20, 0))+sq(max(abs(ray-50)-20, 0)));
        float smallestDist = currDist1;
        if (currDist2 < smallestDist) smallestDist = currDist2;
        if (currDist3 < smallestDist) smallestDist = currDist3;
        if (smallestDist*0.5 < minThreshold) {
          int px = round(screenX(fi, fj, -ray));
          int py = round(screenY(fi, fj, -ray));
          if (px >= 0 && px < width && py >= 0 && py < height) {
            pixVals[px+py*width] += steps;
            pixDivs[px+py*width]++;
          }
          break;
        }
        else {
          ray += smallestDist*0.5;
          steps--;
        }
      }
    }
  }
  float highestVal = 0;
  float[] pixAv = new float[width*height];
  for (int i = 0; i < width*height; i++) {
    pixAv[i] = pixVals[i]/pixDivs[i];
    if (pixAv[i] > highestVal) highestVal = pixAv[i];
  }
  for (int i = 0; i < width*height; i++) {
    int val = (int)(pixAv[i]/highestVal*255);
    pixels[i] = (val<<16)+(val<<8)+val;
  }
  updatePixels();
}

[eiffie]

marius - my mom always said I could be the best at anything I set my mind to, as long as no-one else set their's.

CAUTION: I do not have Processing installed but I compiled it in my graphics prog and made some changes so you may have to fix a few bugs. In your version the rays are still cast orthogonally and then you add perspective in post-processing. This reduces your world-space to a long rectangular box (always 800x600). Try to move a sphere in the backbround and it disappears. The answer is to cast the rays along perspective lines. Imagine a ray starting at your eye and going thru a pixel at (i,j). You need to determine how far the ray travels thru x,y and z (I call these dx,dy,dz) as the ray marchs 1 "unit".

Code:

float minThreshold = 0.0001;
int compact = 1;
int maxDepth = 800;

void setup() {
  size(800, 600, P3D);
  noLoop();
}

void draw() {
  loadPixels();
  for (int i = 0; i < width*compact; i++) {
    for (int j = 0; j < height*compact; j++) {
      float dx=(float)i/(float)width*compact-0.5;//creating a perspective ray direction
      float dy=((float)j/(float)height*compact-0.5)*(float)height/(float)width;//aspect ratio
      float dz=1.0,r=1.0/sqrt(dx*dx+dy*dy+1.0);dx*=r;dy*=r;dz*=r;//normalize the ray
      float x=0,y=0,z=0,rayLen=0; //the camera position, these were called fi,fj and ray before
      int steps = 255;
      while (steps > 0 && rayLen < maxDepth) {
        float currDist1 = sqrt(sq(x)+sq(y+20)+sq(z-200))-20;
        float currDist2 = sqrt(sq(x)+sq(y-10)+sq(z-100))-20;
        float currDist3 = max(abs(x)-5,max(abs(y-15)-5, abs(z-50)-5));//got rid of sqrt(sq(...
        float smallestDist = min(currDist1,min(currDist2,currDist3))*0.5;//just simplified
        if (smallestDist < minThreshold) {
          pixels[i+j*width] = (steps<<16)+(steps<<8)+steps; //no need for pixDiv
          break;
        }
        else {
          rayLen += smallestDist;
          x=dx*rayLen;y=dy*rayLen;z=dz*rayLen;//advance the ray
          steps--;
        }
      }
    }
  }
  updatePixels();
}

But why stop there. Now that you have a Processing ray-marcher try making a fractal.

Code:

const float minThreshold = 0.0001,maxDepth=8.0;//8.0 is about right for most fractals

void setup() {
  size(800, 600, P3D);
  noLoop();
}
float DE(float x, float y, float z){//this is our old friend menger
	int n,iters=5;float t;
	for(n=0;n<iters;n++){
		x=fabs(x);y=fabs(y);z=fabs(z);//fabs is just abs for floats
		if(x<y){t=x;x=y;y=t;}
		if(y<z){t=y;y=z;z=t;}
		if(x<y){t=x;x=y;y=t;}
		x=x*3.0-2.0;y=y*3.0-2.0;z=z*3.0-2.0;
		if(z<-1.0)z+=2.0;
	}
	return (sqrt(x*x+y*y+z*z)-1.5)*pow(3.0,-(float)iters);
}
void draw() {
  loadPixels();
  for (int i = 0; i < width; i++) {
    for (int j = 0; j < height; j++) {
      float dx=(float)i/(float)width-0.5;
      float dy=((float)j/(float)height-0.5)*(float)height/(float)width;
      float dz=1.0,r=1.0/sqrt(dx*dx+dy*dy+1.0);dx*=r;dy*=r;dz*=r;
      float x=0,y=0,z=0,rayLen=0,dist=maxDepth; 
      int steps = 255;
      while (steps-- > 0 && dist>=minThreshold && rayLen < maxDepth) {
        rayLen += dist = DE(x,y,z-4.0)*0.5;
        x=dx*rayLen;y=dy*rayLen;z=dz*rayLen;
      }
      if(dist<minThreshold){//eventually you will want to do more here
        //... like find the surface normal and calculate lighting
        pixels[i+j*width] = (steps<<16)+(steps<<8)+steps; 
      }
    }
  }
  updatePixels();
}

[asimes]

Thank you, that makes sense why it didn't work before now. I'm going to play with your first code until I'm comfortable with it. I tried out your Menger code and it looks amazing but I have never made a Menger Sponge myself so I'll get to figuring that one out in a bit. If I could bother you for one last thing, how would you rotate a cube?

Here is an image of your Sponge as it looks in Processing:

[eiffie]

Put the rotation as the first line of your DE function:

float DERotCube(float x, float y, float z, float size, float angle){
  float x2=x*cos(angle)-z*sin(angle),z2=z*cos(angle)+x*sin(angle);//rotating around the y axis
  return max(fabs(x2),max(fabs(y),fabs(z2)))-size;//return the distance estimate to the cube
}

Eventually as you explore more you will want something faster than Processing. I recommend you look into Fragmentarium at that point. For now tho just have fun!

[asimes]

Thanks, for all the help eiffie, above and beyond the expectations for forum help. Now I just have to get the hang of positioning after the rotation, it is not very intuitive to me yet. It isn't possible to rotate an object from its center as opposed to around an axis is it?

[eiffie]

With ray-marching you have to remember you are rotating the ray not the object so sometimes things seem backwards. But certainly you can translate the center of the object to the origin - rotate around an axis then translate back. Just do everything backwards - lol. If you want to move the object up - move the ray down etc.

[asimes]

I came back to this (it has been a while) and have been experimenting with sphere tracing, this is a test that made use of some of eiffie's code

[eiffie]

The menger looks right but strangely the sphere seems banded. Should be DE=max(SphereDE,MengerDE) right? I haven't done this for a while either.

[asimes]

I see banding mostly just on the ground unless I'm missing something. Yup, it is the max of a sphere and Menger cube