Implicit-Surface Ray-Tracer

 

Bernardo S. Gama

2/20/2006

Lab 2

 

 

Objective

 

The objective of this lab was to write a simplified ray-casting program. The program should be able to read a text file and extract the definition of the scene consisting of light, root extraction iteration and gaussian implicit function parameters.  Once the scene is specified the program should apply numerical algorithms to find the light reflected from the implicit surface(s) and write an image file.

 

Program Overview

 

The core and most CPU time consuming part of the program consists of a root finding and an illumination algorithm. The root finding algorithm returns the root of a sum of implicit gaussian functions specified in a text file and represents the position light rays intercept the implicit surface. This is accomplished by first bracketing the root, i.e. finding an interval where the sign of the function changes, and subsequently applying the bisection and Newton’s method repeatedly.  The second part is an illumination procedure that calculates the angle in which light rays are reflected back by the surface and consequently the intensity of light in each point. This is done by calculating the normal to the surface at each root computed.  The normal is evaluated by taking the gradient of the scalar field that constitutes the surface using central differences to estimate the partial derivatives. The rest of the program consists of a vector class to support the use of vector operations in the previous functions and a procedure to output an image file.

The program was written in C++ exclusively, and compiled using Microsoft Visual C++ 2005 using the compiler full speed optimization switches. A full listing of the program is contained in the appendix.

Windows XP was chosen as the destination platform. The program source code and Windows executable can be found here.

 

Results

 

Below are the images obtained using the sample input files and number of iterations defined at:

http://www.cse.ohio-state.edu/~shareef/cse541/Labs/MetaBalls/

 

 

 

 

The images below were rendered by forcing 5 iterations of the bisection method.

 

 

 

 

The images below were generated by forcing 5 iterations of Newton’s method.

 

 

 

 

Conclusion

 

The images generated with 5 iterations of the bisection method only did not show any sign of deterioration. The images generated with 5 iterations of Newton’s method showed large portions in which the root finding algorithm did not converge. This is rather surprising since Newton’s method is supposed to converge faster than the bisection method. The reason for this is still a mystery at the time.

 

 

Appendix

 

Below are the listings of the principal classes in the program.

 

 

// vector class

// vector functions and operations

 

#pragma once

 

class vector {

 

public:

      double x,y,z;

     

      vector() {

            x = 0.0;

            y = 0.0;

            z = 0.0;

      }

      vector(const vector& v) {     // copy constructor

            x = v.x;

            y = v.y;

            z = v.z;

      }

      vector(double x_val, double y_val, double z_val) {

            x = x_val;

            y = y_val;

            z = z_val;

      }

      vector& operator=(vector &v) {

            x = v.x;

            y = v.y;

            z = v.z;

         return *this; 

      }

      vector operator+(vector &v) { // vector addition

            return vector( x + v.x, y + v.y, z + v.z );

      }

      vector operator-(vector &v){  // vector subtraction

            return vector( x - v.x, y - v.y, z - v.z );

      }

      vector operator*(vector &v) { // cross product

            return vector( y * v.z - z * v.y, z * v.x - x * v.z, x * v.y - y * v.x );

      }

      vector operator*(double c){   // scalar multiplication

            return vector( c * x, c * y, c * z);

      }

      double dot(vector &v) { // dot product

            return ( x * v.x + y * v.y + z * v.z );

      }

      double length2(void) {

            return (x*x + y*y + z*z);

      }

      double length(void) {

            return sqrt(x*x + y*y + z*z);

      }

      vector& norm(void) {    // normalize vector

            double q=sqrt(x*x + y*y + z*z);

            x/=q;

            y/=q;

            z/=q;

            return *this;

      }

// friend functions

      friend vector operator*(double c, vector &v) {

            return vector( c * v.x, c * v.y, c * v.z );

      }

      friend double dot(vector &v1, vector &v2) {     // dot product

            return ( v1.x * v2.x + v1.y * v2.y + v1.z * v2.z );

      }

      friend vector mul(vector &v1, vector &v2) {     // component-wise multiplication

            return vector( v1.x * v2.x, v1.y * v2.y, v1.z * v2.z );

      }

      friend double length2(vector &v) { // length squared

            return ( v.x * v.x + v.y * v.y + v.z * v.z );

      }

      friend double length(vector &v) {   // length of vector

            return sqrt( v.x * v.x + v.y * v.y + v.z * v.z );

      }    

      friend vector norm(vector &v) {     // returns unit vector

            double q=sqrt(v.x * v.x + v.y * v.y + v.z * v.z);

            return vector(v.x / q , v.y / q , v.z/q);

      }

      friend void print(vector &v) {     

            printf ("(%lg, %lg, %lg)\n",v.x,v.y,v.z);

      }    

};

 

// metaball class

