Added compilerbouw & graphics ass 7: raytracing 2

.gitmodules

+[submodule "compilerbouw"]
+	path = compilerbouw
+	url = ssh://sander@vo20.nl/git/compilerbouw.git

compilerbouw @ 28ecce81

+Subproject commit 28ecce816a2ffceb9eb23e05b39c4b403bd57539

graphics/ass7/Makefile

+CC=gcc
+
+WARNING_FLAGS=-Wall -Wextra -Werror-implicit-function-declaration -Wshadow -Wstrict-prototypes -pedantic-errors
+CFLAGS=-g -O3 -std=c99 $(WARNING_FLAGS)
+LDFLAGS=-g -lGL -lglut -lGLU
+
+.c.o:
+	$(CC) -c $(CFLAGS) $<
+
+all: main
+main: bbox.o bvh.o intersection.o main.o perlin.o plymodel.o quat.o rply.o scene.o shaders.o
+	$(CC) $(LDFLAGS) -o main bbox.o bvh.o intersection.o main.o perlin.o plymodel.o quat.o rply.o scene.o shaders.o
+
+clean:
+	rm -f *.o
+	rm -f main
+
+shaders.o      : constants.h intersection.h perlin.h quat.h scene.h shaders.h v3math.h shaders.c
+shaders.o      : types.h shaders.h
+perlin.o       : perlin.h
+main.o         : bvh.h constants.h intersection.h perlin.h plymodel.h scene.h shaders.h v3math.h main.c
+perlin.o       : perlin.h perlin.c
+rply.o         : rply.h rply.c
+v3math.o       : v3math.h
+quat.o         : v3math.h quat.h
+rply.o         : rply.h
+quat.o         : constants.h quat.h quat.c
+scene.o        : bvh.h plymodel.h scene.h scene.c
+scene.o        : types.h scene.h
+types.o        : v3math.h types.h
+bvh.o          : bvh.h scene.h v3math.h bvh.c
+plymodel.o     : plymodel.h rply.h v3math.h plymodel.c
+plymodel.o     : types.h plymodel.h
+bvh.o          : bbox.h bvh.h
+constants.o    : types.h constants.h
+intersection.o : types.h intersection.h
+intersection.o : bvh.h constants.h intersection.h scene.h v3math.h intersection.c
+bbox.o         : bbox.h constants.h v3math.h bbox.c
+bbox.o         : types.h bbox.h

graphics/ass7/bbox.c

+#include "constants.h"
+#include "v3math.h"
+#include "bbox.h"
+#include "intersection.h"
+
+// Create a bounding box (initialized to be empty, i.e. min > max)
+boundingbox
+bbox_create(void)
+{
+    boundingbox bbox;
+
+    bbox.min = v3_create(C_INFINITY, C_INFINITY, C_INFINITY);
+    bbox.max = v3_create(-C_INFINITY, -C_INFINITY, -C_INFINITY);
+
+    return bbox;
+}
+
+// Update the given bounding box to include the point p
+void
+bbox_update(boundingbox* bbox, vec3 p)
+{
+    if (p.x < bbox->min.x) bbox->min.x = p.x;
+    if (p.y < bbox->min.y) bbox->min.y = p.y;
+    if (p.z < bbox->min.z) bbox->min.z = p.z;
+
+    if (p.x > bbox->max.x) bbox->max.x = p.x;
+    if (p.y > bbox->max.y) bbox->max.y = p.y;
+    if (p.z > bbox->max.z) bbox->max.z = p.z;
+}
+
+// Return the bounding box resulting from the combination of the
+// two given bounding boxes
+boundingbox
+bbox_combine(boundingbox bbox1, boundingbox bbox2)
+{
+    boundingbox bbox = bbox1;
+
+    if (bbox2.min.x < bbox.min.x) bbox.min.x = bbox2.min.x;
+    if (bbox2.min.y < bbox.min.y) bbox.min.y = bbox2.min.y;
+    if (bbox2.min.z < bbox.min.z) bbox.min.z = bbox2.min.z;
+
+    if (bbox2.max.x > bbox.max.x) bbox.max.x = bbox2.max.x;
+    if (bbox2.max.y > bbox.max.y) bbox.max.y = bbox2.max.y;
+    if (bbox2.max.z > bbox.max.z) bbox.max.z = bbox2.max.z;
+
+    return bbox;
+}
+
+
+// Intersect a ray with given origin and direction, with the given
+// bounding box.
+// Returns 1 (and sets t_min and t_max) when an intersection is found,
+// returns 0 otherwise
+int
+bbox_intersect(float* t_min, float* t_max, boundingbox bbox,
+    vec3 origin, vec3 direction, float t0, float t1)
+{
+    // This code makes handy use of IEEE floating-point behaviour, such
+    // as division by zero. See
+    // 1) Shirley et.al. Section 10.9.1
+    // 2) Williams et.al., "An Efficient and Robust Ray-Box Intersection Algorithm"
+
+    float   tmin, tmax;
+    float   tymin, tymax, tzmin, tzmax;
+    float   inv;
+
+    num_bboxes_tested++;
+
+    inv = 1.0 / direction.x;
+    if (inv >= 0.0)
+    {
+        tmin = (bbox.min.x - origin.x) * inv;
+        tmax = (bbox.max.x - origin.x) * inv;
+    }
+    else
+    {
+        tmin = (bbox.max.x - origin.x) * inv;
+        tmax = (bbox.min.x - origin.x) * inv;
+    }
+
+    inv = 1.0 / direction.y;
+    if (inv >= 0.0)
+    {
+        tymin = (bbox.min.y - origin.y) * inv;
+        tymax = (bbox.max.y - origin.y) * inv;
+    }
+    else
+    {
+        tymin = (bbox.max.y - origin.y) * inv;
+        tymax = (bbox.min.y - origin.y) * inv;
+    }
+
+    if (tmin > tymax || tymin > tmax)
+        return 0;
+
+    if (tymin > tmin)
+        tmin = tymin;
+    if (tymax < tmax)
+        tmax = tymax;
+
+    inv = 1.0 / direction.z;
+    if (inv >= 0.0)
+    {
+        tzmin = (bbox.min.z - origin.z) * inv;
+        tzmax = (bbox.max.z - origin.z) * inv;
+    }
+    else
+    {
+        tzmin = (bbox.max.z - origin.z) * inv;
+        tzmax = (bbox.min.z - origin.z) * inv;
+    }
+
+    if (tmin > tzmax || tzmin > tmax)
+        return 0;
+    if (tzmin > tmin)
+        tmin = tzmin;
+    if (tzmax < tmax)
+        tmax = tzmax;
+
+    *t_min = tmin;
+    *t_max = tmax;
+
+    return (tmin < t1 && tmax > t0);
+}
+
+float
+bbox_volume(boundingbox bbox)
+{
+    return (bbox.max.x - bbox.min.x) * (bbox.max.y - bbox.min.y) * (bbox.max.z - bbox.min.z);
+}

