Project Amiga Juggler

Step 0

Programming—especially when you are just starting out, learning how to code—endows you with a sense of control, a sense of power. Often, it’s like playing a God game, at least when you do it recreationally. Of course, you need to spell out to extremely high precision exactly what you want the machine to do. But, then it’s forced to carry out your instructions. It’s forced to obey.

You can derive a similar sense of fun from mathematics. You can force math to do amazing things on your behalf. Of course, you need to figure out the right equations, but once you code them into your programs, the formulae are under your control. They too will obey.

To show you what I mean, let’s do some recreational programming. Let’s recreate the Amiga Juggler from scratch using core Java.

In the fall of 1986, Eric Graham discovered that the Amiga 1000 personal computer, running at a meager 7.16 MHz, was powerful enough to ray trace simple 3D scenes. His test animation consisted of an abstract, humanoid robot juggling 3 mirrored balls. It stunned the Amiga community including the executives at Commodore, who later acquired the rights to the animation for promotional purposes.

I want to share some of the pain that Eric Graham went through by coding this from first principles. To accomplish this, I’ll have to develop a mini ray tracer, a program that generates images by simulating the optical properties of light interacting with an environment. I’ll delve into linear algebra, geometry, trigonometry and a bit of calculus. Hopefully, you vaguely remember learning about these topics back in school. Let’s review.

Step 1

A three-dimensional vector is simply a list of 3 numbers. It can represent a point in space, a direction, a velocity, an acceleration, etc. It can even represent a color in RGB color space. For coding simplicity, I’ll represent a vector as an array of doubles of length 3. An object with 3 double fields is probably more efficient, but I’m trying to minimize the amount of code to write.