// defines a gaussian function

 

#pragma once

 

class metaball

{

public:

      double lambda;

      double sigma;

      vector r0;

 

      metaball () {

            lambda=0.0;

            sigma=0.0;

      }

      metaball(double lambda_val, double sigma_val, vector center) {

            lambda=lambda_val;

            sigma=sigma_val;

            r0=center;

      }

      double gaussian(vector r) {

            return (lambda * exp(-sigma * length2(r - r0)));

      }

      double dx_gaussian(vector r) {

            return (-2.0*lambda*sigma * (r.x - r0.x) * exp(-sigma * length2(r - r0)));

      }

      double dy_gaussian(vector r) {

            return (-2.0*lambda*sigma * (r.y - r0.y) * exp(-sigma * length2(r - r0)));

      }

      double dz_gaussian(vector r) {

            return (-2.0*lambda*sigma * (r.z - r0.z) * exp(-sigma * length2(r - r0)));

      }

};

 

// ray class

// contains the root finding and illumnination algorithms

 

#pragma once

 

struct root {

      double val;

      bool found;

};

 

class ray

{

 

private:

      int imgsz;  // image size

      double amb, diff; // ambient and diffuse reflection terms

      vector light; // Light source direction

      int kbr, kbs, knt; // Bracket, Bisection and Newton method iterations

      int nballs; // number of metaballs defined

      double iso; // iso-contour threshold

      metaball* ball; // metaballs

     

public:

      ray(const wchar_t* fname);

      ~ray(void) {delete[] ball;}

      double f(double x, double y, double z);

      double fdx(double x, double y, double z);

      double fdy(double x, double y, double z);

      double fdz(double x, double y, double z);

      double nfdx(double x, double y, double z);

      double nfdy(double x, double y, double z);

      double nfdz(double x, double y, double z);

      bool findRoot(double x, double y, double &z);

      double illuminate(double x, double y, double z);

      void writeImg(const wchar_t *fname, System::Windows::Forms::ProgressBar^ progressBar1);

};

 

// Ray class function declarations

 

#include "StdAfx.h"

 

// open file and initialize ray class

ray::ray(const wchar_t* fname)

{

      FILE *stream;

      char line[100];

 

      _wfopen_s( &stream, fname, L"r" );

     

      fgets(line, 100, stream);

      sscanf_s(line, "%d",&imgsz);

      fgets(line, 100, stream);

      sscanf_s(line, "%lf %lf",&amb, &diff);

      fgets(line, 100, stream);

      sscanf_s(line, "%lf %lf %lf",&light.x, &light.y, &light.z);

      fgets(line, 100, stream);

      sscanf_s(line, "%d", &kbr);

      fgets(line, 100, stream);

      sscanf_s(line, "%d", &kbs);

      fgets(line, 100, stream);

      sscanf_s(line, "%d", &knt);

      fgets(line, 100, stream);

      sscanf_s(line, "%d", &nballs);

      fgets(line, 100, stream);

      sscanf_s(line, "%lf", &iso);

#ifdef DEBUG     

      printf("%d %lf %lf %lf %lf %lf %d %d %d %d %lf\n", imgsz, amb, diff, light.x, light.y, light.z, kbr, kbs, knt, nballs, iso);

#endif

      ball=new metaball[nballs];    // create the metaballs

      for (int i=0; i<nballs; i++) {

            fgets(line, 100, stream);     // and initialize them

            sscanf_s(line, "%lf %lf %lf %lf %lf", &ball[i].r0.x, &ball[i].r0.y, &ball[i].r0.z, &ball[i].lambda, &ball[i].sigma);

#ifdef DEBUG

            printf("%lf %lf %lf %lf %lf \n", ball[i].r0.x, ball[i].r0.y, ball[i].r0.z, ball[i].lambda, ball[i].sigma);

#endif

      }

}

 