graphics/ass7/bbox.h

+#ifndef BBOX_H
+#define BBOX_H
+
+#include "types.h"
+
+//
+// An axis-aligned bounding box (AABB)
+//
+
+typedef struct
+{
+    // The two corners of the bounding box
+    vec3    min;
+    vec3    max;
+}
+boundingbox;
+
+// Create a bounding box (initialized to be empty, i.e. min > max)
+boundingbox bbox_create(void);
+
+// Update the given bounding box to include the point p
+void        bbox_update(boundingbox* bbox, vec3 p);
+
+// Return the bounding box resulting from the combination of the
+// two given bounding boxes
+boundingbox bbox_combine(boundingbox bbox1, boundingbox bbox2);
+
+// XXX ray segment, t0, t1
+// Intersect a ray with given origin and direction (and delimited by
+// t0 and t1), with the given bounding box.
+// Returns 1 (and sets t_min and t_max) when the ray intersects the bbox,
+// returns 0 otherwise
+int         bbox_intersect(float* t_min, float* t_max, boundingbox bbox,
+                vec3 origin, vec3 direction, float t0, float t1);
+
+// Return the volume of the bounding box
+float       bbox_volume(boundingbox bbox);
+
+#endif

graphics/ass7/buddha.scn

+light 9 -4.0 10.0  0.7
+light -10 2.0 15.0  0.3
+
+material 1
+
+sphere 0.65 -0.6 0.4 0.2
+sphere -0.6 0.75 0.4 0.4
+sphere -0.5 0.15 0.2 0.2
+
+ply_file ../scenes/plateau_subdivided.ply
+ply_file ../scenes/kubus.ply
+ply_file ../scenes/tetra.ply
+ply_file ../scenes/buddha.ply

graphics/ass7/bvh.h

+#ifndef BVH_H
+#define BVH_H
+
+#include <assert.h>
+#include "bbox.h"
+
+// Flag that controls if intersection testing is done with the use
+// of the BVH
+int use_bvh;
+
+//
+// One node in the Bounding Volume Hierarchy
+// We use a single data type for this to avoid having to cast
+// pointer values and such.
+//
+
+// We can't create a typedef bvh_node, as a _bvh_node contains
+// pointers to _bvh_node's (its children)
+
+#define bvh_node    struct _bvh_node
+
+struct _bvh_node
+{
+    // Boolean determining if this node is a leaf node or an inner node
+    int is_leaf;
+
+    // The bounding box for this node
+    boundingbox bbox;
+
+    // Depending on the is_leaf flag access .leaf or .inner
+    union
+    {
+        // Leaf node
+        struct
+        {
+            int         num_triangles;
+            triangle    *triangles;
+        }
+        leaf;
+
+        // Inner node
+        struct
+        {
+            // Neither of these will be NULL
+            bvh_node    *left_child, *right_child;
+        }
+        inner;
+    }
+    u;
+};
+
+// Build a BVH for the triangles in the scene.
+// This function is called from read_scene().
+void        bvh_build(void);
+
+// The root node of the hierarchy.
+// Set by bvh_build()
+bvh_node    *bvh_root;
+
+
+// Return a child node for an inner node
+// (either left or right child)
+
+static inline bvh_node*
+inner_node_left_child(const bvh_node *n)
+{
+    assert (!n->is_leaf && "Node n is not an inner node!");
+    return n->u.inner.left_child;
+}
+
+static inline bvh_node*
+inner_node_right_child(const bvh_node *n)
+{
+    assert (!n->is_leaf && "Node n is not an inner node!");
+    return n->u.inner.right_child;
+}
+
+// Return the (number of) triangles for a leaf node
+
+static inline int
+leaf_node_num_triangles(const bvh_node *n)
+{
+    assert (n->is_leaf && "Node n is not a leaf node!");
+    return n->u.leaf.num_triangles;
+}
+
+static inline triangle*
+leaf_node_triangles(const bvh_node *n)
+{
+    assert (n->is_leaf && "Node n is not a leaf node!");
+    return n->u.leaf.triangles;
+}
+
+#endif

graphics/ass7/constants.h

+#ifndef CONSTANTS_H
+#define CONSTANTS_H
+
+#include <math.h>
+#include "types.h"
+
+#define C_EPSILON 	1.0e-6
+#define C_INFINITY	1.0e6
+
+#ifndef M_PI
+#define M_PI 3.1415926535897932384626
+#endif
+
+#endif

graphics/ass7/cow.scn

+light 9 -4.0 10.0  0.7
+light -10 2.0 15.0  0.3
+
+material 1
+
+sphere 0.65 -0.6 0.4 0.2
+sphere -0.6 0.75 0.4 0.4
+sphere -0.5 0.15 0.2 0.2
+
+ply_file ../scenes/plateau_subdivided.ply
+ply_file ../scenes/kubus.ply
+ply_file ../scenes/tetra.ply
+
+ply_file ../scenes/cow.ply

graphics/ass7/intersection.c