Many vector functions act on the individual components independently. For example, to add 2 vectors together, you simply need to add each pair of components together.

  // a = b + c
  public static void add(double[] a, double[] b, double[] c) {
    for(int i = 0; i < 3; i++) {
      a[i= b[i+ c[i];
    }
  }

The magnitude of a vector (the length of a vector) is computing using the 3-dimensional version of the Pythagorean Theorem.

  // |v|
  public static double magnitude(double[] v) {
    return Math.sqrt(v[0* v[0+ v[1* v[1+ v[2* v[2]);
  }

You can scale a vector (change its length) by multiplying it by a constant (a.k.a. a scalar).

  // a = b * s
  public static void scale(double[] a, double[] b, double s) {
    for(int i = 0; i < 3; i++) {
      a[i= b[i* s;
    }
  }

You can negate a vector by scaling by −1, reversing its direction. If you divide a vector by its magnitude (scale by 1/magnitude), you end up with a unit vector (a vector of length 1).

Vectors can do cools things. Consider an astronaut floating in space at some origin point o that fires a gun in direction d (a unit vector). The bullet leaves the gun and it continues to travel indefinitely in a straight line at a constant velocity (a unit of distance for every unit of time). Its position at any time is given by the ray equation.

Above, a vector is scaled by a constant and the result is added to another vector.

  // a = b + c * s
  public static void ray(double[] a, double[] b, double[] c, double s) {
    for(int i = 0; i < 3; i++) {
      a[i= b[i+ c[i* s;
    }
  }

Most vector functions are intuitive, but 2 are mysterious: dot-product and cross-product.

You may recall this trick involving a determinant that is helpful for remembering how to perform cross product:

The notation used in the trick actually reveals a little bit about the historic origin of dot product and cross product. They come from a number system invented by William Hamilton in 1843 called quaternions. Hamilton attempted to extend complex numbers by introducing a new imaginary-like constant j, but he could not work out a self-consistent system. According to legend, while out on a walk in Dublin by the Brougham Bridge, formulae that would work popped into his head. He was so overcome by excitement in his revelation that he impulsively vandalized the bridge by scratching the equations into the stonework (WTF?).

Hamilton’s solution involved introducing 2 new imaginary-like constants:

You can multiply pairs of those constants in 6 different ways. Below, note that the commutative property of multiplication does not work here (order matters).

Using those pairs, it’s possible to multiply 2 quaternions together. When the real components are both 0 (absent), you end up with this:

The real part on the left is the negative dot product and vector part on the right is the cross product. As mathematicians used quaternions in the decades after their discovery, they began to recognize that they were often using them to compute exactly those values. To make life easier, the vector product operators were introduced.

Dot product is an operation that takes place on the plane for which the vectors reside. The image below shows a view of that plane.

These vectors contain only direction information; you can think of their tails as touching the origin. Vector u defines a line and vector v casts a perpendicular shadow onto that line. The problem is find vector p.

You can redraw it as 2 right triangles that share a common side and height.

From there, it’s just a matter of applying the Pythagorean Theorem twice.

In addition to finding p, I wrote down a useful equation that shows that dot product is proportional to the magnitude of the vectors involved and the cosine of the angle between them. Note if the angle is 90 degrees, then the dot product is 0. If the angle is less-than 90, then the dot product is positive. And, if the angle is greater-than 90, then the dot product is negative.

But, what if vector u was a lot shorter than vector v?

As you can see, you end up with exactly the same answer. If the angle is greater than 90 degrees, you can negate one of the vectors to end up with something similar to one of the pictures above. Again, you’ll get the same answer.

Cross product is used to find a vector, w, perpendicular to the plane that u and v reside.

The angle between u and v is not necessarily 90, but the angles between u and w and v and w are both 90. From dot product, that gives us:

If u and v reside on a plane, as opposed to a line, then they cannot be parallel. But, if they were parallel, then they would only be distinguishable by some constant k.

If they are not parallel, then at least one of the following relationships cannot exist.

Suppose that the first or the second relationship is false. Then, we can rewrite the 2 original equations like this:

Next, we solve for g and h using Cramer’s Rule. Note that if the first and the second of the aforementioned relationships were both true, then we would be dividing by 0.

We end up with this:

Now, you have some freedom here. You can actually set i to any value of your choosing and compute g and h accordingly. The resulting vector, w, will be perpendicular to the uv-plane. But, you may happen to notice that the (non-zero) denominators of both fractions in the formulae are the same. If you set i to that value, the denominators vanish and the result is the definition of cross product.

If the first and the second of the aforementioned relationships were true, then you could make a similar argument using the second and the third instead since at least one of them must be false assuming the vectors are not parallel. When they are parallel, the cross product is a null vector.

One important aspect of cross product (not proved here) is that it obeys the right-hand rule.

Finally, using all these vector ideas, I put together the following utility class.

public final class Vec {

  private Vec() {
  }

  // a += b
  public static void add(double[] a, double[] b) {
    for(int i = 0; i < 3; i++) {
      a[i+= b[i];
    }
  }

  // a = b + c
  public static void add(double[] a, double[] b, double[] c) {
    for(int i = 0; i < 3; i++) {
      a[i= b[i+ c[i];
    }
  }

  // a -= b
  public static void subtract(double[] a, double[] b) {
    for(int i = 0; i < 3; i++) {
      a[i-= b[i];
    }
  }

  // a = b - c
  public static void subtract(double[] a, double[] b, double[] c) {
    for(int i = 0; i < 3; i++) {
      a[i= b[i- c[i];
    }
  }

  // a *= s
  public static void scale(double[] a, double s) {
    for(int i = 0; i < 3; i++) {
      a[i*= s;
    }
  }

  // a = b * s
  public static void scale(double[] a, double[] b, double s) {
    for(int i = 0; i < 3; i++) {
      a[i= b[i* s;
    }
  }

  // a /= s
  public static void divide(double[] a, double s) {
    double inverse = 1.0 / s;
    for(int i = 0; i < 3; i++) {
      a[i*= inverse;
    }
  }

  // a = b / s
  public static void divide(double[] a, double[] b, double s) {
    double inverse = 1.0 / s;
    for(int i = 0; i < 3; i++) {
      a[i= b[i* inverse;
    }
  }

  // a = b
  public static void assign(double[] a, double[] b) {
    for(int i = 0; i < 3; i++) {
      a[i= b[i];
    }
  }

  // a = (x, y, z)
  public static void assign(double[] a, double x, double y, double z) {
    a[0= x;
    a[1= y;
    a[2= z;
  }

  // |v|
  public static double magnitude(double[] v) {
    return Math.sqrt(v[0* v[0+ v[1* v[1+ v[2* v[2]);
  }

  // |v|^2
  public static double magnitude2(double[] v) {
    return v[0* v[0+ v[1* v[1+ v[2* v[2];
  }

  // |a - b|
  public static double distance(double[] a, double[] b) {
    double x = a[0- b[0];
    double y = a[1- b[1];
    double z = a[2- b[2];
    return Math.sqrt(x * x + y * y + z * z);
  }

  // |a - b|^2
  public static double distance2(double[] a, double[] b) {
    double x = a[0- b[0];
    double y = a[1- b[1];
    double z = a[2- b[2];
    return x * x + y * y + z * z;
  }

  // a = a / |a|
  public void normalize(double[] a) {
    double s = 1.0 / Math.sqrt(a[0* a[0+ a[1* a[1+ a[2* a[2]);
    for(int i = 0; i < 3; i++) {
      a[i*= s;
    }
  }

  // a = b / |b|
  public void normalize(double[] a, double[] b) {
    double s = 1.0 / Math.sqrt(b[0* b[0+ b[1* b[1+ b[2* b[2]);
    for(int i = 0; i < 3; i++) {
      a[i= b[i* s;
    }
  }

  // a = -a
  public static void negate(double[] a) {
    for(int i = 0; i < 3; i++) {
      a[i= -a[i];
    }
  }

  // a = -b
  public static void negate(double[] a, double[] b) {
    for(int i = 0; i < 3; i++) {
      a[i= -b[i];
    }
  }

  // a += b * s
  public static void ray(double[] a, double[] b, double s) {
    for(int i = 0; i < 3; i++) {
      a[i+= b[i* s;
    }
  }

  // a = b + c * s
  public static void ray(double[] a, double[] b, double[] c, double s) {
    for(int i = 0; i < 3; i++) {
      a[i= b[i+ c[i* s;
    }
  }

  // a . b
  public static double dot(double[] a, double[] b) {
    return a[0* b[0+ a[1* b[1+ a[2* b[2];
  }

  // b = b x c
  public static void cross(double[] b, double[] c) {
    double x = b[1* c[2- b[2* c[1];
    double y = b[2* c[0- b[0* c[2];
    double z = b[0* c[1- b[1* c[0];
    b[0= x;
    b[1= y;
    b[2= z;
  }

  // a = b x c
  public static void cross(double[] a, double[] b, double[] c) {
    double x = b[1* c[2- b[2* c[1];
    double y = b[2* c[0- b[0* c[2];
    double z = b[0* c[1- b[1* c[0];
    a[0= x;
    a[1= y;
    a[2= z;  
  }
}

It’s unfortunate that Java does not support operator overloading. Instead, this is an all static class like the Math class.

Step 2

The original Amiga Juggler animation was a 1 second loop consisting of 24 frames rendered at 320×200 resolution. I am aiming for 30 frames rendered at the HD resolution of 1920×1080. However, during testing, I’ll render 1 frame at a third of that resolution (640×360). To that end, I need some boilerplate code to save an image.

import java.awt.image.*;
import javax.imageio.*;
import java.io.*;
import java.util.*;

public class Main {

  public static final int SQRT_SAMPLES = 1;

  public static final int IMAGE_SCALE = 3;
  public static final int WIDTH = 1920 / IMAGE_SCALE;
  public static final int HEIGHT = 1080 / IMAGE_SCALE;

  public static final String IMAGE_TYPE = "jpg";
  public static final String OUTPUT_FILE = "output." + IMAGE_TYPE;

  public static final double GAMMA = 2.2;

  public static final int SAMPLES = SQRT_SAMPLES * SQRT_SAMPLES;
  public static final double INVERSE_SAMPLES = 1.0 / SAMPLES;
  public static final double INVERSE_GAMMA = 1.0 / GAMMA;

  public static final long SECOND_MILLIS = 1000L;
  public static final long MINUTE_MILLIS = 60 * SECOND_MILLIS;
  public static final long HOUR_MILLIS = 60 * MINUTE_MILLIS;

  private double[][][] pixels = new double[HEIGHT][WIDTH][3];

  public void launch() throws Throwable {
    saveImage();
  }

  private void saveImage() throws Throwable {
    int[] data = new int[WIDTH * HEIGHT];

    for(int y = 0, k = 0; y < HEIGHT; y++) {
      for(int x = 0; x < WIDTH; x++, k++) {
        int value = 0;
        for(int i = 0; i < 3; i++) {
          int intensity = (int)Math.round(255
              * Math.pow(pixels[y][x][i* INVERSE_SAMPLES, INVERSE_GAMMA));
          if (intensity < 0) {
            intensity = 0;
          else if (intensity > 255) {
            intensity = 255;
          }
          value <<= 8;
          value |= intensity;
        }
        data[k= value;
      }
    }

    BufferedImage image = new BufferedImage(
        WIDTH, HEIGHT, BufferedImage.TYPE_INT_RGB);
    image.setRGB(00, WIDTH, HEIGHT, data, 0, WIDTH);
    ImageIO.write(image, IMAGE_TYPE, new File(OUTPUT_FILE));
  }

  public static void main(String... argsthrows Throwable {

    long startTime = System.currentTimeMillis();

    Main main = new Main();
    main.launch();

    long interval = System.currentTimeMillis() - startTime;
    long hours = interval / HOUR_MILLIS;
    interval %= HOUR_MILLIS;
    long minutes = interval / MINUTE_MILLIS;
    interval %= MINUTE_MILLIS;
    long seconds = interval / SECOND_MILLIS;
    interval %= SECOND_MILLIS;
    System.out.format("%d hour%s, %d minute%s, %d second%s, %d millisecond%s%n",
        hours, hours == "" "s",
        minutes, minutes == "" "s",
        seconds, seconds == "" "s",
        interval, interval == "" "s");
  }

}

0 hours, 0 minutes, 0 seconds, 191 milliseconds
source

As you can see, it stores an all black image.

The main method contains code that outputs the rendering time (ray tracing is a slow process).

The pixels array will store the RGB values of each pixel in the image as it is rendered. Each color intensity component is in the range of 0.0 (darkest) to 1.0 (lightest). However, to compute the color value of a pixel, the code will statistically sample that part of the 3D scene a specified number of times (the constant SAMPLES) and it will average the samples together. The pixels array stores the sum of all the samples for each pixel. The saveImage method divides by the number of samples to convert the sum into an average.

For saving, the pixels array is temporarily converted into a BufferedImage. The BufferedImage.setRGB method expects an integer array where color intensities are in the range of 0 to 255. Unfortunately, we cannot simply multiply the pixels array values by 255 because the relationship between them is not linear. The actual brightness of a monitor pixel increases exponentially with the size of that integer value by a power of a constant called gamma (typically 2.2 on PCs and 1.8 on Macs). To compensate, the pixels array values are raised to the inverse of gamma before multiplying.

Step 3

To push this machine to the max, I added code to render in parallel executing threads.

  private int runningCount;
  private int rowIndex;

  public void launch() throws Throwable {
    
    int processors = Runtime.getRuntime().availableProcessors();
    updateRunningCount(processors);
    for(int i = 0; i < processors; i++) {
      new Thread(Integer.toString(i)) {
        @Override
        public void run() {
          render();
          updateRunningCount(-1);
        }
      }.start();
    }

    synchronized(this) {
      while(runningCount != 0) {
        wait();
      }
    }
    
    saveImage();
  }
  
  private void render() {
    while(true) {

      int y = getNextRowIndex();
      if (y >= HEIGHT) {
        return;
      }

      for(int x = 0; x < WIDTH; x++) {
        
      
    }
  }

  private synchronized void updateRunningCount(int dx) {
    runningCount += dx;
    if (runningCount == 0) {
      notifyAll();
    }
  }

  private synchronized int getNextRowIndex() {
    return rowIndex++;
  }

This is a Dell Studio XPS 435T with an Intel Core i7 CPU 920 at 2.67 GHz with 64-bit Vista. Runtime.getRuntime().availableProcessors() returns 8 and I confirmed that it takes 8 threads before the Task Manager reports 100% CPU usage.

Each thread will work on separate rows of the image. The getNextRowIndex method returns the next row of the image yet to be rendered.

Step 4

I decided it would be nice to see the output as it is slowly generated. To make that happen, I created a panel capable of displaying images.

import java.awt.*;
import javax.swing.*;
import java.awt.image.*;

public class ImagePanel extends JPanel {

  private BufferedImage image;

  public ImagePanel(BufferedImage image) {
    this.image = image;
    setPreferredSize(new Dimension(image.getWidth(), image.getHeight()));
  }

  @Override
  protected void paintComponent(Graphics g) {
    g.drawImage(image, 00this);
  }
}

Next, I created a frame to display the panel. The constructor captures a BufferedImage and calling the imageUpdated method forces it to repaint.

import java.awt.*;
import javax.swing.*;
import java.awt.image.*;

public class RenderFrame extends JFrame {

  private ImagePanel imagePanel;

  public RenderFrame(BufferedImage image) {
    setTitle("Amiga Juggler");
    setContentPane(imagePanel = new ImagePanel(image));
    setDefaultCloseOperation(JFrame.EXIT_ON_CLOSE);
    setResizable(false);
    pack();
    setLocationRelativeTo(null);
    setVisible(true);
  }

  public void imageUpdated() {
    EventQueue.invokeLater(new Runnable() {
      @Override
      public void run() {
        imagePanel.repaint();
      }
    });
  }

  @Override
  public void setTitle(final String title) {
    EventQueue.invokeLater(new Runnable() {
      @Override
      public void run() {
        RenderFrame.super.setTitle(title);
      }
    });
  }
}

To support the frame, I had to move some stuff around. There is no longer a global pixels array. Instead, the render method allocates an int array representing color intensities of a single row with RGB values compatible with BufferedImage. The render method does the conversion itself from double intensities to int intensities by dividing by the number of samples to compute the average, applying gamma correction and then multiplying by 255. When it finishes rendering a row, it invokes rowCompleted. That method updates the BufferedImage and it forces the RenderFrame to repaint.

import java.awt.image.*;
import javax.imageio.*;
import java.io.*;
import java.util.*;

public class Main {

  public static final int SQRT_SAMPLES = 1;

  public static final int IMAGE_SCALE = 3;
  public static final int WIDTH = 1920 / IMAGE_SCALE;
  public static final int HEIGHT = 1080 / IMAGE_SCALE;

  public static final String IMAGE_TYPE = "jpg";
  public static final String OUTPUT_FILE = "output." + IMAGE_TYPE;
  public static final boolean RENDER_IN_WINDOW = true;

  public static final double GAMMA = 2.2;

  public static final int SAMPLES = SQRT_SAMPLES * SQRT_SAMPLES;
  public static final double INVERSE_SAMPLES = 1.0 / SAMPLES;
  public static final double INVERSE_GAMMA = 1.0 / GAMMA;

  public static final long SECOND_MILLIS = 1000L;
  public static final long MINUTE_MILLIS = 60 * SECOND_MILLIS;
  public static final long HOUR_MILLIS = 60 * MINUTE_MILLIS;

  private RenderFrame renderFrame;
  private BufferedImage image;
  private int runningCount;
  private int rowIndex;

  public void launch() throws Throwable {

    image = new BufferedImage(WIDTH, HEIGHT, BufferedImage.TYPE_INT_RGB);
    if (RENDER_IN_WINDOW) {
      renderFrame = new RenderFrame(image);
    }
    
    int processors = Runtime.getRuntime().availableProcessors();
    updateRunningCount(processors);
    for(int i = 0; i < processors; i++) {
      new Thread(Integer.toString(i)) {
        @Override
        public void run() {
          render();
          updateRunningCount(-1);
        }
      }.start();
    }

    synchronized(this) {
      while(runningCount != 0) {
        wait();
      }
    }
    
    saveImage();

    if (RENDER_IN_WINDOW) {
      renderFrame.setTitle("Amiga Juggler [DONE]");
    }
  }
  
  private void render() {

    int[] pixels = new int[WIDTH];
    double[] pixel = new double[3];

    while(true) {

      int y = getNextRowIndex();
      if (y >= HEIGHT) {
        return;
      }

      for(int x = 0; x < WIDTH; x++) {

        // -- SIMULATE SLOW DRAWING ---------------------------
        pixel[01;
        try {
          Thread.sleep(1);
        catch(InterruptedException e) {
        }
        // ----------------------------------------------------

        int value = 0;
        for(int i = 0; i < 3; i++) {
          int intensity = (int)Math.round(255
              * Math.pow(pixel[i* INVERSE_SAMPLES, INVERSE_GAMMA));
          if (intensity < 0) {
            intensity = 0;
          else if (intensity > 255) {
            intensity = 255;
          }
          value <<= 8;
          value |= intensity;
        }

        pixels[x= value;
      }

      rowCompleted(y, pixels);
    }
  }

  private synchronized void updateRunningCount(int dx) {
    runningCount += dx;
    if (runningCount == 0) {
      notifyAll();
    }
  }

  private synchronized int getNextRowIndex() {
    return rowIndex++;
  }

  private synchronized void rowCompleted(int rowIndex, int[] pixels) {
    image.setRGB(0, rowIndex, WIDTH, 1, pixels, 0, WIDTH);
    if (RENDER_IN_WINDOW) {
      renderFrame.imageUpdated();
    }
  }

  private void saveImage() throws Throwable {
    ImageIO.write(image, IMAGE_TYPE, new File(OUTPUT_FILE));
  }

  public static void main(String... argsthrows Throwable {

    long startTime = System.currentTimeMillis();

    Main main = new Main();
    main.launch();

    long interval = System.currentTimeMillis() - startTime;
    long hours = interval / HOUR_MILLIS;
    interval %= HOUR_MILLIS;
    long minutes = interval / MINUTE_MILLIS;
    interval %= MINUTE_MILLIS;
    long seconds = interval / SECOND_MILLIS;
    interval %= SECOND_MILLIS;
    System.out.format("%d hour%s, %d minute%s, %d second%s, %d millisecond%s%n",
        hours, hours == "" "s",
        minutes, minutes == "" "s",
        seconds, seconds == "" "s",
        interval, interval == "" "s");
  }

}

To test it out, the inner loop in the render method contains code to slowly set pixels red.


source

Step 5

The virtual camera is based on the camera obscura and the pinhole camera. Surprisingly, a real camera does not actually need a lens to function. A pinhole-sized aperture ensures that each spot on the film is only exposed to light arriving from one particular direction. Unfortunately, in the real world, the smaller the aperture, the longer the exposure time.

Instead of putting virtual film behind the pinhole, a virtual screen is placed in front of it. Geometrically, it works the same way. Rather than a pinhole, consider it as the pupil of a virtual eye staring at the center of the virtual screen.

The ray tracer will follow the pathway of light backwards from the eye, through the virtual screen and out into the scene. In this way, the virtual screen acts like a window. Each pixel in the image has a corresponding point on the virtual screen in 3D space. Before it can cast a ray out into the scene, it has to figure out the direction of the ray.

I should mention that the coordinate system used looks like this:

The position and the orientation of the virtual screen is determined by 3 things: the location of the eye, e; a point that the eye is looking at, l; and, the distance, D, between the eye and the center of the virtual screen, c.

We need to construct an orthonormal basis (ONB) situated at the center of the screen, c.

The unit vectors, u, v and w, are perpendicular to each other. Collectively, the point c and the ONB define a right-handed coordinate system. The point c acts like the origin. The vector u serves as the x-axis of the virtual screen. The vector v acts as the y-axis of the virtual screen. And, the vector w points toward the eye.

If you had some point p = (a, b, e) and you wanted to know where it was located in that coordinate system, you could position it using a variation of the ray equation.

In our case, p represents a point on the virtual screen; so, e = 0.

For the general case, I added 2 new methods to the Vec utility class.

  // p = o + p[0] u + p[1] v + p[2] w
  public static void transform(
      double[] p, double[] o, double[] u, double[] v, double[] w) {
    double x = o[0+ p[0* u[0+ p[1* v[0+ p[2* w[0];
    double y = o[1+ p[0* u[1+ p[1* v[1+ p[2* w[1];
    double z = o[2+ p[0* u[2+ p[1* v[2+ p[2* w[2];
    p[0= x;
    p[1= y;
    p[2= z;    
  }  
  
  // q = o + p[0] u + p[1] v + p[2] w
  public static void transform(
      double[] q, double[] p, double[] o, double[] u, double[] v, double[] w) {
    double x = o[0+ p[0* u[0+ p[1* v[0+ p[2* w[0];
    double y = o[1+ p[0* u[1+ p[1* v[1+ p[2* w[1];
    double z = o[2+ p[0* u[2+ p[1* v[2+ p[2* w[2];
    q[0= x;
    q[1= y;
    q[2= z;
  }

The unit vector w is the easiest of the 3 to obtain:

The formula inspired this addition to the Vec utility:

  // v = ( head - tail ) / | head - tail |
  public static void constructUnitVector(
      double[] v, double[] head, double[] tail) {
    for(int i = 0; i < 3; i++) {
      v[i= head[i- tail[i];
    }
    double s = 1.0 / Math.sqrt(v[0* v[0+ v[1* v[1+ v[2* v[2]);
    for(int i = 0; i < 3; i++) {
      v[i= v[i* s;
    }
  }

There is already a subtract method to construct a non-unit vector from 2 point vectors.

Now that we have w, we can find c by starting at e and traveling a distance of −D along w using the ray equation (w points toward the e).

On the picture of the virtual screen above, observe points r and s. Imagine that the virtual screen were a camera that you were holding with your hands over r and s. You can rotate around (shoot the photo in any direction) and you can tilt the camera upwards and downwards. But, you would typically keep the line between r and s level. In this way, when you take the picture, the horizon will be parallel to the bottom of the photograph. Rarely do you roll the camera because it produces strange photographs.

We’ll construct the remainder of the ONB with the assumption that there is no roll. If you do need to roll the virtual camera, you can rotate u and v around w after the ONB is fully constructed. No camera roll is used in this project.

If there is no roll, then u is perpendicular to the normal of the ground plane that the juggler stands on, n = (0, 1, 0). And, by the definition of the ONB, u is also perpendicular to w. Hence, we can find u using:

Again from the definition of the ONB, we can now solve for v.

You can verify that applying cross products in this sequence does produce a right-handed coordinate system.

There is only 1 catch with the above construction: If w is parallel to the ground normal, n, then taking the cross product produces a null vector. Meaning, if you are looking straight up or straight down, then it will fail. But, it’s clear what u and v should be in those cases.

I added these methods to Vec:

  public static final double TINY = 1e-9;

  // a == b
  public static boolean equals(double[] a, double[] b) {
    for(int i = 0; i < 3; i++) {
      if (Math.abs(a[i- b[i]) > TINY) {
        return false;
      }
    }
    return true;
  }

  // v == (x, y, z)
  public static boolean equals(double[] v, double x, double y, double z) {
    return Math.abs(v[0- x<= TINY
        && Math.abs(v[1- y<= TINY
        && Math.abs(v[2- z<= TINY;
  }

  // w -> u, v
  public static void onb(double[] u, double[] v, double[] w) {
    if (equals(w, 010)) {
      assign(u, 100);
      assign(v, 00, -1);
    else if (equals(w, 0, -10)) {
      assign(u, 100);
      assign(v, 001);
    else {
      assign(u, w[2]0, -w[0]);
      cross(v, w, u);
    }
  }

The input of the onb method is w and the output is u and v.

Putting it all together, I added this to the top of the render method:

    double[] u = new double[3];
    double[] v = new double[3];
    double[] w = new double[3];
    double[] c = new double[3];
    
    Vec.constructUnitVector(w, EYE, LOOK);
    Vec.ray(c, EYE, w, -DISTANCE_TO_VIRTUAL_SCREEN);
    Vec.onb(u, v, w);

Step 6

We normally think of addressing pixels in an image or on the screen in terms of row and column index. The upper-left corner is origin and the y-axis points downwards. Consider this 4×4 image:

The width and height are 4, but we use the indices 0 to 3 (inclusive). Instead of row and column indices, consider the upper-left corner of each pixel as the origin of that pixel. For instance, the upper-left corner of the lower-rightmost pixel is (3, 3). The lower-right corner of that pixel is (4, 4). In this way, it is possible to refer to points anywhere on the image, even within pixels.

Now, consider some arbitrary point on the image that we want to find the equivalent point in 3D space on the virtual screen. The image has a width and height of W and H respectively.

First, without modifying the image, I’ll move the origin to the lower-left and I’ll reverse the direction of the y-axis. To address the same point on the image, the coordinates need to change to:

Again, without modifying the image, I’ll move the origin to the center of the image and adjust the points coordinates to compensate.

Now, let’s assume that the virtual screen has the same aspect ratio as the image. We’ll specify the width of the virtual screen as a constant in the program. For this discussion, I’ll refer to the virtual screen width as A and the virtual screen height as B.

That gives us a way to find the virtual screen height:

If we had some arbitrary point (x, y) in the image aspect ratio and we needed to convert it to the virtual screen aspect ratio, we could accomplish the conversion with this:

Combining the last 2 concepts leaves us with this:

Finally, we have a way to transform some arbitrary 2D point (a, b) on the image to a 3D point on the virtual screen.

That formula inspired these additions to Vec:

  // p = o + p[0] u + p[1] v
  public static void map(
      double[] p, double[] o, double[] u, double[] v) {
    double x = o[0+ p[0* u[0+ p[1* v[0];
    double y = o[1+ p[0* u[1+ p[1* v[1];
    double z = o[2+ p[0* u[2+ p[1* v[2];
    p[0= x;
    p[1= y;
    p[2= z;
  }

  // q = o + p[0] u + p[1] v
  public static void map(
      double[] q, double[] p, double[] o, double[] u, double[] v) {
    double x = o[0+ p[0* u[0+ p[1* v[0];
    double y = o[1+ p[0* u[1+ p[1* v[1];
    double z = o[2+ p[0* u[2+ p[1* v[2];
    q[0= x;
    q[1= y;
    q[2= z;
  }
  
  // q = o + x u + y v
  public static void map(
      double[] q, double[] o, double[] u, double[] v, double x, double y) {
    double i = o[0+ x * u[0+ y * v[0];
    double j = o[1+ x * u[1+ y * v[1];
    double k = o[2+ x * u[2+ y * v[2];
    q[0= i;
    q[1= j;
    q[2= k;
  }

Since the point, p, on the virtual screen is known, we can finally compute the direction of the ray, d:

Step 7

As mentioned above, the ray tracer will fire several rays through each pixel at slightly different coordinates to sample the 3D scene and then it will assign that pixel to the average of the results. To accomplish this, an image pixel is uniformly divided into a grid. Within each grid cell, a random point is chosen. This is known as jittered sampling. There are superior sampling techniques, but this technique works well enough and it’s easy to code. Below is an example of a sampled image pixel.

Generating pseudorandom values may be a computationally expensive operation. To guard against that possibility, I created a utility class that stores a large cache of random values.

public class RandomDoubles {

  private static volatile double[] values;

  private int index;

  public RandomDoubles() {
    synchronized(RandomDoubles.class) {
      if (values == null) {
        values = new double[70001];
        for(int i = values.length - 1; i >= 0; i--) {
          values[i= Math.random();
        }
      }
      index = (int)(Math.random() * values.length);
    }
  }

  public double random() {
    if (index >= values.length) {
      index = 0;
    }
    return values[index++];
  }
}

The cache is initialized when the first instance of RandomDoubles is created. 70001 is a prime number. I chose a prime to reduce the chances of strange patterns from showing up in the output image as a consequence of cyclically using the same set of random values.

Here is the render method after introducing random sampling:

  private void render() {

    double[] u = new double[3];
    double[] v = new double[3];
    double[] w = new double[3];
    double[] c = new double[3];
    double[] p = new double[3];
    double[] d = new double[3];

    Vec.constructUnitVector(w, EYE, LOOK);
    Vec.ray(c, EYE, w, -DISTANCE_TO_VIRTUAL_SCREEN);
    Vec.onb(u, v, w);

    RandomDoubles randomDoubles = new RandomDoubles();
    int[] pixels = new int[WIDTH];
    double[] pixel = new double[3];

    while(true) {

      int y = getNextRowIndex();
      if (y >= HEIGHT) {
        return;
      }

      for(int x = 0; x < WIDTH; x++) {

        for(int i = 0; i < SQRT_SAMPLES; i++) {
          double b = y + INVERSE_SQRT_SAMPLES * i;
          if (SQRT_SAMPLES == 1) {
            b += 0.5;
          else {
            b += randomDoubles.random();
          }

          b = VIRTUAL_SCREEN_RATIO * (HALF_HEIGHT - b);

          for(int j = 0; j < SQRT_SAMPLES; j++) {
            double a = x + INVERSE_SQRT_SAMPLES * j;
            if (SQRT_SAMPLES == 1) {
              a += 0.5;
            else {
              a += randomDoubles.random();
            }

            a = VIRTUAL_SCREEN_RATIO * (a - HALF_WIDTH);

            Vec.map(p, EYE, u, v, a, b);
            Vec.constructUnitVector(d, p, EYE);

            // ... USE d ...
          }
        }

        int value = 0;
        for(int i = 0; i < 3; i++) {
          int intensity = (int)Math.round(255
              * Math.pow(pixel[i* INVERSE_SAMPLES, INVERSE_GAMMA));
          if (intensity < 0) {
            intensity = 0;
          else if (intensity > 255) {
            intensity = 255;
          }
          value <<= 8;
          value |= intensity;
        }

        pixels[x= value;
      }

      rowCompleted(y, pixels);
    }
  }

The square-root of the number of samples, SQRT_SAMPLES, can be increased to improve image quality. Notice that when it is 1, instead of using a random sample, the sample is taken in the center of the pixel.

Step 8

The juggling robot is constructed entirely out of spheres. The sky also appears to be a hemisphere. So, it is pertinent to work out how to intersect a ray and a sphere.

Below, a sphere of radius R is centered about p. A ray originating from o, fired in direction d, strikes the sphere at some point (x, y, z).

Now we have the amount of time it takes for the ray to reach a sphere from the virtual eye. We will apply this formula to every sphere in the scene. If the value within the square-root is negative, then the ray does not intersect the sphere. Of the ones that do intersect, the one with the minimal time is the one that we see. Well almost. If a sphere is behind the eye, this formula will output a negative time value. Intersection only took place if time is positive.

We actually need to constrain time a bit more. The eye is not the only point that we will be casting rays from. For example, when a ray bounces off one of the mirrored spheres, we will cast a new ray off its surface. To prevent the possibility of intersection between a ray and a surface that emitted the ray, we need to make sure that time is greater than a tiny positive value that we’ll call epsilon. Also, we need to make sure that time is smaller than some max value. For example, to determine if a surface is illuminated, we’ll cast a ray back to a light source. In that case, time cannot be greater than the time it takes to reach the light source. Any intersections that occur beyond the source of light cannot be considered.

Once we locate the intersection with minimal time, we’ll compute the color of the hit point using Phong Shading. Phong Shading takes into account the following vectors:

All of the vectors are unit vectors. The vector d is the direction of the incident ray. The negation of d points back to the ray’s origin, back to the virtual eye. The vector n is the surface normal. In case of a sphere, it is the vector between the hit point and the center of the sphere, normalized. However, it is possible that we may need to negate the normal to face the viewer. The angle between n and −d must be less-than-or-equal-to 90 degrees. Meaning, if n dot d is positive, then we reverse the normal. The sky is a good example of this. It is the interior surface of a hemisphere and the normal must point inwards. The vector l is a unit vector pointing toward a light source. And, finally, the vector r is the mirror reflection of the vector l. It represents light originating from the light source and bouncing off of the surface.

Since Phong Shading uses a mirror reflection vector, let’s consider bouncing a ray off of one of the mirror spheres. Below, a ray traveling in direction d strikes a surface at point h with normal n and then it reflects in direction r.

Vector p is the projection of −d onto n (vector projection was illustrated in step 1).

But, for Phong Shading, vector l actually points toward the light source; it’s in the opposite direction that the light ray travels. If we negate the ray direction, we end up with:

Note, if l dot n < 0 (greater than 90 degrees), then the light source is on the opposite side of the surface and none of its light reaches the side that we are viewing. Also, we need to cast a ray back to the light source to determine if there are any objects in the way. If the path to the light source is interrupted, then the surface is in shadow and only ambient light comes into play.

Finally, to shade the surface, we apply the Phong reflection model.

L represents radiance. We need to apply this formula 3 times, once for each RGB color intensity component. There is only one light source in this project. If there were multiple light sources, we would have to apply it to each of them and add them together.

There are 2 dot products in the formula. If one of them returns a negative value, change it to 0. Meaning, if the light source and the normal are on opposite sides, then the light does not contribute to the color of the surface. And, if the reflected light ray cannot be viewed seen, then it does not contribute to the color either.

The 4 terms represent ambient, diffuse, specular and reflected light. Each term is weighted by a k constant that is typically between 0 and 1. Collectively, the constants enable you to change the behavior of the surface.

I created a class to describe the material of a surface.

public class Material {
  public double ambientWeight;
  public double diffuseWeight;
  public double specularWeight;
  public double reflectionWeight;  
  public double shininess;
  public double[] diffuseColor;
  public double[] highlightColor;
  public double[] reflectionColor;

  public Material(
      double ambientWeight,
      double diffuseWeight,
      double specularWeight,
      double reflectionWeight,
      double shininess,
      double[] diffuseColor,
      double[] highlightColor,
      double[] reflectionColor) {
    this.ambientWeight = ambientWeight;
    this.diffuseWeight = diffuseWeight;
    this.specularWeight = specularWeight;
    this.reflectionWeight = reflectionWeight;
    this.shininess = shininess;
    this.diffuseColor = diffuseColor;
    this.highlightColor = highlightColor;
    this.reflectionColor = reflectionColor;    
  }
}

Next, I created a class to house different types of materials. Right now, there is only red plastic. It’s a bright red material that has a white specular highlight.

public final class Materials {

  private Materials() {
  }

  public static final Material RED_PLASTIC = new Material(
    0.1,
    2.0,
    1.0,
    0.0,
    10.0,
    new double[] { 10},
    new double[] { 11},
    new double[] { 00}
  );
}

For ray-object intersection, I created an object called Intersection that represents the result. It contains the hit time, the hit point, the surface normal and the surface material at the hit point.

public class Intersection {
  public double time;
  public double[] hit = new double[3];
  public double[] normal = new double[3];
  public Material material;
}

3D objects in the scene are represented by the following interface:

public interface IObject {
  public boolean intersect(
      double[] o, double[] d, boolean primaryRay, double maxTime,
      double[][] temps,
      Intersection intersection);
}

The intersect method returns true if the ray intersects the 3D object. If the primaryRay parameter is set to true, then additional information about the intersection will be stored in the intersection parameter. The parameters o and d are the origin and direction of the ray respectively. The maxTime parameter is the intersection time upper-bound. The epsilon lower-bound constant is defined in Main. Finally, the temps parameter are temporary vectors that can be used during the intersection calculation. The render method passes in 16 of them.

Here is an implementation using the ray-sphere intersection formula from above. Notice that if primaryRay is not set to true, it does not bother to compute the hit point and the surface normal. When casting shadow rays and when casting ambient occlusion rays, we only care about the return value.

public class Sphere implements IObject {

  public double[] center = new double[3];
  public double radius = 1;
  public Material material;

  public Sphere(
      double x, double y, double z, double radius, Material material) {
    Vec.assign(center, x, y, z);
    this.radius = radius;
    this.material = material;
  }

  public boolean intersect(
      double[] o, double[] d, boolean primaryRay, double maxTime,
      double[][] temps,
      Intersection intersection) {

    Vec.subtract(temps[0], o, center);
    double B = 2.0 * Vec.dot(d, temps[0]);
    double C = Vec.magnitude2(temps[0]) - radius * radius;
    double square = B * B - * C;
    if (square >= 0) {
      double sqrt = Math.sqrt(square);
      double t1 = 0.5 (-B - sqrt);
      double t2 = 0.5 (-B + sqrt);
      boolean intersected = false;
      if (t1 >= Main.EPSILON && t1 <= maxTime) {
        intersected = true;
        intersection.time = t1;
      else if (t2 >= Main.EPSILON && t2 < maxTime) {
        intersected = true;
        intersection.time = t2;
      }
      if (primaryRay && intersected) {
        Vec.ray(intersection.hit, o, d, intersection.time);
        Vec.constructUnitVector(intersection.normal, intersection.hit, center);
        intersection.material = material;
      }
      if (intersected) {
        return true;
      }
    }

    return false;
  }
}

Finally, I updated the render method using the Phong Shading model and I tested it with a scene containing a single sphere.

  private volatile IObject[] scene
      new Sphere(00050, Materials.RED_PLASTIC) };

  private void render() {

    double[] u = new double[3];
    double[] v = new double[3];
    double[] w = new double[3];
    double[] c = new double[3];
    double[] p = new double[3];
    double[] d = new double[3];
    double[] l = new double[3];
    double[] r = new double[3];
    double[][] temps = new double[16][3];
    Intersection intersection = new Intersection();
    Intersection bestIntersection = new Intersection();

    Vec.constructUnitVector(w, EYE, LOOK);
    Vec.ray(c, EYE, w, -DISTANCE_TO_VIRTUAL_SCREEN);
    Vec.onb(u, v, w);

    RandomDoubles randomDoubles = new RandomDoubles();
    int[] pixels = new int[WIDTH];
    double[] pixel = new double[3];

    while(true) {

      int y = getNextRowIndex();
      if (y >= HEIGHT) {
        return;
      }

      for(int x = 0; x < WIDTH; x++) {

        Vec.assign(pixel, 000);

        for(int i = 0; i < SQRT_SAMPLES; i++) {
          double b = y + INVERSE_SQRT_SAMPLES * i;
          if (SQRT_SAMPLES == 1) {
            b += 0.5;
          else {
            b += randomDoubles.random();
          }

          b = VIRTUAL_SCREEN_RATIO * (HALF_HEIGHT - b);

          for(int j = 0; j < SQRT_SAMPLES; j++) {
            double a = x + INVERSE_SQRT_SAMPLES * j;
            if (SQRT_SAMPLES == 1) {
              a += 0.5;
            else {
              a += randomDoubles.random();
            }

            a = VIRTUAL_SCREEN_RATIO * (a - HALF_WIDTH);

            Vec.map(p, c, u, v, a, b);
            Vec.constructUnitVector(d, p, EYE);

            boolean hit = false;
            bestIntersection.time = Double.POSITIVE_INFINITY;
            for(int k = scene.length - 1; k >= 0; k--) {
              IObject object = scene[k];
              if (object.intersect(EYE, d, true, Double.POSITIVE_INFINITY,
                  temps, intersection)) {
                if (intersection.time < bestIntersection.time) {
                  hit = true;
                  bestIntersection.time = intersection.time;
                  Vec.assign(bestIntersection.normal, intersection.normal);
                  Vec.assign(bestIntersection.hit, intersection.hit);
                  bestIntersection.material = intersection.material;
                }
              }
            }

            if (hit) {

              if (bestIntersection.material.ambientWeight > 0) {
                for(int k = 0; k < 3; k++) {
                  pixel[k+= bestIntersection.material.ambientWeight
                      * bestIntersection.material.diffuseColor[k]
                      * RADIANCE_SCALE
                      * AMBIENT_COLOR[k];
                }
              }

              if (Vec.dot(bestIntersection.normal, d>= 0) {
                Vec.negate(bestIntersection.normal);
              }

              Vec.constructUnitVector(l, LIGHT, bestIntersection.hit);

              double nDotl = Vec.dot(l, bestIntersection.normal);
              if (nDotl <= 0) {
                continue;
              }

              double maxTime = Vec.distance(bestIntersection.hit, l);

              boolean illuminated = true;
              for(int k = scene.length - 1; k >= 0; k--) {
                IObject object = scene[k];
                if (object.intersect(bestIntersection.hit, l, false,
                    maxTime, temps, intersection)) {
                  illuminated = false;
                  break;
                }
              }

              if (illuminated) {

                if (bestIntersection.material.diffuseWeight > 0) {
                  for(int k = 0; k < 3; k++) {
                    pixel[k+= bestIntersection.material.diffuseWeight
                        * bestIntersection.material.diffuseColor[k]
                        * INVERSE_PI
                        * RADIANCE_SCALE
                        * LIGHT_COLOR[k]
                        * nDotl;
                  }
                }

                if (bestIntersection.material.specularWeight > 0) {

                  Vec.scale(r, bestIntersection.normal, 2.0 * nDotl);
                  Vec.subtract(r, l);
                  double rDotMd = Vec.dotNegative(r, d);
                  if (rDotMd > 0) {
                    for(int k = 0; k < 3; k++) {
                      pixel[k+= bestIntersection.material.specularWeight
                          * Math.pow(rDotMd,
                              bestIntersection.material.shininess)
                          * RADIANCE_SCALE
                          * LIGHT_COLOR[k]
                          * nDotl
                          * bestIntersection.material.highlightColor[k];
                    }
                  }
                }
              }
            }
          }
        }

        int value = 0;
        for(int i = 0; i < 3; i++) {
          int intensity = (int)Math.round(255
              * Math.pow(pixel[i* INVERSE_SAMPLES, INVERSE_GAMMA));
          if (intensity < 0) {
            intensity = 0;
          else if (intensity > 255) {
            intensity = 255;
          }
          value <<= 8;
          value |= intensity;
        }

        pixels[x= value;
      }

      rowCompleted(y, pixels);
    }
  }

Here is the result with 256 samples per pixel:


0 hours, 0 minutes, 3 seconds, 552 milliseconds
source

If you study the code, you’ll notice that I did not plug in the reflection term of the Phong Shading model yet.

Step 9

The ground is a checkerboard covered plane. Below is a plane containing point p with a surface normal of n. The point (x, y, z) is some other point on the plane. Note that the vector between that point and p is perpendicular to n.

Since those vectors are perpendicular, we can take advantage of dot product to derive the equation for the plane.

But, for this project, the ground is just the xz-plane. Meaning, p = (0, 0, 0) and n = (0, 1, 0). That leaves us with just y = 0. Nevertheless, we can still derive the generic equation for ray-plane intersection.

Our case is a little simpler:

Above, y = 0, simplifying it further.

Below is an aerial view of the xz-plane looking straight down the y-axis. The checkerboard consists of squares with side length s. For some point on a square, if we divide each coordinate by s and then floor it, we get a pair of integer coordinates of the upper-left corner of the square.

Notice that diagonals always have the same color. We can represent diagonals going from the lower-left to the upper-right by the following equation where t and a are integers.

Above, a is the z-intercept. Adjusting it affects where the diagonal line crosses the z-axis. If we add the components together, we get:

2t is an even number. Even plus even is even and even plus odd is odd. The parity (even/odd) of the sum is completely determined by a. This means we can determine the color of any point on the plane by dividing each coordinate component by s, flooring them and adding the resulting integer components together. The parity of the sum determines the color.

Based off the above discussion, I created the Ground:

public class Ground implements IObject {

  public double squareSize;
  public double inverseSquareSize;
  public Material material1;
  public Material material2;

  public Ground(double squareSide, Material material1, Material material2) {
    this.squareSize = squareSide;
    this.material1 = material1;
    this.material2 = material2;

    inverseSquareSize = 1.0 / squareSize;
  }

  public boolean intersect(double[] o, double[] d, boolean primaryRay,
      double maxTime, double[][] temps, Intersection intersection) {

    if (d[1== 0) {
      return false;
    }

    intersection.time = -o[1/ d[1];

    if (intersection.time >= Main.EPSILON && intersection.time <= maxTime) {

      if (!primaryRay) {
        return true;
      }

      Vec.assign(intersection.hit, o[0+ d[0* intersection.time, 0,
          o[2+ d[2* intersection.time);
      Vec.assign(intersection.normal, 010);

      long a = ((long)Math.floor(intersection.hit[0* inverseSquareSize)) 1;
      long b = ((long)Math.floor(intersection.hit[2* inverseSquareSize)) 1;
      intersection.material = (((a + b1== 0? material1 : material2;

      return true;
    }

    return false;
  }
}

I was worried that the sum of the floors might cause the long to overflow; so, I and-ed it with 1 to keep only the least-significant bit. The sum of those bits is all that matters in determining the parity of the sum.

I added 2 matte materials (materials without specular highlights and without reflections) to the Materials collection.

  public static final Material GREEN_MATTE = new Material(
    0.1,
    2.0,
    0.0,
    0.0,
    0.0,
    new double[] { 01},
    new double[] { 00},
    new double[] { 00}
  );

  public static final Material YELLOW_MATTE = new Material(
    0.1,
    2.0,
    0.0,
    0.0,
    0.0,
    new double[] { 11},
    new double[] { 00},
    new double[] { 00}
  );

Finally, I created a scene containing the ground to test it out. Below is the result with 256 samples per pixel.


0 hours, 0 minutes, 3 seconds, 244 milliseconds
source

Step 10

As mentioned above, the juggler lives under a giant dome also known as the sky.

To generate the sky, I first updated the Materials class to make it easier to define materials.

public final class Materials {

  private Materials() {
  }

  public static final Material RED_PLASTIC = createPlastic(100);

  public static final Material BLUE_MATTE = createMatte(001);
  public static final Material GREEN_MATTE = createMatte(010);
  public static final Material YELLOW_MATTE = createMatte(110);

  public static Material createPlastic(double... color) {
    return new Material(
      0.1,
      2.0,
      1.0,
      0.0,
      10.0,
      color,
      new double[] { 11},
      new double[] { 00}
    );
  }

  public static Material createMatte(double... color) {
    return new Material(
      0.1,
      2.0,
      0.0,
      0.0,
      0.0,
      color,
      new double[] { 00},
      new double[] { 00}
    );
  }
}

Next, I defined this scene in Main:

  private volatile IObject[] scene = {
    new Ground(25, Materials.GREEN_MATTE, Materials.YELLOW_MATTE),
    new Sphere(0001E6, Materials.BLUE_MATTE)
  };

Here’s the result (256 samples per pixel):


0 hours, 0 minutes, 5 seconds, 959 milliseconds
source

Step 11

I added a white matte material to Materials and then I defined a scene consisting of the ground, the sky and a sphere.

  private volatile IObject[] scene = {
    new Ground(25, Materials.WHITE_MATTE, Materials.WHITE_MATTE),
    new Sphere(0001E6, Materials.WHITE_MATTE),
    new Sphere(050050, Materials.WHITE_MATTE)
  };

Here is the result:

The shadow on the ground is exactly the same color as the dark side of the sphere. Consequentially, there is no visual indication of where one ends and the other begins. This effect is the result of the Phong Shading model. It says if the light source cannot reach a surface, then the surface is colored based off ambient light resulting in a uniform dark color. We’ll enhance the image using ambient occlusion.

Ambient occlusion estimates the amount of ambient light reaching a surface by firing random rays off of the surface originating from the hit point. If the random ray strikes an object within some threshold distance, then there is no ambient light coming from that direction. The possible random rays exist within a hemisphere around the surface normal.

Similar to RandomDoubles, I created a class called RandomRays. It pre-calculates a large number of random unit vectors.

public class RandomRays {

  private static volatile double[][] rays;

  private int index;

  public RandomRays() {
    synchronized(RandomRays.class) {
      if (rays == null) {
        rays = new double[24001][3];
        for(int i = rays.length - 1; i >= 0; i--) {
          do {
            for(int j = 0; j < 3; j++) {
              rays[i][j2.0 * Math.random() 1.0;
            }
          while(Vec.magnitude2(rays[i]) 1.0);
          Vec.normalize(rays[i]);
        }
      }
      index = (int)(Math.random() * rays.length);
    }
  }

  public void randomRay(double[] ray) {
    if (index >= rays.length) {
      index = 0;
    }
    Vec.assign(ray, rays[index++]);
  }
}

To generate the random unit vectors in a uniform way, the inner loop creates a random vector within a within a cube bounded between (−1, −1, −1)–(1, 1, 1). A unit sphere fits snuggly within the cube.

If the vector is outside of the unit sphere, it is thrown away. Otherwise, it is normalized, extending it to the surface of the unit sphere.

Next, the Material object was modified to include an ambientOcclusionPercent field. This is the percentage (0 to 1) of the ambient term for which we’ll apply the occlusion test. If the hit point is fully occluded, we’ll still give it a bit of ambient light equal to 1 minus that field.

public class Material {
  public double ambientWeight;
  public double ambientOcclusionPercent;
  public double diffuseWeight;
  public double specularWeight;
  public double reflectionWeight;  
  public double shininess;
  public double[] diffuseColor;
  public double[] highlightColor;
  public double[] reflectionColor;

  public Material(
      double ambientWeight,
      double ambientOcclusionPercent,
      double diffuseWeight,
      double specularWeight,
      double reflectionWeight,
      double shininess,
      double[] diffuseColor,
      double[] highlightColor,
      double[] reflectionColor) {
    this.ambientWeight = ambientWeight;
    this.ambientOcclusionPercent = ambientOcclusionPercent;
    this.diffuseWeight = diffuseWeight;
    this.specularWeight = specularWeight;
    this.reflectionWeight = reflectionWeight;
    this.shininess = shininess;
    this.diffuseColor = diffuseColor;
    this.highlightColor = highlightColor;
    this.reflectionColor = reflectionColor;    
  }
}

The render method will not fire multiple occlusion rays off of the hit point. Instead, for each pixel sample, it will fire only 1. That means that the number of pixel samples needs to be large to enable many occlusion tests to average out. I added a boolean constant to Main turn on/off ambient occlusion. For test images, I’ll shut it off when the number of samples is below some threshold.

The occlusion test ray needs to be within a hemisphere around the surface normal. This is accomplished by computing the dot product between the normal and the random occlusion ray. If the dot product is negative, then the occlusion ray is negated forcing the angle between the normal and the occlusion ray to be less-than 90.

            if (hit) {

              if (Vec.dot(bestIntersection.normal, d>= 0) {
                Vec.negate(bestIntersection.normal);
              }

              if (bestIntersection.material.ambientWeight > 0) {
                double ambientPercent = 1.0;

                if (AMBIENT_OCCLUSION
                    && bestIntersection.material.ambientOcclusionPercent > 0) {

                  randomRays.randomRay(ar);
                  if (Vec.dot(ar, bestIntersection.normal0) {
                    Vec.negate(ar);
                  }

                  for(int k = scene.length - 1; k >= 0; k--) {
                    IObject object = scene[k];
                    if (object.intersect(bestIntersection.hit, ar, false,
                        MAX_OCCLUSION_DISTANCE, temps, tempIntersection)) {
                      ambientPercent = 1.0
                          - bestIntersection.material.ambientOcclusionPercent;
                      break;
                    }
                  }
                }

                for(int k = 0; k < 3; k++) {
                  pixel[k+= bestIntersection.material.ambientWeight
                      * bestIntersection.material.diffuseColor[k]
                      * RADIANCE_SCALE
                      * AMBIENT_COLOR[k]
                      * ambientPercent;
                }
              }

Below you can compare before and after (1024 samples per pixel).



0 hours, 0 minutes, 31 seconds, 340 milliseconds
source

Step 12

The shadows currently have hard edges. This is because the light originates from a point-source. To alleviate that situation, we’ll use a spherical light source.

I added a constant, LIGHT_RADIUS, to Main that specifies the radius of our single light source. I modified the render method to pick a random point on the sphere. That point will act as a point-light-source for the pixel sample. In the average, we will get the effect of a spherical light source. The random point needs to be on the hemisphere around the vector defined between the hit point and the light center point. Similar to what was done for ambient occlusion, if the dot product between that vector and the random vector is negative, the random vector is negated.

              if (LIGHT_RADIUS > 0) {
                Vec.subtract(l2, bestIntersection.hit, LIGHT);
                randomRays.randomRay(l);
                if (Vec.dot(l, l20) {
                  Vec.negate(l);
                }
                Vec.ray(l, LIGHT, l, LIGHT_RADIUS);
                Vec.constructUnitVector(l, bestIntersection.hit);
              else {
                Vec.constructUnitVector(l, LIGHT, bestIntersection.hit);
              }

Below you can see the result with 1024 samples per pixel. The hazy part of the shadow is called the “penumbra.”


0 hours, 0 minutes, 33 seconds, 868 milliseconds
source

Here’s the same scene, but colorful materials are used:

Step 13

In a typical ray tracer, reflection is generated by recursively applying the Phong Shading formula (see the fourth term in the sum). When a ray strikes a translucent material, it splits into a reflected ray and a transmitted ray. Recursion is used to follow those branching rays similar to traversing a tree data structure. However, in this project, there are no translucent materials, only reflective ones. The ray never splits. This means that a for-loop can be used to effectively expand the Phong Shading formula out for a specified number of iterations.

As a first step, I dropped Pi from the Phong Shading diffuse shading term. This simplifies specifying the constants when creating materials. Next, I added a method to Materials to make metals.

  public static final Material YELLOW_METAL = createMetal(110);

  public static Material createMetal(double... color) {
    return new Material(
      0.1,
      1.0,
      0.7,
      1.0,
      0.2,
      20.0,
      color,
      color,
      color
    );
  }

Note that the diffuse, specular and reflective colors are all the same.

Next, I defined a color vector f within the render method to act as a filter. It is initialized to the radiance scaling factor (ls in the Phong Shading model) and I replaced the uses of the radiance scaling factor with the components of f. Meaning, each term of the Phong Shading model is scaled by f instead of ls.

I modified the render method to apply the Phong Shading formula up to MAX_DEPTH times.

Vec.assign(f, RADIANCE_SCALE, RADIANCE_SCALE, RADIANCE_SCALE);

reflectionLoop: for(int m = 0; m < MAX_DEPTH; m++) {

At the very bottom of the loop, the code tests if the hit material is reflective. If not, it breaks out of the loop. If the material is reflective, the code reduces the values of the filter f in accordance with the Phong Shading model; it scales them down by the reflection weight and reflection color. If all 3 components of the filter drop below a specified threshold, then it breaks out of the loop because the filter would wipe out any effect of successive iterations. Finally, the direction vector is reflected around the surface normal using the formula discussed above.

                if (bestIntersection.material.reflectionWeight > 0) {
                  boolean aboveMinColorIntensity = false;
                  for(int k = 0; k < 3; k++) {
                    f[k*= bestIntersection.material.reflectionWeight
                        * bestIntersection.material.reflectionColor[k];
                    if (f[k>= MIN_COLOR_INTENSITY) {
                      aboveMinColorIntensity = true;
                    }
                  }
                  if (!aboveMinColorIntensity) {
                    break reflectionLoop;
                  }

                  Vec.ray(d, bestIntersection.normal,
                      -* Vec.dot(bestIntersection.normal, d));
                else {
                  break reflectionLoop;
                }

This reflection loop is nested within the sampling loop. Meaning, the reflected rays contribute to the same pixel sample.

Here are the results with 256 samples per pixel. Note that the sphere does not reflect the sky because the reflection color, yellow, does not contain a blue component. However, yellow does contain a green component; hence, it is able to reflect both colors of the checkered ground.


0 hours, 0 minutes, 12 seconds, 573 milliseconds
source

To create a mirrored sphere, the ambient and the diffuse weights are set to 0. The reflective weight is set to 1 and the reflective color is set to white. The specular weight should also be set to 0, but I set it to 0.5 to create a specular highlight on the mirrored sphere. In the real world, a specular highlight is a blurry reflection of a light source. On a perfectly mirrored sphere, instead of a specular highlight, you would see a light bulb or the sun. In fact, the specular term of the Phong Shading model is essentially a hack. It can be replaced by modifying the reflective term to handle glossy reflections. A glossy reflection of a glowing object would appear as a specular highlight. But, for this project, I am not going to implement glossy reflections or glowing objects. I’ll stick with the hack and just add a specular highlight to the mirrored sphere.

Below is a mirror sphere computed with 1024 samples per pixel.


0 hours, 0 minutes, 50 seconds, 14 milliseconds
source

Here is a red metallic sphere, a mirrored sphere and a yellow metallic sphere rendered at 256 samples per pixel. If you carefully compare the image below against the 2 previous images above, you’ll notice that there is something slightly wrong with the 2 previous images. Above, the reflective spheres do not reflect their own shadows. This is due to a bug that was fixed before generating the image below. Consider this document as a progress log. I am writing it in one pass as I develop the code and I failed to notice the problem when I attached the previous images.


0 hours, 0 minutes, 20 seconds, 15 milliseconds
source

Now, to answer the all important question, what happens if we put a few colored metallic spheres, a light source and the virtual camera within a mirrored sphere?


0 hours, 1 minute, 17 seconds, 142 milliseconds
source

Step 14

To simulate motion blur, the scene must be updated between each pixel sample in accordance with the motion occurring in that animation frame. This animation will consist of 30 frames indexed from 0 to 29 (inclusive), but time will go from 0 to 30 time units where each time unit is 1/30th of a second. Some arbitrary frame, x, represents the time period x (inclusive) to x + 1 (exclusive). In the timeline below, the vertical lines are points in time. Below the timeline are the frame indices.

A unit of time needs to be sampled similarly to how a pixel is sampled. The time unit will be divided by the number of samples and within each division, a random sample will be selected.

Each pixel sample will have an associated time sample from 0 (inclusive) to 1 (exclusive). However, pixel sampling is occurring within each pixel by scanning from left-to-right, top-to-bottom. I want to make sure that there is no correlation between the pixel sample and the time sample; that is, I do not want the upper-left corner of a pixel to always have a small time sample and the lower-right corner of a pixel to have a high time sample. To that end, I created a class called RandomTimes. After it divides up the unit time and selects random times within each division, it shuffles the order of the samples.

import java.util.*;

public class RandomTimes {

  private static volatile double[][] times;

  private int index;

  public RandomTimes() {
    synchronized(RandomTimes.class) {

      ArrayList<Double> doubles = new ArrayList<Double>();

      if (times == null) {
        times = new double[503][Main.SAMPLES];
        for(int i = 0; i < times.length; i++) {
          doubles.clear();
          for(int j = 0; j < Main.SAMPLES; j++) {
            doubles.add(Main.INVERSE_SAMPLES * (j + Math.random()));
          }
          for(int j = 0; j < Main.SAMPLES; j++) {
            times[i][j= doubles.remove((int)(Math.random() * doubles.size()));
          }
        }
      }
      index = (int)(Math.random() * times.length);
    }
  }

  public double[] random() {
    if (index >= times.length) {
      index = 0;
    }
    return times[index++];
  }
}

Since each thread will update the scene between casting rays, each thread needs its own copy of the scene. To generate it, I added a createScene method. To update the scene, I introduced an updateScene method that takes as parameters: the frame index, the time sample within the frame (0 to 1), and the scene to be updated. As a test, I created a scene containing the ground, the sky and a red metallic sphere. The updateScene method moves the red metallic sphere linearly based on the time sample.

  private IObject[] createScene() {
    return new IObject[] {
      new Ground(25, Materials.YELLOW_MATTE, Materials.GREEN_MATTE),
      new Sphere(0001E6, Materials.BLUE_MATTE),
      new Sphere(-2050050, Materials.RED_METAL),
    };
  }

  private void updateScene(int frame, double t, IObject[] scene) {
    Sphere sphere = (Sphere)scene[2];
    sphere.center[0= -20 40 * t;
  }

The nested loops within render now looks like this:

    while(true) {

      int y = getNextRowIndex();
      if (y >= HEIGHT) {
        return;
      }

      for(int x = 0; x < WIDTH; x++) {

        Vec.assign(pixel, 000);

        double[] times = randomTimes.random();
        int timesIndex = 0;

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

          double b = y + INVERSE_SQRT_SAMPLES * i;
          if (SQRT_SAMPLES == 1) {
            b += 0.5;
          else {
            b += randomDoubles.random();
          }

          b = VIRTUAL_SCREEN_RATIO * (HALF_HEIGHT - b);

          for(int j = 0; j < SQRT_SAMPLES; j++) {
            double a = x + INVERSE_SQRT_SAMPLES * j;
            if (SQRT_SAMPLES == 1) {
              a += 0.5;
            else {
              a += randomDoubles.random();
            }

            a = VIRTUAL_SCREEN_RATIO * (a - HALF_WIDTH);

            Vec.map(p, c, u, v, a, b);
            Vec.constructUnitVector(d, p, EYE);

            Vec.assign(o, EYE);
            Vec.assign(f, RADIANCE_SCALE, RADIANCE_SCALE, RADIANCE_SCALE);

            updateScene(frameIndex, times[timesIndex++], scene);

            reflectionLoop: for(int m = 0; m < MAX_DEPTH; m++) {

Above, for each pixel, the code obtains a new set of time samples. Just before it enters the reflection loop, it updates the scene with a different time sample. The results look like this (256 samples/pixel):


0 hours, 0 minutes, 13 seconds, 401 milliseconds
source

Step 15

Typically, when a person learns to juggle, the first thing he or she picks up is a figure-eight pattern known as the “cascade.” But, the Amiga Juggler tosses balls in a circular pattern called the “shower.” The pattern consists of 2 arcs. It tosses balls on a high arc from his left hand to his right hand. It returns the balls in a low arc directly below the high arc.

When a ball is tossed from hand-to-hand it follows a ballistic trajectory.

The x and the y components of the trajectory are independent. The ball moves linearly in the x-direction starting from an initial value of x0 and travels at a constant rate v0. In this way, the x component is governed by the ray equation. We know the amount of time it takes for the ball to traverse the arc and that provides a way to compute the velocity in the x-direction.

The maximum velocity in the y-direction, v, is a function of the height of the arc, h, and the acceleration due to gravity, g. Since all the potential energy stored in the ball at the apex is converted into kinetic energy when it reaches the bottom of the arc, we can express the relationship between v, h and g by equating the formulae for both types of energy.

Movement in the y-direction is governed by the parabolic formula below where y0 is the initial height of the ball, v0 is the initial velocity in the y-direction, g is the acceleration due to gravity (a positive constant) and t is time.

Due to the symmetry of the arc, the starting and the ending heights are the same. Also, the magnitudes of the starting and ending velocities in the y-direction are the same, but opposite in sign.

The derivation below shows that the height of the arc and the time it takes the ball to traverse the arc determines the maximum velocity of the ball in the negative y-direction. Also, from above, the acceleration due to gravity is computable from the maximum velocity and the arc height.

To test out these formulas, I modified the launch method to generate 30 animation frames.

  public void launch() throws Throwable {

    image = new BufferedImage(WIDTH, HEIGHT, BufferedImage.TYPE_INT_RGB);
    if (RENDER_IN_WINDOW) {
      renderFrame = new RenderFrame(image);
    }
    
    int processors = Runtime.getRuntime().availableProcessors();    

    for(int frameIndex = 0; frameIndex < 30; frameIndex++) {

      if (RENDER_IN_WINDOW) {
        renderFrame.setTitle("Amiga Juggler [" + frameIndex + "]");
      }

      resetRowIndex();
      updateRunningCount(processors);
      for(int i = 0; i < processors; i++) {
        final int _frameIndex = frameIndex;
        new Thread(Integer.toString(i)) {
          @Override
          public void run() {
            render(_frameIndex);
            updateRunningCount(-1);
          }
        }.start();
      }

      synchronized(this) {
        while(runningCount != 0) {
          wait();
        }
      }

      saveImage("balls" ((frameIndex < 10"0" ""+ frameIndex);
    }
  }

Next, I added a series of constants to control the arcs. Essentially, you can specify the coordinates of the hands and the height of the upper arc. The rest of the constants are derived from these values and the fact that it takes 30 frames to traverse the lower arc and 60 frames to traverse the upper arc.

  public static final double JUGGLE_X0 = -75;
  public static final double JUGGLE_X1 = 75;
  public static final double JUGGLE_Y0 = 10;
  public static final double JUGGLE_H_Y = 90;

  public static final double JUGGLE_H_VX = (JUGGLE_X0 - JUGGLE_X160.0;
  public static final double JUGGLE_L_VX = (JUGGLE_X1 - JUGGLE_X030.0;

  public static final double JUGGLE_H_H = JUGGLE_H_Y - JUGGLE_Y0;
  public static final double JUGGLE_H_VY = 4.0 * JUGGLE_H_H / 60.0;
  public static final double JUGGLE_G = JUGGLE_H_VY * JUGGLE_H_VY
      (2.0 * JUGGLE_H_H);

  public static final double JUGGLE_L_VY = 0.5 * JUGGLE_G * 30.0;

Then it was just a matter of modifying the createScene and updateScene methods to incorporate the parabolic arc formulae.

  private IObject[] createScene() {

    return new IObject[] {
      new Ground(25, Materials.YELLOW_MATTE, Materials.GREEN_MATTE)// 0
      new Sphere(0001E6, Materials.BLUE_MATTE),                // 1
      
      new Sphere(00010, Materials.RED_PLASTIC)// 2
      new Sphere(00010, Materials.RED_PLASTIC)// 3
      new Sphere(00010, Materials.RED_PLASTIC)// 4
    };
  }

  private void updateScene(int frameIndex, double t, IObject[] scene) {
    double T = (30 + frameIndex + t);
    Sphere sphere = (Sphere)scene[3];
    sphere.center[0= JUGGLE_X1 + JUGGLE_H_VX * T;
    sphere.center[1= JUGGLE_Y0 + (JUGGLE_H_VY - 0.5 * JUGGLE_G * T* T;

    T = frameIndex + t;
    sphere = (Sphere)scene[4];
    sphere.center[0= JUGGLE_X1 + JUGGLE_H_VX * T;
    sphere.center[1= JUGGLE_Y0 + (JUGGLE_H_VY - 0.5 * JUGGLE_G * T* T;

    sphere = (Sphere)scene[2];
    sphere.center[0= JUGGLE_X0 + JUGGLE_L_VX * T;
    sphere.center[1= JUGGLE_Y0 + (JUGGLE_L_VY - 0.5 * JUGGLE_G * T* T;
  }

The program outputs 30 gifs at 320×180 resolution with 256 samples/pixel. I used UnFREEz to combine the images into an animated gif. Below, the background seems to flicker a bit because the animated gif format is crummy. That effect does not exist in the individual images.


0 hours, 1 minute, 14 seconds, 209 milliseconds
source

Step 16

The arms and legs of the juggler consist of spheres packed into lines. As far as I can tell, the juggler stands with his right foot forward and his left foot back. As he juggles, the hips move up and down. In one animation frame, the left leg looks almost straight; where as the right knee is slightly bent.

   

Each appendage consists of 2 line segments connected by a joint. To make appendages easy to control, I established the following model: The lengths of the line segments are A and B. The appendage exists in a plane with normal w. The end points of the appendage are p and q. The normalized difference between those points establishes a vector v. The cross product of v and w generates u, completing an orthonormal basis (ONB) for the appendage. I need to find the location of the joint, (x, y), which appears in the direction of positive u. The point (x, y) is in the coordinate system defined by the ONB where p is the origin and z = 0.

 

You can redraw the appendage as 2 right-triangles that share a common height.

From there it’s just a matter of applying the Pythagorean Theorem twice and solving for x and y.

If the value within the square-root is negative, then it is not possible to construct an appendage with the specified parameters.

To test the concept, I created a scene containing an appendage. The appendage consists of 2 segments, each composed of 8 spheres, plus an extra sphere for the joint. In total, that’s 17 spheres.

  private IObject[] createScene() {

    IObject[] scene = new IObject[8];
    scene[0new Ground(25, Materials.YELLOW_MATTE, Materials.GREEN_MATTE);
    scene[1new Sphere(0001E6, Materials.BLUE_MATTE);
    for(int i = 2; i < scene.length; i++) {
      scene[inew Sphere(00010, Materials.RED_PLASTIC);
    }

    return scene;
  }

Below is a method I added to update an appendage based on the parameters discussed above. sceneIndex is the offset into scene where the appendage is stored. The parameters p, q, w, A and B control the appendage as described above. countA and countB are the number of spheres in each segment. Neither count includes the joint; the total number of spheres in the appendage is countA + 1 + countB. The spheres must be in the following order within scene: segment A spheres from p to j, the joint sphere at j, and segment B spheres from q to j. The temps parameter is a set of temporary vectors used for intermediate calculations.

  private void updateAppendage(
      IObject[] scene, int sceneIndex,
      double[] p, double[] q, double[] w, double A, double B,
      int countA, int countB, double[][] temps) {

    double[] U = temps[0];
    double[] V = temps[1];
    double[] W = temps[2];
    double[] j = temps[3];
    double[] d = temps[4];

    Vec.subtract(V, q, p);
    double D = Vec.magnitude(V);
    double inverseD = 1.0 / D;
    Vec.scale(V, inverseD);

    Vec.normalize(W, w);
    Vec.cross(U, V, W);

    double A2 = A * A;

    double y = 0.5 * inverseD * (A2 - B * B + D * D);
    double x = Math.sqrt(A2 - y * y);

    Vec.map(j, p, U, V, x, y);

    Vec.subtract(d, j, p);
    Vec.scale(d, 1.0 / countA);
    for(int i = 0; i <= countA; i++) {
      double[] center = ((Sphere)scene[sceneIndex + i]).center;
      Vec.ray(center, p, d, i);
    }

    Vec.subtract(d, j, q);
    Vec.scale(d, 1.0 / countB);
    for(int i = 0; i < countB; i++) {
      double[] center = ((Sphere)scene[countA + + sceneIndex + i]).center;
      Vec.ray(center, q, d, i);      
    }
  }

Finally, I modified updateScene to move the appendage. The appendage exists within the xy-plane, parallel to the virtual screen. The code anchors point p and it moves point q in a circular motion.

  private void updateScene(
      int frameIndex, double t, IObject[] scene, double[][] temps) {

    double[] w = temps[5];
    w[00;
    w[10;
    w[21;

    double[] p = temps[6];
    p[0= -75;
    p[150;
    p[20;

    double angle = 2.0 * Math.PI * (frameIndex + t30.0;
    double[] q = temps[7];
    q[075 * Math.cos(angle);
    q[150 * Math.sin(angle);
    q[20;

    updateAppendage(scene, 2, p, q, w, 808088, temps);
  }

Here’s the result (320×180, 256 samples/pixel):


0 hours, 3 minutes, 36 seconds, 200 milliseconds
source

Step 17

Here is a frame from the Amiga Juggler animation. I added black bars to the left and right to force the image into the aspect ratio that I am interested in and then I scaled it to the size of my other test images.

The horizon line is slightly above center, which means that the camera is tilted somewhat downwards.

To help me match the original animation, I modified the program to load and blend this frame with my own. I introduced the following method to load the image into a double[][][] array. Note that it applies gamma correction.

  private void loadImage(String namethrows Throwable {
    BufferedImage blendImage = ImageIO.read(new File(name));
    blendPixels = new double[blendImage.getHeight()][blendImage.getWidth()][3];
    int[] row = new int[blendImage.getWidth()];
    for(int i = 0; i < blendImage.getHeight(); i++) {
      blendImage.getRGB(0, i, blendImage.getWidth()1, row, 0,
          blendImage.getWidth());
      for(int j = 0; j < blendImage.getWidth(); j++) {
        for(int k = 2; k >= 0; k--) {
          int intensity = row[j0xFF;
          row[j>>= 8;
          blendPixels[i][j][k= Math.pow(intensity / 255.0, GAMMA);
        }
      }
    }
  }

Blending occurs within the render method using a weighted average. Currently, it’s equally weighted.

pixel[i0.5 * pixel[i0.5 * blendPixels[y][x][i];

The launch method loads the image and then it attaches a key listener to the frame. By pressing the up or down arrow-keys, I can control the y-coordinate of the look point. This enabled me to find the camera tilt. I’ll reuse this framework in successive steps to further match the original animation.

  public void launch() throws Throwable {
    
    loadImage("frame00.png");

    image = new BufferedImage(WIDTH, HEIGHT, BufferedImage.TYPE_INT_RGB);
    if (RENDER_IN_WINDOW) {
      renderFrame = new RenderFrame(image);
      addKeyListener();
    }
    
    int processors = Runtime.getRuntime().availableProcessors();    

    while(true) {

      resetRowIndex();
      updateRunningCount(processors);
      for(int i = 0; i < processors; i++) {
        final int _frameIndex = 0;
        new Thread(Integer.toString(i)) {
          @Override
          public void run() {
            render(_frameIndex);
            updateRunningCount(-1);
          }
        }.start();
      }

      synchronized(this) {
        while(runningCount != 0) {
          wait();
        }
      }

      synchronized(this) {
        while(!keyPressed) {          
          wait();
        }

        keyPressed = false;
        switch(keyCode) {
          case KeyEvent.VK_UP:
            LOOK[1+= 1;
            break;
          case KeyEvent.VK_DOWN:
            LOOK[1-= 1;
            break;
        }

        System.out.format("%f%n", LOOK[1]);
      }      
    }
  }

  private void addKeyListener() {
    renderFrame.addKeyListener(new KeyAdapter() {
      @Override
      public void keyPressed(KeyEvent e) {
        synchronized(Main.this) {
          keyCode = e.getKeyCode();
          keyPressed = true;
          Main.this.notifyAll();
        }
      }
    });
  }

source

Step 18

It dawned on me while I was trying to work out the camera angle that the Amiga 1000 does not display square pixels. To compensate, I scaled the animation frame horizontally by 86%.

Here is my attempt at reproducing the same view. It’s not a perfect match, but it’s close enough.


source

Step 19

Using the shadows, the square tiles and the spheres in the original animation frame as a guide, I worked out key coordinates in space that will help me align everything else. Below is my attempt to orient the mirrored spheres. Now all I need is the juggler.


0 hours, 1 minute, 10 seconds, 138 milliseconds
source

Step 20

Next, I sized and positioned the head, neck and torso spheres (64 samples/pixel).


0 hours, 0 minutes, 40 seconds, 529 milliseconds
source

It’s unclear if the juggler in the original animation has a large, white button on his chest or if that is just a large specular highlight as seen in the above animation.

Step 21

The torso definitely bobs up and down during juggling. I measured the range of the bob and I used cosine to simulate the effect. Below is the result with 64 samples/pixel.


0 hours, 0 minutes, 42 seconds, 385 milliseconds
source

Step 22

The head also sways side-to-side. To achieve that, I constructed an orthonormal basis centered at the crotch point between the hips. In the image below, the unit vector v follows a circular path. Using v and the ray equation is sufficient to position the torso, the neck and the head. The unit vector u will be used with v for the shoulders, the hips and the appendages.

Here is the result with 64 samples per pixel.


0 hours, 0 minutes, 45 seconds, 459 milliseconds
source

Step 23

Using the ONB created in the prior step, I attached the legs. Below you can see the result with 64 samples/pixel. The knees bend independently because one foot is forward and the other is back and because the hips are tilting from side-to-side. The thigh is 80% shorter than the calf.


0 hours, 2 minutes, 15 seconds, 349 milliseconds
source

Step 24

Again using the crotch-centered coordinate system, I attached the arms to the shoulder points. The normal vector of each arm plane is based on the angle between the elbow point and the shoulder point. The hand points rotate around an arbitrary point in space such that the angle is based on cosine of time.

Below is the result with 16 samples/pixel.


0 hours, 0 minutes, 54 seconds, 850 milliseconds
source

Step 25

To help me match colors sampled from frames of the original animation, I added 2 new methods to the Materials class.

  public static double[] createColor(int hexColor, double scale) {
    double[] color = new double[3];
    createColor(color, hexColor, scale);
    return color;
  }

  public static void createColor(
      double[] color, int hexColor, double scale) {
    for(int i = 2; i >= 0; i--) {
      int intensity = hexColor & 0xFF;
      hexColor >>= 8;
      color[i= Math.pow(scale * intensity / 255.0, Main.GAMMA);
    }
  }

They accept a color in the form of an RGB hex value, the format readily available from photo-editing software. Gamma correction is applied to convert the color into the form necessary for the ray tracer. The second parameter provides a convenient way to scale the result.

In the original animation, the yellow and green squares closest to the viewer are pure, max yellow and green respectively. On the ground, within the juggler’s shadow, the squares are reduced to 40% of their original color intensities. This means that the ambient weight in the Phong Shading Model is set to 0.4 for the ground material. For the spheres, it is closer to 0.15.

It was relatively easy to measure colors and to derive materials from them to color the ground and the spheres; however, the sky proved to be a bigger challenge. Close to the horizon the sky color is washed out with a high saturation, but high above, the sky is dark blue as revealed by the reflection in the mirrored spheres. That’s an effect that cannot be produced by diffuse and specular lighting with a single point light source.

To overcome the problem, I added a method called update to the Material class. The method is an empty stub.

  public void update(double[] hitPoint) {
  }

I modified the render method to call update with the ray-intersection hit point just prior to using the material in the Phong Shading Model. Then, I introduced a class that extends Material and overrides update to color the sky based on elevation.

public class SkyMaterial extends Material {

  private double[] minColor = Materials.createColor(0xBDBDFF1.0);
  private double[] maxColor = Materials.createColor(0x2223F61.0);

  public SkyMaterial() {
    super(100000new double[3]new double[3]new double[3]);
  }

  @Override
  public void update(double[] hitPoint) {
    Vec.interpolate(diffuseColor, minColor, maxColor, hitPoint[11E-6);
  }  
}

Notice that the sky material uses only ambient lighting, ignoring the other terms of the Phong Shading Model. The interpolate method used within update is essentially the ray equation. It performs linear interpolation between 2 specified points.

  // a = b + (c - b) * s   [s = 0 is b, s = 1 is c]
  public static void interpolate(double[] a, double[] b, double[] c, double s) {
    for(int i = 0; i < 3; i++) {
      a[i= b[i(c[i- b[i]) * s;
    }
  }

  // a = a + (b - a) * s   [s = 0 is a, s = 1 is b]
  public static void interpolate(double[] a, double[] b, double s) {
    for(int i = 0; i < 3; i++) {
      a[i= a[i(b[i- a[i]) * s;
    }
  }

Here is the result with 256 samples per pixel.


0 hours, 1 minute, 41 seconds, 322 milliseconds
source

On a side note, in the original animation, I found a frame where one of the balls casts a shadow over the juggler’s chest. The “white button” vanishes in the shadow, which demonstrates that it is indeed just a specular highlight.

Step 26

Next, I worked out the position of the eyes and the hair.

   

Below is the result with 49 samples/pixel.


0 hours, 2 minutes, 26 seconds, 609 milliseconds
source

Here’s a frame rendered with 1024 samples/pixel.

Step 27

Now, it’s time to render. But, before we do, let’s take a look at the final versions of createScene and updateScene.

  private IObject[] createScene() {

    // juggling balls 2 -- 4
    // body (hips to chest) 5 -- 12
    // head 13
    // neck 14
    // left leg (15 -- 31)
    // right leg (32 -- 48)
    // left arm (49 -- 65)
    // right arm (66 -- 82)
    // left eye 83
    // right eye 34
    // hair 85

    IObject[] scene = new IObject[86];
    scene[0new Ground(Materials.YELLOW_MATTE, Materials.GREEN_MATTE);
    scene[1new Sphere(0001E6, new SkyMaterial());
    for(int i = 2; i <= 4; i++) {
      scene[inew Sphere(1100014, Materials.MIRROR);
    }
    for(int i = 5; i <= 12; i++) {
      double percent = (i - 57.0;
      scene[inew Sphere(15185 32 * percent, -15116 * percent,
          Materials.TORSO);
    }
    scene[13new Sphere(151155, -15114, Materials.SKIN);
    scene[14new Sphere(151140, -1515, Materials.SKIN);

    for(int i = 15; i <= 22; i++) {
      for(int j = 0; j < 4; j++) {
        scene[i + 17 * jnew Sphere(0002.5 2.5 (i - 157.0,
            Materials.SKIN);
      }
    }
    for(int i = 23; i <= 31; i++) {
      for(int j = 0; j < 4; j++) {
        scene[i + 17 * jnew Sphere(0005, Materials.SKIN);
      }
    }
    scene[83new Sphere(142154, -1444, Materials.EYE);
    scene[84new Sphere(142154, -1444, Materials.EYE);
    scene[85new Sphere(152156, -15114, Materials.HAIR);

    return scene;
  }

  private void updateScene(
      int frameIndex, double t, IObject[] scene, double[][] temps) {

    // mirrored juggling balls (2 -- 4)
    double T = (30 + frameIndex + t);
    Sphere sphere = (Sphere)scene[3];
    sphere.center[2= JUGGLE_X1 + JUGGLE_H_VX * T;
    sphere.center[1= JUGGLE_Y0 + (JUGGLE_H_VY - 0.5 * JUGGLE_G * T* T;

    T = frameIndex + t;
    sphere = (Sphere)scene[4];
    sphere.center[2= JUGGLE_X1 + JUGGLE_H_VX * T;
    sphere.center[1= JUGGLE_Y0 + (JUGGLE_H_VY - 0.5 * JUGGLE_G * T* T;

    sphere = (Sphere)scene[2];
    sphere.center[2= JUGGLE_X0 + JUGGLE_L_VX * T;
    sphere.center[1= JUGGLE_Y0 + (JUGGLE_L_VY - 0.5 * JUGGLE_G * T* T;

    // body (hips to chest) 5 -- 12
    double angle = HIPS_ANGLE_MULTIPLIER * T;
    double oscillation = 0.5 (1.0 + Math.cos(angle));

    double[] o = temps[5];
    o[0151;
    o[1= HIPS_MIN_Y + (HIPS_MAX_Y - HIPS_MIN_Y* oscillation;
    o[2= -151;
    
    double[] u = temps[6];
    double[] v = temps[7];
    double[] w = temps[8];

    v[00;
    v[170;
    v[2(HIPS_MIN_Y - HIPS_MAX_Y* Math.sin(angle);
    Vec.normalize(v);

    u[00;
    u[1= v[2];
    u[2= -v[1];

    w[01;
    w[10;
    w[20;

    for(int i = 5; i <= 12; i++) {
      double percent = (i - 57.0;
      sphere = (Sphere)scene[i];
      Vec.ray(sphere.center, o, v, 32 * percent);
    }
    sphere = (Sphere)scene[13];
    Vec.ray(sphere.center, o, v, 70);
    sphere = (Sphere)scene[14];
    Vec.ray(sphere.center, o, v, 55);

    // left leg (15 -- 32)
    double[] p = temps[9];
    p[0159;
    p[12.5;
    p[2= -133;
    double[] q = temps[10];
    Vec.mapYZ(q, o, v, u, -9, -16);
    updateAppendage(scene, 15, p, q, u, 42.5834.0788, temps);

    // right leg (32 -- 48)
    p[0139;
    p[12.5;
    p[2= -164;
    Vec.mapYZ(q, o, v, u, -916);
    updateAppendage(scene, 32, p, q, u, 42.5834.0788, temps);

    // left arm (49 -- 65)
    double[] n = temps[11];
    double armAngle = -0.35 * oscillation;
    p[069 41 * Math.cos(armAngle);
    p[160 41 * Math.sin(armAngle);
    p[2= -108;
    Vec.mapYZ(q, o, v, u, 45, -19);    
    Vec.mapYZ(n, o, v, u, 45.41217, -19.91111);
    Vec.subtract(n, q);
    updateAppendage(scene, 49, p, q, n, 44.29446.09888, temps);

    // right arm (66 -- 82)
    p[2= -182;
    Vec.mapYZ(q, o, v, u, 4519);
    Vec.mapYZ(n, o, v, u, 45.4121719.91111);
    Vec.subtract(n, q, n);
    updateAppendage(scene, 66, p, q, n, 44.29446.09888, temps);

    // left eye (83)
    sphere = (Sphere)scene[83];
    Vec.mapYZ(sphere.center, o, v, u, 69, -7);
    sphere.center[0142;

    // right eye (84)
    sphere = (Sphere)scene[84];
    Vec.mapYZ(sphere.center, o, v, u, 697);
    sphere.center[0142

    // hair (85)
    sphere = (Sphere)scene[85];
    Vec.ray(sphere.center, o, v, 71);
    sphere.center[0152;
  }

As you can see, the scene consists of 85 spheres and the ground plane. Any initial sphere positions specified in createScene are actually overridden by the formulae in updateScene. The createScene method only determines the sizes and the materials of the spheres.

I rendered the frames at 640×360 and at 1920×1080 with 1024 samples/pixel. At the larger (HD) resolution, it took approximately 1 hour to render each frame. Ironically, according to Eric Graham’s documentation, each frame of the original animation also took about an hour to render on the Amiga 1000. But, the original animation was rendered at 320×200 and the renderer only used 1 sample/pixel. That’s 64,000 rays versus 2,123,366,400 rays for each image. It’s also unlikely that it applied ambient occlusion or that it used double-precision floating-point operations. Being conservative, it would have taken about 114 years to generate the complete HD animation on the Amiga 1000. If you launched the program back in 1985 when the Amiga 1000 was introduced, on 5 of those machines, they would have collectively completed the job in 2008.

I merged the frames together using VirtualDub. At 30 frames/sec, there is only 1 second of animation. I wanted it to loop 30 times, but I could not get VirtualDub to import the images more than once. So, I simply wrote a program to duplicate the files and then I imported all 900.

package jugglerfilenames;

import java.io.*;

public class Main {

  private static final String path
      "C:/NetBeansProjects/JugglerFileNames/files/";
  private static byte[] buffer = new byte[1024 1024];

  public static void main(String[] argsthrows Throwable {
    for(int i = 0; i < 30; i++) {
      String inputFile = String.format("%sjuggler%02d.png", path, i);
      for(int j = 0; j < 30; j++) {
        String outputFile = String.format(
            "%sjuggler%03d.png", path, j * 30 + i);
        copy(inputFile, outputFile);
      }
    }
  }

  private static void copy(String file1, String file2throws Throwable {
    FileInputStream in = new FileInputStream(file1);
    FileOutputStream out = new FileOutputStream(file2);

    while(true) {
      int bytes = in.read(buffer);
      if (bytes < 0) {
        break;
      }
      out.write(buffer, 0, bytes);
    }

    in.close();
    out.close();
  }

}

The original animation was accompanied by a clank sound every time the mirrored spheres contacted the robot’s arms. I decided to add some royalty free music reminiscent of the demoscene. I downloaded "Double O" by Kevin MacLeod (incompetech.com), I edited it with Audacity and then I imported it into VirtualDub.

VirtualDub produced an uncompressed AVI file. I used SUPER © to compress it into an MP4 before uploading it to YouTube.

Here is the final version of the source used to render the HD frames. Below are results.


3 hours, 19 minutes, 28 seconds, 616 milliseconds
http://www.youtube.com/watch?v=hI2J27wNAB0&feature=channel


29 hours, 55 minutes, 17 seconds, 544 milliseconds
http://www.youtube.com/watch?v=Vdh-fPX1gSk

2010.05.01