bool ray::findRoot(double x, double y, double &z ) {

     

      double a=-1.0;    // a,b endpoints

      double b=1.0;

      const double epsilon=2.0/(double)(imgsz-1); // root finding tolerance

           

      int i;

     

      // search for root (bracket the root)

      if (kbr) {

     

            double h=(b-a)/kbr;

            double y1=f(x,y,z=a);

            double y2;

            for ( i=0; i<kbr; i++) {

                  //if (y1==0.0) {

                  //    return true;      // exact root found

                  //}

                  y2=f(x,y,z+=h);

                  //if (y2==0.0) {

                  //    return true;      // exact root found

                  //}

                  if (y1*y2<0.0) {

                        a=z-h;

                        b=z;

#ifdef DEBUG

                        printf("a=%lf b=%lf\n",a,b);

#endif

                        break; // root bracketed

                  }

#ifdef DEBUG

                  printf("i=%d y1=%lf y2=%lf\n",i,y1,y2);

#endif

                  y1=y2;

            }

      }

     

      double u=f(x,y,a);

      double v=f(x,y,b);

      if (u * v > 0.0) return false; // unable to bracket root

 

           

      // bisecting (expects root is bracketed)

 

      if (kbs) {

            double w;

            double error=b-a;

            for (i=0; i<kbs; i++) {

                  error/=2.0;

                  z=a+error;

                  w=f(x,y,z);

#ifdef DEBUG

                  printf("i=%d\tz=%lg\tf(x)=%lg\te=%lg\n",i,z,w,error);

#endif

                  if (fabs(error)<epsilon) {

                        return true;      // root converged

                  }

                  if (w * u < 0.0) {

                        b=z;

                        v=w;

                  }

                  else {

                        a=z;

                        u=w;

                  }

            }

      }

      // and finally, use sir Isac Newton method

      if (knt) {

            double d;

            for (i=0; i<knt; i++) {

                  d=f(x,y,z)/fdz(x,y,z);

                  z-=d;

#ifdef DEBUG

                  printf("i=%d z=%lg\n",i,z);

#endif

                  if (fabs(d)<epsilon) {

                        return true;      // root converged

                  }

            }

      }

     

      return false; // root did not converge

}

 

//The implicit function f

double ray::f(double x, double y, double z) {

 

      int i;

      double val=0.0;

     

      for (i=0; i<nballs; i++) {

                  val += ball[i].gaussian(vector(x,y,z));

      }

      return (val - iso);

}

// df/dx

double ray::fdx(double x, double y, double z) {

      int i;

      double val=0.0;

     

      for (i=0; i<nballs; i++) {

                  val += ball[i].dx_gaussian(vector(x,y,z));

      }

      return (val);

}

 

 

// df/dy

double ray::fdy(double x, double y, double z) {

      int i;

      double val=0.0;

     

      for (i=0; i<nballs; i++) {

                  val += ball[i].dy_gaussian(vector(x,y,z));

      }

      return (val);

}

 

// df/dz

double ray::fdz(double x, double y, double z) {

      int i;

      double val=0.0;

     

      for (i=0; i<nballs; i++) {

                  val += ball[i].dz_gaussian(vector(x,y,z));

      }

      return (val);

}

 

// numeric differentiation by central differences

 

double ray::nfdx(double x, double y, double z) {

      const double h=0.001;

      return ((f(x+h, y, z) - f(x-h, y, z)) / (2.0 * h));

}

 

double ray::nfdy(double x, double y, double z) {

      const double h=0.001;

      return ((f(x, y+h, z) - f(x, y-h, z)) / (2.0 * h));

}

 

double ray::nfdz(double x, double y, double z) {

      const double h=0.001;

      return ((f(x, y, z+h) - f(x, y, z-h)) / (2.0 * h));

}

 

double ray::illuminate(double x, double y, double z) {

 

      vector gradient(nfdx(x,y,z), nfdy(x,y,z), nfdz(x,y,z));

      gradient.norm();  // convert gradient to unit length

      double illum=gradient.dot(light);   // find cos of the angle

      if (illum<0.0) illum=0.0;     // normal points away

      illum=amb + diff * illum;     // add the constant ambient to the diffuse;

      //illum=0.5*(1.0-z);    // to see the root location;

      if (illum>1.0) illum=1.0;     // error chk

      if (illum<0.0) illum=0.0;     // error chk

      return illum;

}

 

// method for writing the image

void ray::writeImg(const wchar_t* fname, System::Windows::Forms::ProgressBar^ progressBar1) {

 

      const int nx=imgsz;

      const int ny=imgsz;

      double x,y,z;

      double dx=2.0/(double)(imgsz-1);

      double dy=dx;

      progressBar1->Maximum=imgsz;  // indicates current progress

 

      FILE *stream;

      _wfopen_s(&stream,fname,L"wb");

      fprintf(stream,"P6\n%d %d\n255\n",nx, ny);

      const unsigned char background[3] ={0,0,255};

      light.norm(); // normalize light vector

      y=-1.0;

      for(int iy = 0; iy < ny; iy++, y+=dy) {

            x=-1.0;

            progressBar1->Value=iy; // update progress bar

            for (int ix = 0; ix < nx; ix++, x+=dx) {

                  if (findRoot(x, y, z)) {      // root found ?

                        double intensity = illuminate(x, y, z); // yes, calculate light itensity at root

                        unsigned char gray[3];

                        gray[0]=gray[1]=gray[2]=(int) (255.0 * intensity);

                        fwrite(gray, 1, 3, stream);

                  }

                  else {

                        fwrite(background, 1, 3, stream);   // root not found, use background color

                  }

            }

      }

      fclose(stream);

}