+/* Computer Graphics, Assignment, Ray-tracing 2
+ *
+ * Student name ....
+ * Student email ...
+ * Collegekaart ....
+ * Date ............
+ * Comments ........
+ *
+ *
+ * (always fill in these fields before submitting!!)
+ */
+
+#include <math.h>
+#include <stdio.h>
+#include "intersection.h"
+#include "v3math.h"
+#include "constants.h"
+#include "scene.h"
+#include "bvh.h"
+
+// A few counters for gathering statistics on the number and types
+// of ray shot
+
+// The total number of rays
+int num_rays_shot = 0;
+
+// Number of shadow rays
+int num_shadow_rays_shot = 0;
+
+// Number of triangles tested for intersection
+int num_triangles_tested = 0;
+
+// Number of bounding boxes tested for intersection
+int num_bboxes_tested = 0;
+
+// Forward declarations
+
+static int  find_first_intersected_bvh_triangle(intersection_point* ip,
+                vec3 ray_origin, vec3 ray_direction);
+
+// Checks if the given triangle is intersected by ray with given
+// origin and direction.
+//
+// Returns 1 if there is an intersection, or 0 otherwise.
+//
+// When an intersection is found the fields of 'ip' will be filled in
+// with the relevant values.
+//
+// Note: this routine does NOT return an intersection for triangles
+// whose back side faces the ray (by definition a triangle normal
+// points to the triangle's front side).
+// I.e. we do back-face culling here ...
+//
+// Code based on Moller & Trumbore, 1997, "Fast, minimum storage
+// ray/triangle intersection"
+
+static int
+ray_intersects_triangle(intersection_point* ip, triangle tri,
+    vec3 ray_origin, vec3 ray_direction)
+{
+    vec3    edge1, edge2;
+    vec3    tvec, pvec, qvec;
+    double  det, inv_det;
+    double  t, u, v;        // u, v are barycentric coordinates
+    // t is ray parameter
+
+    num_triangles_tested++;
+
+    edge1 = v3_subtract(scene_vertices[tri.v[1]], scene_vertices[tri.v[0]]);
+    edge2 = v3_subtract(scene_vertices[tri.v[2]], scene_vertices[tri.v[0]]);
+
+    pvec = v3_crossprod(ray_direction, edge2);
+
+    det = v3_dotprod(edge1, pvec);
+
+    if (det < 1.0e-6)
+        return 0;
+
+    tvec = v3_subtract(ray_origin, scene_vertices[tri.v[0]]);
+
+    u = v3_dotprod(tvec, pvec);
+    if (u < 0.0 || u > det)
+        return 0;
+
+    qvec = v3_crossprod(tvec, edge1);
+
+    v = v3_dotprod(ray_direction, qvec);
+    if (v < 0.0 || u+v > det)
+        return 0;
+
+    t = v3_dotprod(edge2, qvec);
+
+    if (t < 0.0)
+        return 0;
+
+    inv_det = 1.0 / det;
+    t *= inv_det;
+    u *= inv_det;
+    v *= inv_det;
+
+    // We have a triangle intersection!
+    // Return the relevant intersection values.
+
+    // Compute the actual intersection point
+    ip->t = t;
+    ip->p = v3_add(ray_origin, v3_multiply(ray_direction, t));
+
+    // Compute an interpolated normal for this intersection point, i.e.
+    // we use the barycentric coordinates as weights for the vertex normals
+    ip->n = v3_normalize(v3_add(
+        v3_add(
+            v3_multiply(tri.vn[0], 1.0-u-v),
+            v3_multiply(tri.vn[1], u)
+        ),
+        v3_multiply(tri.vn[2], v)));
+
+    ip->i = v3_normalize(v3_negate(ray_direction));
+    ip->material = tri.material;
+
+    return 1;
+}
+
+// Check if the given sphere is intersected by the given ray.
+// See Shirley et.al., section 10.3.1
+// Returns 1 if there is an intersection (and sets the appropriate
+// fields of ip), or 0 otherwise.
+static int
+ray_intersects_sphere(intersection_point* ip, sphere sph,
+    vec3 ray_origin, vec3 ray_direction)
+{
+    float   A, B, C, D;
+    vec3    diff;
+    float   t_hit;
+
+    A = v3_dotprod(ray_direction, ray_direction);
+
+    diff = v3_subtract(ray_origin, sph.center);
+    B = 2.0 * v3_dotprod(diff, ray_direction);
+    C = v3_dotprod(diff, diff) - sph.radius * sph.radius;
+
+    D = B*B - 4*A*C;
+
+    if (D < 0.0)
+        return 0;
+
+    D = sqrt(D);
+
+    // We're only interested in the first hit, i.e. the one with
+    // the smallest t_hit, so we check -B-D first, followed by -B+D
+
+    t_hit = (-B - D)/(2*A);
+
+    if (t_hit < 0.0)
+    {
+        t_hit = (-B + D)/(2*A);
+        if (t_hit < 0.0)
+            return 0;
+    }
+
+    ip->t = t_hit;
+    ip->p = v3_add(ray_origin, v3_multiply(ray_direction, t_hit));
+    ip->n = v3_normalize(v3_subtract(ip->p, sph.center));
+    ip->i = v3_normalize(v3_negate(ray_direction));
+    ip->material = sph.material;
+
+    return 1;
+}
+
+// Checks for an intersection of the given ray with the triangles
+// stored in the BVH.
+//
+// Returns 1 if there is an intersection. The fields of 'ip' will be
+// set to the relevant values. The intersection returned
+// will be the one closest to the ray origin.
+//
+// Returns 0 if there are no intersections
+
+static int
+find_first_intersected_bvh_triangle(intersection_point* ip,
+    vec3 ray_origin, vec3 ray_direction)
+{
+    return 0;
+}
+
+// Returns the nearest hit of the given ray with objects in the scene
+// (either a sphere or a triangle).
+//
+// Returns 1 and sets the intersection point values if there
+// is an intersection, returns 0 otherwise.
+int
+find_first_intersection(intersection_point *ip, vec3 ray_origin, vec3 ray_direction)
+{
+    int     have_hit;
+    float   t_nearest;
+    intersection_point  ip2;
+
+    num_rays_shot++;
+
+    // We have found no hit yet
+    t_nearest = C_INFINITY;
+    have_hit = 0;
+
+    // First check against spheres in the scene
+    for (int s = 0; s < scene_num_spheres; s++)
+    {
+        // We need a second set of p and n variables, as there's the
+        // possibility that we'll overwrite a closer intersection already
+        // found
+        if (ray_intersects_sphere(&ip2, scene_spheres[s], ray_origin, ray_direction))
+        {
+            if (ip2.t < t_nearest)
+            {
+                *ip = ip2;
+                t_nearest = ip2.t;
+                have_hit = 1;
+            }
+        }
+    }
+
+    // Then check against triangles in the scene
+
+    if (use_bvh)
+    {
+        // Use the BVH to speed up intersection testing
+        if (find_first_intersected_bvh_triangle(&ip2, ray_origin, ray_direction))
+        {
+            if (ip2.t < t_nearest)
+            {
+                *ip = ip2;
+                t_nearest = ip2.t;
+                have_hit = 1;
+            }
+        }
+    }
+    else
+    {
+        // Simply iterate over all the triangles in the scene and check for intersection
+        for (int t = 0; t < scene_num_triangles; t++)
+        {
+            if (ray_intersects_triangle(&ip2, scene_triangles[t], ray_origin, ray_direction))
+            {
+                if (ip2.t < t_nearest)
+                {
+                    *ip = ip2;
+                    t_nearest = ip2.t;
+                    have_hit = 1;
+                }
+            }
+        }
+    }
+
+    return have_hit;
+}
+
+// Optimized routine for tracing a shadow ray.
+//
+// This routine doesn't return the nearest intersection, but simply
+// checks if there is any intersection.
+int
+shadow_check(vec3 ray_origin, vec3 ray_direction)
+{
+    intersection_point  ip;
+
+    num_rays_shot++;
+    num_shadow_rays_shot++;
+
+    for (int s = 0; s < scene_num_spheres; s++)
+    {
+        if (ray_intersects_sphere(&ip, scene_spheres[s], ray_origin, ray_direction))
+            return 1;
+    }
+
+    if (use_bvh)
+    {
+        // Use the BVH for speedy intersection testing
+        if (find_first_intersected_bvh_triangle(&ip, ray_origin, ray_direction))
+            return 1;
+    }
+    else
+    {
+        // Simply iterate over all the triangles in the scene and check for intersection
+        for (int t = 0; t < scene_num_triangles; t++)
+        {
+            if (ray_intersects_triangle(&ip, scene_triangles[t], ray_origin, ray_direction))
+                return 1;
+        }
+    }
+
+    return 0;
+}
+

graphics/ass7/intersection.h

+#ifndef INTERSECTION_H
+#define INTERSECTION_H
+
+#include "types.h"
+
+int     num_rays_shot;
+int     num_shadow_rays_shot;
+int     num_triangles_tested;
+int     num_bboxes_tested;
+
+int     find_first_intersection(intersection_point *ip,
+            vec3 ray_origin, vec3 ray_direction);
+
+int     shadow_check(vec3 ray_origin, vec3 ray_direction);
+
+#endif

graphics/ass7/matte.scn

+light 9 -4.0 10.0  0.7
+light -10 2.0 15.0  0.3
+
+material 1
+
+sphere 0.65 -0.6 0.4 0.2
+sphere -0.6 0.75 0.4 0.4
+sphere -0.5 0.15 0.2 0.2
+
+ply_file ../scenes/plateau.ply
+ply_file ../scenes/kubus.ply
+ply_file ../scenes/tetra.ply
+

graphics/ass7/perlin.c

+/* Implemention of Perlin's improved noise; See "Improving Noise", SIGGRAPH 2002 */
+#include <math.h>
+#include <stdio.h>
+#include "perlin.h"
+
+static int p[512];
+
+static int permutation[] = {
+    151,160,137,91,90,15,
+    131,13,201,95,96,53,194,233,7,225,140,36,103,30,69,142,8,99,37,240,
+    21,10,23,190,6,148,247,120,234,75,0,26,197,62,94,252,219,203,117,
+    35,11,32,57,177,33,88,237,149,56,87,174,20,125,136,171,168,68,175,
+    74,165,71,134,139,48,27,166,77,146,158,231,83,111,229,122,60,211,133,
+    230,220,105,92,41,55,46,245,40,244,102,143,54,65,25,63,161,1,216,
+    80,73,209,76,132,187,208,89,18,169,200,196,135,130,116,188,159,86,
+    164,100,109,198,173,186,3,64,52,217,226,250,124,123,5,202,38,147,
+    118,126,255,82,85,212,207,206,59,227,47,16,58,17,182,189,28,42,223,
+    183,170,213,119,248,152,2,44,154,163,70,221,153,101,155,167,43,
+    172,9,129,22,39,253,19,98,108,110,79,113,224,232,178,185,112,104,
+    218,246,97,228,251,34,242,193,238,210,144,12,191,179,162,241,81,51,
+    145,235,249,14,239,107,49,192,214,31,181,199,106,157,184,84,204,176,
+    115,121,50,45,127,4,150,254,138,236,205,93,222,114,67,29,24,72,243,
+    141,128,195,78,66,215,61,156,180
+};
+
+/* Function declarations */
+
+static double   fade(double t);
+static double   lerp(double t, double a, double b);
+static double   grad(int hash, double x, double y, double z);
+
+void init_noise(void)
+{
+    int i;
+    for(i = 0; i < 256 ; i++)
+        p[256+i] = p[i] = permutation[i];
+}
+
+// Signed noise
+// Result is in [-1,1]
+double snoise(double x, double y, double z)
+{
+    int   X = (int)floor(x) & 255,             /* FIND UNIT CUBE THAT */
+          Y = (int)floor(y) & 255,             /* CONTAINS POINT.     */
+          Z = (int)floor(z) & 255;
+    x -= floor(x);                             /* FIND RELATIVE X,Y,Z */
+    y -= floor(y);                             /* OF POINT IN CUBE.   */
+    z -= floor(z);
+    double  u = fade(x),                       /* COMPUTE FADE CURVES */
+            v = fade(y),                       /* FOR EACH OF X,Y,Z.  */
+            w = fade(z);
+    int  A = p[X]+Y,
+         AA = p[A]+Z,
+         AB = p[A+1]+Z, /* HASH COORDINATES OF */
+         B = p[X+1]+Y,
+         BA = p[B]+Z,
+         BB = p[B+1]+Z; /* THE 8 CUBE CORNERS, */
+
+    return lerp(w,lerp(v,lerp(u, grad(p[AA  ], x, y, z),   /* AND ADD */
+                         grad(p[BA  ], x-1, y, z)),        /* BLENDED */
+                 lerp(u, grad(p[AB  ], x, y-1, z),         /* RESULTS */
+                         grad(p[BB  ], x-1, y-1, z))),     /* FROM  8 */
+                 lerp(v, lerp(u, grad(p[AA+1], x, y, z-1 ),/* CORNERS */
+                         grad(p[BA+1], x-1, y, z-1)),      /* OF CUBE */
+                 lerp(u, grad(p[AB+1], x, y-1, z-1),
+                         grad(p[BB+1], x-1, y-1, z-1))));
+}
+
+
+// Result is in [0,1]
+double noise(double x, double y, double z)
+{
+    return 0.5 * (1.0 + snoise(x, y, z));
+}
+
+static double
+fade(double t)
+{
+    return t * t * t * (t * (t * 6 - 15) + 10);
+}
+
+static double
+lerp(double t, double a, double b)
+{
+    return a + t * (b - a);
+}
+
+static double
+grad(int hash, double x, double y, double z)
+{
+    int     h = hash & 15;       /* CONVERT LO 4 BITS OF HASH CODE */
+    double  u = h < 8 ? x : y,   /* INTO 12 GRADIENT DIRECTIONS.   */
+            v = h < 4 ? y : h==12||h==14 ? x : z;
+    return ((h&1) == 0 ? u : -u) + ((h&2) == 0 ? v : -v);
+}

graphics/ass7/perlin.h

+#ifndef PERLIN_H
+#define PERLIN_H
+
+void    init_noise(void);
+double  noise(double x, double y, double z);
+double  snoise(double x, double y, double z);
+
+#endif

graphics/ass7/plymodel.c

+#include <sys/time.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include "rply.h"
+#include "plymodel.h"
+#include "v3math.h"
+
+static int  vertex_cb(p_ply_argument argument);
+static int  normal_cb(p_ply_argument argument);
+static int  face_cb(p_ply_argument argument);
+
+static void dump(void);
+
+static int  vertex_element_idx, normal_element_idx, face_element_idx;
+
+void
+read_ply_model(const char *fname)
+{
+    int     i, j, k;
+
+    p_ply   ply = ply_open(fname, NULL);
+
+    // XXX handle better
+    if (!ply) return;
+    if (!ply_read_header(ply)) return;
+
+    ply_num_vertices = ply_set_read_cb(ply, "vertex", "x", vertex_cb, NULL, 0);
+    ply_set_read_cb(ply, "vertex", "y", vertex_cb, NULL, 0);
+    ply_set_read_cb(ply, "vertex", "z", vertex_cb, NULL, 0);
+
+    ply_set_read_cb(ply, "vertex", "nx", normal_cb, NULL, 0);
+    ply_set_read_cb(ply, "vertex", "ny", normal_cb, NULL, 0);
+    ply_set_read_cb(ply, "vertex", "nz", normal_cb, NULL, 0);
+
+    int ply_num_faces = ply_set_read_cb(ply, "face", "vertex_indices", face_cb, NULL, 0);
+
+    fprintf(stderr, "Reading %s\n", fname);
+
+    ply_vertices = malloc(ply_num_vertices * sizeof(vec3));
+    ply_normals = malloc(ply_num_vertices * sizeof(vec3));
+
+    // XXX this will overallocate if some of the faces aren't triangles
+    ply_triangles = malloc(ply_num_faces * sizeof(triangle));
+
+    // will be updated by the callbacks
+    ply_num_vertices = 0;
+    ply_num_normals = 0;
+    ply_num_triangles = 0;
+
+    vertex_element_idx = 0;
+    normal_element_idx = 0;
+    face_element_idx = 0;
+
+    // XXX
+    if (!ply_read(ply))
+    {
+        fprintf(stderr, "\nFATAL: Could not read .ply file %s!\n", fname);
+        exit(-1);
+    }
+
+    fprintf(stderr, "... %d vertices, %d normals, %d triangles\n",
+        ply_num_vertices, ply_num_normals, ply_num_triangles);
+
+    // Compute triangle normals
+
+    fprintf(stderr, "... Computing triangle normals ... ");
+
+    for (int t = 0; t < ply_num_triangles; t++)
+    {
+        i = ply_triangles[t].v[0];
+        j = ply_triangles[t].v[1];
+        k = ply_triangles[t].v[2];
+
+        // Calculate normal and store it with the triangle
+        ply_triangles[t].n = v3_normalize(
+            v3_crossprod(
+                v3_subtract(ply_vertices[j], ply_vertices[i]),
+                v3_subtract(ply_vertices[k], ply_vertices[i])));
+    }
+
+    if (ply_num_normals == ply_num_vertices)
+    {
+        // Simply copy triangle vertex normals
+        for (int t = 0; t < ply_num_triangles; t++)
+        {
+            for (i = 0; i < 3; i++)
+                ply_triangles[t].vn[i] = ply_normals[ply_triangles[t].v[i]];
+        }
+    }
+    else
+    {
+        // Incorrect number of vertex normals provided.
+        // Set triangle's vertex normals to triangle normal, i.e.
+        // flat shading
+        for (int t = 0; t < ply_num_triangles; t++)
+        {
+            for (i = 0; i < 3; i++)
+                ply_triangles[t].vn[i] = ply_triangles[t].n;
+        }
+    }
+
+    fprintf(stderr, "done\n");
+
+    // all done...
+
+    ply_close(ply);
+}
+
+int
+vertex_cb(p_ply_argument argument)
+{
+    double      v = ply_get_argument_value(argument);
+
+    switch (vertex_element_idx)
+    {
+        case 0:
+            ply_vertices[ply_num_vertices].x = v;
+            break;
+        case 1:
+            ply_vertices[ply_num_vertices].y = v;
+            break;
+        case 2:
+            ply_vertices[ply_num_vertices].z = v;
+            break;
+    }
+
+    if (vertex_element_idx == 2)
+    {
+        ply_num_vertices++;
+        vertex_element_idx = 0;
+    }
+    else
+        vertex_element_idx++;
+
+    return 1;
+}
+
+int
+normal_cb(p_ply_argument argument)
+{
+    double      v = ply_get_argument_value(argument);
+
+    switch (normal_element_idx)
+    {
+        case 0:
+            ply_normals[ply_num_normals].x = v;
+            break;
+        case 1:
+            ply_normals[ply_num_normals].y = v;
+            break;
+        case 2:
+            ply_normals[ply_num_normals].z = v;
+            break;
+    }
+
+    if (normal_element_idx == 2)
+    {
+        ply_num_normals++;
+        normal_element_idx = 0;
+    }
+    else
+        normal_element_idx++;
+
+    return 1;
+}
+
+int
+face_cb(p_ply_argument argument)
+{
+    int32        length, value_index;
+
+    ply_get_argument_property(argument, NULL, &length, &value_index);
+
+    // Skip any polygon that isn't a triangle
+    if (length == 3 && value_index >= 0)
+    {
+        int idx = (int)ply_get_argument_value(argument);
+        ply_triangles[ply_num_triangles].v[face_element_idx] = idx;
+
+        if (face_element_idx == 2)
+        {
+            ply_num_triangles++;
+            face_element_idx = 0;
+        }
+        else
+            face_element_idx++;
+    }
+
+    return 1;
+}
+
+static void
+dump(void)
+{
+    int i;
+
+    for (i = 0; i < ply_num_vertices; i++)
+    {
+        if (ply_num_normals == ply_num_vertices)
+        {
+            printf("[%d] %.6f, %.6f, %.6f  (%.6f, %.6f, %.6f)\n", i,
+                ply_vertices[i].x, ply_vertices[i].y, ply_vertices[i].z,
+                ply_normals[i].x, ply_normals[i].y, ply_normals[i].z);
+        }
+        else
+        {
+            printf("[%d] %.6f, %.6f, %.6f\n", i, ply_vertices[i].x,
+                ply_vertices[i].y, ply_vertices[i].z);
+        }
+    }
+
+    for (i = 0; i < ply_num_triangles; i++)
+    {
+        printf("[%d] %d %d %d\n", i, ply_triangles[i].v[0],
+            ply_triangles[i].v[1], ply_triangles[i].v[2]);
+    }
+}
+

graphics/ass7/plymodel.h

+#ifndef PLYMODEL_H
+#define PLYMODEL_H
+
+#include "types.h"
+
+int         ply_num_triangles, ply_num_vertices, ply_num_normals;
+
+triangle    *ply_triangles;
+vec3        *ply_vertices;
+vec3        *ply_normals;
+
+void        read_ply_model(const char *fname);
+
+#endif

graphics/ass7/quat.c

+#include <math.h>
+#include "quat.h"
+#include "constants.h"
+
+// Create a quaternion that represents a rotation of
+// angle degrees around the given axis.
+quat
+quat_create_rotation(vec3 axis, float angle)
+{
+    float   sin_a;
+    quat    res;
+
+    axis = v3_normalize(axis);
+
+    sin_a = sin(0.5 * (angle/180.0*M_PI));
+
+    res.x = axis.x * sin_a;
+    res.y = axis.y * sin_a;
+    res.z = axis.z * sin_a;
+    res.w = cos(0.5 * (angle/180.0*M_PI));
+
+    return quat_normalize(res);
+}
+
+// Create a quaternion representation of point p
+quat
+quat_create_point(vec3 p)
+{
+    quat    res;
+
+    res.x = p.x;
+    res.y = p.y;
+    res.z = p.z;
+    res.w = 0.0;
+
+    return res;
+}
+
+
+// Calculate quaternion magnitude
+float
+quat_magnitude(quat q)
+{
+    return sqrt(q.x*q.x + q.y*q.y + q.z*q.z + q.w*q.w);
+}
+
+// Normalize a quaternion
+quat
+quat_normalize(quat q)
+{
+    quat    res = q;
+    float   inv_mag = 1.0 / quat_magnitude(q);
+
+    res.x *= inv_mag;
+    res.y *= inv_mag;
+    res.z *= inv_mag;
+    res.w *= inv_mag;
+
+    return res;
+}
+
+
+quat
+quat_conjugate(quat q)
+{
+    quat    res = q;
+
+    res.x = -q.x;
+    res.y = -q.y;
+    res.z = -q.z;
+
+    return res;
+}
+
+quat
+quat_inverse(quat q)
+{
+    float   invlen = 1.0 / quat_magnitude(q);
+    quat    res = quat_conjugate(q);
+
+    invlen = invlen * invlen;
+
+    res.x *= invlen;
+    res.y *= invlen;
+    res.z *= invlen;
+    res.w *= invlen;
+
+    return res;
+}
+
+
+quat
+quat_multiply(quat q, quat r)
+{
+    quat    res;
+
+	res.w = q.w*r.w - q.x*r.x - q.y*r.y - q.z*r.z;
+	res.x = q.w*r.x + q.x*r.w + q.y*r.z - q.z*r.y;
+	res.y = q.w*r.y + q.y*r.w + q.z*r.x - q.x*r.z;
+	res.z = q.w*r.z + q.z*r.w + q.x*r.y - q.y*r.x;
+
+    return res;
+}
+
+

graphics/ass7/quat.h

+#ifndef QUAT_H
+#define QUAT_H
+
+#include "v3math.h"
+
+// A quaternion, used to represent rotations around arbitrary angles
+typedef struct
+{
+    float x, y, z, w;
+}
+quat;
+
+// Create a quaternion that represents a rotation of
+// angle degrees around the given axis.
+quat    quat_create_rotation(vec3 axis, float angle);
+
+// Create a quaternion representation of point p
+quat    quat_create_point(vec3 p);
+
+// Calculate quaternion magnitude
+float   quat_magnitude(quat q);
+
+// Normalize a quaternion
+quat    quat_normalize(quat q);
+
+// Calculate the conjugate (i.e. (-v, w) for quaternion (v, w))
+quat    quat_conjugate(quat q);
+
+// Calculate the inverse of a quaterion
+quat    quat_inverse(quat q);
+
+// Transform
+vec3    quat_transform_vector(quat q, vec3 v);
+
+quat    quat_multiply(quat q, quat r);
+
+#endif

graphics/ass7/scene.c

+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include "scene.h"
+#include "plymodel.h"
+#include "bvh.h"
+
+int         scene_num_triangles, scene_num_vertices;
+triangle    *scene_triangles=NULL;
+vec3        *scene_vertices=NULL;
+
+int     	scene_num_spheres;
+sphere  	*scene_spheres=NULL;
+
+int         scene_num_lights;
+light       *scene_lights=NULL;
+float       scene_ambient_light = 0.05;
+
+vec3    	scene_background_color = { 1.0, 1.0, 1.0 };
+
+vec3        scene_camera_position, scene_camera_lookat;
+
+void
+read_scene(const char *fname)
+{
+	FILE 	*f;
+	char	line[256];
+	int		i;
+	float	x, y, z, r;
+	char	s[256];
+	int		material;
+
+	f = fopen(fname, "rt");
+
+	if (!f)
+	{
+        fprintf(stderr, "\nFATAL: Could not open scene file %s!\n", fname);
+        exit(-1);
+    }
+
+    printf("Reading scene %s\n", fname);
+
+	scene_num_triangles = scene_num_vertices = 0;
+	scene_num_spheres = 0;
+
+	material = 0;
+
+	// Bare bones, but usable, parsing
+	fgets(line, 255, f);
+
+	while (!feof(f))
+	{
+        if (line[0] == '/' && line[1] == '/')
+        {
+            // Comment, skip
+        }
+        else if (sscanf(line, "material %d\n", &material) == 1)
+		{
+			// Nothing to do, value already stored...
+		}
+        else if (sscanf(line, "light %f %f %f %f\n", &x, &y, &z, &r) == 4)
+        {
+			scene_lights = realloc(scene_lights, (scene_num_lights+1)*sizeof(light));
+			scene_lights[scene_num_lights].position.x = x;
+			scene_lights[scene_num_lights].position.y = y;
+			scene_lights[scene_num_lights].position.z = z;
+			scene_lights[scene_num_lights].intensity = r;
+			scene_num_lights++;
+        }
+		else if (sscanf(line, "sphere %f %f %f %f\n", &x, &y, &z, &r) == 4)
+		{
+			scene_spheres = realloc(scene_spheres, (scene_num_spheres+1)*sizeof(sphere));
+			scene_spheres[scene_num_spheres].center.x = x;
+			scene_spheres[scene_num_spheres].center.y = y;
+			scene_spheres[scene_num_spheres].center.z = z;
+			scene_spheres[scene_num_spheres].radius = r;
+			scene_spheres[scene_num_spheres].material = material;
+			scene_num_spheres++;
+		}
+		else if (sscanf(line, "ply_file %s\n", s) == 1)
+		{
+            // Read triangle model from a .ply file, add the triangles
+            // in the file to the global list of triangles in the scene
+
+            printf("[%s]\n", s);
+
+			read_ply_model(s);
+
+			// Allocate memory to hold the new vertices and triangles
+
+			scene_vertices = realloc(scene_vertices,
+                (scene_num_vertices + ply_num_vertices)*sizeof(vec3));
+			scene_triangles = realloc(scene_triangles,
+                (scene_num_triangles + ply_num_triangles)*sizeof(triangle));
+
+			// Copy the vertices
+			memcpy(scene_vertices + scene_num_vertices,
+                ply_vertices, ply_num_vertices*sizeof(vec3));
+
+			// Copy the triangles, AND RENUMBER THEIR VERTEX INDICES
+            // (as we're going to add the model triangles and vertices
+            // to a scene-wide list of triangles and vertices)
+			for (i = 0; i < ply_num_triangles; i++)
+			{
+				scene_triangles[scene_num_triangles+i] = ply_triangles[i];
+
+				// Renumber the vertex indices
+				scene_triangles[scene_num_triangles+i].v[0] += scene_num_vertices;
+				scene_triangles[scene_num_triangles+i].v[1] += scene_num_vertices;
+				scene_triangles[scene_num_triangles+i].v[2] += scene_num_vertices;
+
+				// Assign the material
+				scene_triangles[scene_num_triangles+i].material = material;
+			}
+
+			// Update the total numer of vertices in the scene
+			scene_num_vertices += ply_num_vertices;
+			scene_num_triangles += ply_num_triangles;
+		}
+        else if (line[0] != '\n')
+        {
+            fprintf(stderr, "Ignoring unknown line '%s'\n", line);
+        }
+
+		fgets(line, 255, f);
+	}
+
+	fclose(f);
+
+    // Build the BVH for the triangles read
+    bvh_build();
+}

graphics/ass7/scene.h

+#ifndef SCENE_H
+#define SCENE_H
+
+#include "types.h"
+
+extern int      scene_num_triangles, scene_num_vertices;
+extern triangle *scene_triangles;
+extern vec3     *scene_vertices;
+
+extern int     	scene_num_spheres;
+extern sphere   *scene_spheres;
+
+extern int      scene_num_lights;
+extern light    *scene_lights;
+
+extern float    scene_ambient_light;
+extern vec3    	scene_background_color;
+
+extern vec3     scene_camera_position;
+extern vec3     scene_camera_lookat;
+
+void		    read_scene(const char *fname);
+
+#endif

graphics/ass7/shaders.c

+#include <math.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include "shaders.h"
+#include "perlin.h"
+#include "v3math.h"
+#include "intersection.h"
+#include "scene.h"
+#include "quat.h"
+#include "constants.h"
+
+// Always return the same color. This shader does no real computations
+// based on normal, light position, etc. As such, it merely creates
+// a "silhouette" of an object.
+
+vec3
+shade_constant(intersection_point ip)
+{
+    return v3_create(1, 0, 0);
+}
+
+vec3
+shade_matte(intersection_point ip)
+{
+    vec3   	l_dir;
+    float   f;
+    float   res;
+
+    res = scene_ambient_light;
+
+    for (int l = 0; l < scene_num_lights; l++)
+    {
+        l_dir = v3_subtract(scene_lights[l].position, ip.p);
+
+        // Shoot a shadow ray. If the point to shade is in shadow w.r.t.
+        // this light then there is no contribution from the light
+        if (shadow_check(v3_add(ip.p, v3_multiply(l_dir, 0.0001)), l_dir))
+            continue;
+
+        // Compute and add the contribution from this light
+        l_dir = v3_normalize(l_dir);
+
+        f = v3_dotprod(ip.n, l_dir);
+        if (f > 0.0)
+            res += f * scene_lights[l].intensity;
+    }
+
+    if (res > 1.0)
+        res = 1.0;
+
+    return v3_create(res, res, res);
+}
+
+// Returns the shaded color for the given point to shade.
+// Calls the relevant shading function based on the material index.
+vec3
+shade(intersection_point ip)
+{
+	switch (ip.material)
+	{
+        case 0:
+            return shade_constant(ip);
+		case 1:
+			return shade_matte(ip);
+        default:
+            return shade_constant(ip);
+	}
+}
+
+// Determine the surface color for the first object intersected by
+// the given ray, or return the scene background color when no
+// intersection is found
+vec3
+ray_color(int level, vec3 ray_origin, vec3 ray_direction)
+{
+    intersection_point  ip;
+
+    // If this ray has been reflected too many times, simply
+    // return the background color.
+    if (level >= 3)
+        return scene_background_color;
+
+    // Check if the ray intersects anything in the scene
+    if (find_first_intersection(&ip, ray_origin, ray_direction))
+    {
+        // Shade the found intersection point
+        ip.ray_level = level;
+        return shade(ip);
+    }
+
+    // Nothing was hit, return background color
+    return scene_background_color;
+}

graphics/ass7/shaders.h

+#ifndef SHADERS_H
+#define SHADERS_H
+
+#include "types.h"
+
+// Determine the surface color for the first object intersected by
+// the given ray, or return the scene background color when no
+// intersection is found
+vec3    ray_color(int level, vec3 ray_origin, vec3 ray_direction);
+
+#endif

graphics/ass7/silhouette.scn

+light 9 -4.0 10.0  0.7
+light -10 2.0 15.0  0.3
+
+material 0
+
+sphere 0.65 -0.6 0.4 0.2
+sphere -0.6 0.75 0.4 0.4
+sphere -0.5 0.15 0.2 0.2
+
+ply_file ../scenes/plateau.ply
+ply_file ../scenes/kubus.ply
+ply_file ../scenes/tetra.ply
+

graphics/ass7/simple.scn

+light 9 -4.0 10.0  0.7
+light -10 2.0 15.0  0.3
+
+material 1
+
+ply_file ../scenes/plateau.ply
+ply_file ../scenes/kubus.ply
+ply_file ../scenes/tetra.ply
+ply_file ../scenes/cone.ply

graphics/ass7/types.h

+#ifndef TYPES_H
+#define TYPES_H
+
+#include "v3math.h"
+
+typedef unsigned char   byte;
+
+typedef struct
+{
+    // indices of vertices
+    int     v[3];
+
+	// the triangle's normal
+    vec3    n;
+
+    // the triangle's vertex normals
+    vec3    vn[3];
+
+	// the material index for this triangle
+    int     material;
+}
+triangle;
+
+typedef struct
+{
+    vec3    center;
+    float   radius;
+
+	int     material;
+}
+sphere;
+
+typedef struct
+{
+    // The ray t-value at which intersection occurred, used to find
+    // the closest hit
+    float   t;
+
+    // intersection point
+    vec3    p;
+    // surface normal at p
+    vec3    n;
+    // direction from which the ray responsible for this intersection
+    // came. points away from the surface!
+    vec3    i;
+
+    // the intersected surface's material index
+    int     material;
+
+    // level of the ray that caused this intersection;
+    // i.e. level 0 is a camera ray, level 1 is is a ray that has
+    // been reflected once, etc
+    int     ray_level;
+}
+intersection_point;
+
+typedef struct
+{
+    vec3    position;
+    float   intensity;      // 0 - 1
+}
+light;
+
+typedef struct
+{
+    float   azimuth;
+    float   rot_z;
+    float   distance;
+}
+viewpoint;
+
+#endif

graphics/ass7/v3math.h

+#ifndef V3MATH_H
+#define V3MATH_H
+
+#include <math.h>
+
+typedef struct
+{
+    float   x, y, z;
+}
+vec3;
+
+// Create a new 3-vector of floats
+static inline vec3
+v3_create(float x, float y, float z)
+{
+    vec3 res;
+    res.x = x;
+    res.y = y;
+    res.z = z;
+    return res;
+}
+
+// Return -a
+static inline vec3
+v3_negate(vec3 a)
+{
+    vec3 res;
+    res.x = - a.x;
+    res.y = - a.y;
+    res.z = - a.z;
+    return res;
+}
+
+// Return a + b
+static inline vec3
+v3_add(vec3 a,vec3 b)
+{
+	vec3 res;
+	res.x = a.x+b.x;
+	res.y = a.y+b.y;
+	res.z = a.z+b.z;
+	return res;
+}
+
+// Return a - b
+static inline vec3
+v3_subtract(vec3 a, vec3 b)
+{
+	vec3 res;
+	res.x = a.x-b.x;
+	res.y = a.y-b.y;
+	res.z = a.z-b.z;
+	return res;
+}
+
+// Return a / |a|
+static inline vec3
+v3_normalize(vec3 a)
+{
+	vec3 res;
+	double l = sqrt(a.x*a.x + a.y*a.y + a.z*a.z);
+
+	res = a;
+	res.x /= l;
+	res.y /= l;
+	res.z /= l;
+
+	return res;
+}
+
+// Return a ^ b
+static inline vec3
+v3_crossprod(vec3 a, vec3 b)
+{
+	vec3 res;
+	res.x = a.y*b.z - a.z*b.y;
+	res.y = a.z*b.x - a.x*b.z;
+	res.z = a.x*b.y - a.y*b.x;
+	return res;
+}
+
+// Return a * b
+static inline float
+v3_dotprod(vec3 a, vec3 b)
+{
+	return a.x*b.x + a.y*b.y + a.z*b.z;
+}
+
+// Return a*s
+static inline vec3
+v3_multiply(vec3 a, float s)
+{
+	vec3 res;
+	res.x = a.x*s;
+	res.y = a.y*s;
+	res.z = a.z*s;
+	return res;
+}
+
+// Return |a|
+static inline float
+v3_length(vec3 a)
+{
+	return sqrt(a.x*a.x + a.y*a.y + a.z*a.z);
+}
+
+// Return the i-th component of a, i.e. for i=0
+// this returns a.x
+static inline float
+v3_component(vec3 a, int i)
+{
+    if (i == 0)
+        return a.x;
+    else if (i == 1)
+        return a.y;
+    else
+        return a.z;
+}
+
+// Set the i-th component of a to v
+static inline vec3
+v3_set_component(vec3 a, int i, float v)
+{
+    vec3    res = a;
+
+    if (i == 0)
+        res.x = v;
+    else if (i == 1)
+        res.y = v;
+    else
+        res.z = v;
+
+    return res;
+}
+
+#endif