Added Graphics ass. 10 and Linear Algebra ass. 3

graphics/ass10/Makefile

+CC=gcc
+
+WARNING_FLAGS=-Wall -Wextra -Werror-implicit-function-declaration -Wshadow -Wstrict-prototypes -pedantic-errors
+CFLAGS=-g -O2 -std=c99 $(WARNING_FLAGS)
+LDFLAGS=-g -lGL -lglut -lGLU
+
+.c.o:
+	$(CC) -c $(CFLAGS) $<
+
+all: main
+main: geometry.o main.o polys.o ppmio.o
+	$(CC) $(LDFLAGS) -o main geometry.o main.o polys.o ppmio.o
+
+clean:
+	rm -f *.o
+	rm -f main
+
+geometry.o : geometry.h polys.h geometry.c
+ppmio.o    : ppmio.c
+geometry.o : polys.h geometry.h
+ppmio.o    : ppmio.h
+main.o     : geometry.h polys.h ppmio.h main.c
+polys.o    : polys.h polys.c

graphics/ass10/geometry.c

+/* Computer Graphics, Assignment, Texture Mapping
+ *
+ * Filename ........ geometry.c
+ * Description .....
+ * Date ............ 22.11.2007
+ * Created by ...... Paul Melis
+ *
+ * Student name ....
+ * Student email ...
+ * Collegekaart ....
+ * Date ............
+ * Comments ........
+ *
+ * (always fill in these fields before submitting!!)
+ */
+
+#include <math.h>
+#include <stdio.h>
+#include <stdlib.h>
+
+#include "geometry.h"
+#include "polys.h"
+
+#ifndef M_PI
+#define M_PI 3.141592653589793
+#endif
+
+static void
+setSpherePoint(vec3 *p, vec3* n, vec3* t, int latitude, int longitude,
+    double sx, double sy, double sz, double ox, double oy, double oz)
+{
+    double  dToR = M_PI / 180.0;
+    double  len;
+
+    // Set vertex position
+    p->x = ox + sin(longitude * dToR) * cos(latitude * dToR) * sx;
+    p->y = oy + sin(latitude * dToR) * sy;
+    p->z = oz + cos(longitude * dToR) * cos(latitude * dToR) * sz;
+
+    // Set texture coordinate
+    t->x = 0.0;
+    t->y = 0.0;
+
+    // calculate normal, this actually doesn't take the sphere size
+    // per axis into account, but should still be usable
+    n->x = p->x - ox;
+    n->y = p->y - oy;
+    n->z = p->z - oz;
+
+    len = n->x*n->x + n->y*n->y + n->z*n->z;
+    n->x /= len;
+    n->y /= len;
+    n->z /= len;
+}
+
+void
+createSphere(polys * list, double sx, double sy, double sz, double ox,
+           double oy, double oz, double r, double g, double b)
+{
+    int     latitude, longitude;
+    poly    p;
+
+    // prepare poly variable, as the color values are the same for all generated polys
+    p.points = 4;
+    p.color[0] = r;
+    p.color[1] = g;
+    p.color[2] = b;
+    p.color[3] = 0;
+
+    // sweep over sphere's surface and generate polys
+    for (latitude = -90; latitude < 90; latitude += 10)
+    {
+        for (longitude = 0; longitude < 360; longitude += 10)
+        {
+            setSpherePoint(&(p.pts[0]), &(p.normal[0]), &(p.tcoord[0]),
+                latitude, longitude, sx, sy, sz, ox, oy, oz);
+            setSpherePoint(&(p.pts[1]), &(p.normal[1]), &(p.tcoord[1]),
+                latitude+10, longitude, sx, sy, sz, ox, oy, oz);
+            setSpherePoint(&(p.pts[2]), &(p.normal[2]), &(p.tcoord[2]),
+                latitude+10, longitude+10, sx, sy, sz, ox, oy, oz);
+            setSpherePoint(&(p.pts[3]), &(p.normal[3]), &(p.tcoord[3]),
+                latitude, longitude+10, sx, sy, sz, ox, oy, oz);
+
+            AddPolyToPolylist(list, p);
+        }
+    }
+}
+
+// Calculate the coordinates, normal vector and texture coordinates for
+// a hemisphere point at the given latitude and longitude (in degrees).
+// The radius of the hemisphere is s, the center point (at the base of the
+// hemisphere) is ox,oy,oz.
+static void
+setHemispherePoint(vec3 *p, vec3* n, vec3* t, int latitude, int longitude,
+    double s, double ox, double oy, double oz)
+{
+    double  dToR = M_PI / 180.0;
+    double  len;
+
+    // Set vertex position
+    p->x = ox + sin(longitude * dToR) * cos(latitude * dToR) * s;
+    p->y = oy + sin(latitude * dToR) * s;
+    p->z = oz + cos(longitude * dToR) * cos(latitude * dToR) * s;
+
+    // Set texture coordinate
+    t->x = 0.0;
+    t->y = 0.0;
+
+    // calculate normal
+    n->x = p->x - ox;
+    n->y = p->y - oy;
+    n->z = p->z - oz;
+
+    len = n->x*n->x + n->y*n->y + n->z*n->z;
+    n->x /= len;
+    n->y /= len;
+    n->z /= len;
+}
+
+void
+createHemisphere(polys * list, double s, double ox, double oy, double oz,
+    double r, double g, double b)
+{
+    int     latitude, longitude;
+    poly    p;
+
+    // prepare poly, as these values are the same for all generated polys
+    p.points = 4;
+    p.color[0] = r;
+    p.color[1] = g;
+    p.color[2] = b;
+    p.color[3] = 0;
+
+    // sweep over sphere's surface and generate polys
+    for (latitude = 0; latitude < 90; latitude += 10)
+    {
+        for (longitude = 0; longitude < 360; longitude += 10)
+        {
+            setHemispherePoint(&(p.pts[0]), &(p.normal[0]), &(p.tcoord[0]),
+                latitude, longitude, s, ox, oy, oz);
+            setHemispherePoint(&(p.pts[1]), &(p.normal[1]), &(p.tcoord[1]),
+                latitude, longitude+10, s, ox, oy, oz);
+            setHemispherePoint(&(p.pts[2]), &(p.normal[2]), &(p.tcoord[2]),
+                latitude+10, longitude+10, s, ox, oy, oz);
+            setHemispherePoint(&(p.pts[3]), &(p.normal[3]), &(p.tcoord[3]),
+                latitude+10, longitude, s, ox, oy, oz);
+
+            AddPolyToPolylist(list, p);
+        }
+    }
+}
+
+// Create a cylinder along the Y axis whose base center point is
+// at (ox, oy, oz), having the given radius and height.
+// Use the given color for the generated polygons.
+void
+createCylinder(polys * list, double radius, double height,
+    double ox, double oy, double oz,
+    double r, double g, double b)
+{
+    int     longitude;
+    GLfloat x, z;
+    GLfloat x2, z2;
+    double  dToR = M_PI / 180.0;
+    double  len;
+    poly    p;
+    int     i;
+
+    // prepare poly datastructure, as these values are the same for all generated polys
+    p.points = 4;
+    p.color[0] = r;
+    p.color[1] = g;
+    p.color[2] = b;
+    p.color[3] = 0;
+
+    // sweep around cylinder axis
+    for (longitude = 0; longitude < 360; longitude += 10)
+    {
+        x = ox + sin(longitude * dToR) * radius;
+        z = oz + cos(longitude * dToR) * radius;
+
+        x2 = ox + sin((longitude+10) * dToR) * radius;
+        z2 = oz + cos((longitude+10) * dToR) * radius;
+
+        p.pts[0].x = x;  p.pts[0].y = oy;           p.pts[0].z = z;
+        p.pts[1].x = x;  p.pts[1].y = oy+height;    p.pts[1].z = z;
+        p.pts[2].x = x2; p.pts[2].y = oy+height;    p.pts[2].z = z2;
+        p.pts[3].x = x2; p.pts[3].y = oy;           p.pts[3].z = z2;
+
+        for (i = 0; i < 4; i++)
+        {
+            p.normal[i].x = p.pts[i].x - ox;
+            p.normal[i].y = 0.0;
+            p.normal[i].z = p.pts[i].z - oz;
+
+            len = p.normal[i].x*p.normal[i].x + p.normal[i].z*p.normal[i].z;
+            p.normal[i].x /= len;
+            p.normal[i].z /= len;
+
+            // Set texture coordinate
+            p.tcoord[i].x = p.tcoord[i].y = 0.0;
+        }
+
+        AddPolyToPolylist(list, p);
+    }
+}
+
+// Read the vertex data for a single polygon and store it in 'p'.
+// Called from loadPolygonalObject() after a 'p ...' line (+ color)
+// has been read.
+static void
+readPolyVertices(poly *p, FILE *f, const vec3 *vertex_list, int num_vertices_to_read,
+    int max_vertex_index, GLuint texture_name)
+{
+    char    line[1024];
+    int     num_read;
+    int     pidx;
+    int     vidx;
+    float   s, t;
+
+    p->points = num_vertices_to_read;
+
+    num_read = 0;
+    pidx = 0;
+    while (num_read < num_vertices_to_read)
+    {
+        if (fgets(line, 1024, f) == NULL)
+        {
+            printf("Unexpected end-of-file reading poly vertices!\n");
+            exit(-1);
+        }
+
+        if (sscanf(line, "%d %f %f\n", &vidx, &s, &t) == 3)
+        {
+            if (vidx > max_vertex_index)
+            {
+                printf("While reading %d polygon vertices, line:\n", num_vertices_to_read);
+                printf("%s\n", line);
+                printf("Vertex index %d is out of range (should be <= %d)\n", vidx, max_vertex_index);
+                exit(-1);
+            }
+
+            //printf("%d %f %f %f %f %f\n",  vidx,  nx,  ny,  nz,  s,  t);
+            p->pts[pidx] = vertex_list[vidx];
+
+            p->texture_id = texture_name;
+
+            p->tcoord[pidx].x = s;
+            p->tcoord[pidx].y = t;
+
+            num_read++;
+            pidx++;
+        }
+        else if (line[0] && line[0] != '#')
+        {
+            printf("Expected line declaring polygon vertex, got instead:\n");
+            printf("%s\n", line);
+            exit(-1);
+        }
+
+    }
+
+}
+
+// Read a polygonal object from a .OBJ file.
+// Scale the input coordinates uniformly with a factor s followed by
+// a translation (tx,ty,tz).
+void
+loadPolygonalObject(polys* list, const char *objfile, GLuint *texture_names,
+    double s, double tx, double ty, double tz)
+{
+    FILE    *f;
+    char    line[1024];
+    float   x, y, z;
+    vec3    *vertices;
+    int     num_vertices;
+    int     num_vertices_allocated;
+    poly    p;
+    float   colr, colg, colb;
+    int     num_poly_vertices, texture;
+    int     i;
+    float   len;
+    vec3    v, w, n;
+
+    f = fopen(objfile, "rt");
+
+    if (!f)
+    {
+        fprintf(stderr, "loadPolygonalObject(): Could not open file '%s' for reading!\n", objfile);
+        exit(-1);
+    }
+
+    printf("Reading %s\n", objfile);
+
+    num_vertices = 0;
+    num_vertices_allocated = 1024;
+    vertices = malloc(num_vertices_allocated * sizeof(vec3));
+
+    while (fgets(line, 1024, f))
+    {
+        if (sscanf(line, "v %f %f %f\n", &x, &y, &z))
+        {
+            // vertex
+            //printf("vertex: %f %f %f\n", x, y, z);
+
+            // current vertex list is full, add more space
+            if (num_vertices == num_vertices_allocated)
+            {
+                num_vertices_allocated = (int)(1.5 * num_vertices_allocated);
+                vertices = realloc(vertices, num_vertices_allocated * sizeof(vec3));
+                if (vertices == NULL)
+                {
+                    printf("No memory to hold vertices!\n");
+                    exit(-1);
+                }
+            }
+
+            // store vertex, after scaling and translation
+            vertices[num_vertices].x = tx + s*x;
+            vertices[num_vertices].y = ty + s*y;
+            vertices[num_vertices].z = tz + s*z;
+
+            num_vertices++;
+        }
+        else if (sscanf(line, "p %d %d\n", &num_poly_vertices, &texture))
+        {
+            if (num_poly_vertices < 3)
+            {
+                printf("In %s, line:\n", objfile);
+                printf("%s\n", line);
+                printf("Need at least 3 vertices to define a polygon!\n");
+                exit(-1);
+            }
+            else if (num_poly_vertices > 4)
+            {
+                printf("In %s, line:\n", objfile);
+                printf("%s\n", line);
+                printf("Polygons may only consist of 3 or 4 vertices, sorry\n");
+                printf("(this is due to OpenGL not handling concave polygons very well)\n");
+                exit(-1);
+            }
+
+            // read polygon color
+            if (fscanf(f, "%f %f %f\n", &colr, &colg, &colb) == EOF)
+            {
+                printf("Could not read polygon color!\n");
+                exit(-1);
+            }
+
+            p.color[0] = colr;
+            p.color[1] = colg;
+            p.color[2] = colb;
+            p.color[3] = 0.0;
+
+            // read vertices making up poly
+            readPolyVertices(&p, f, vertices, num_poly_vertices, num_vertices-1, texture_names[texture]);
+
+            // Use first 3 polygon vertices to calculate the polygon normal vector
+            // (i.e. assumes polygon is planar)
+
+            v.x = p.pts[2].x - p.pts[1].x;
+            v.y = p.pts[2].y - p.pts[1].y;
+            v.z = p.pts[2].z - p.pts[1].z;
+
+            w.x = p.pts[0].x - p.pts[1].x;
+            w.y = p.pts[0].y - p.pts[1].y;
+            w.z = p.pts[0].z - p.pts[1].z;
+
+            n.x = v.y * w.z - v.z * w.y;
+            n.y = v.z * w.x - v.x * w.z;
+            n.z = v.x * w.y - v.y * w.x;
+
+            len = sqrt(n.x * n.x + n.y * n.y + n.z * n.z);
+            n.x /= len;
+            n.y /= len;
+            n.z /= len;
+
+            // Set vertex normals to polygon normal
+
+            for (i = 0; i < num_poly_vertices; i++)
+                p.normal[i] = n;
+
+            // Add polygon
+
+            AddPolyToPolylist(list, p);
+        }
+    }
+
+    //printf("%d vertices read\n", num_vertices);
+
+    fclose(f);
+
+    free(vertices);
+}
+

graphics/ass10/geometry.h

+#ifndef GEOMETRY_H
+#define GEOMETRY_H
+
+#include "polys.h"
+
+// Create a sphere centered at (ox, oy, oz), having sizes
+// sx, sy and sz in X, Y and Z respectively. Use the given color.
+void
+createSphere(polys * list, double sx, double sy, double sz, double ox,
+             double oy, double oz, double r, double g, double b);
+
+// Create a hemisphere whose base point is at (ox, oy, oz), having radius s.
+// Use the given color.
+void
+createHemisphere(polys * list, double s, double ox, double oy, double oz,
+    double r, double g, double b);
+
+// Create a cylinder along the Y axis whose base center point is
+// at (ox, oy, oz), having the given radius and height.
+// Use the given color for the generated polygons.
+void
+createCylinder(polys * list, double radius, double height,
+    double ox, double oy, double oz,
+    double r, double g, double b);
+
+// Read a polygonal object from a .OBJ file.
+// Scale the input coordinates uniformly with a factor s followed by
+// a translation (tx,ty,tz).
+void
+loadPolygonalObject(polys * list, const char *objfile, GLuint *texture_names,
+                    double s, double tx, double ty, double tz);
+
+
+#endif

graphics/ass10/ground.obj

+v -100 0 -50
+v 100 0 -50
+v 100 0 50
+v -100 0 50
+
+# (0=grass texture)
+p 4 0
+0.5 1.0 0.5
+0  0.0 0.0
+1  0.0 0.0
+2  0.0 0.0
+3  0.0 0.0

graphics/ass10/house.obj

+# vertices on the ground, forming the rectangular base of the house
+v 2.5 0.0 1.0
+v 2.5 0.0 -1.0
+v -2.5 0.0 -1.0
+v -2.5 0.0 1.0
+
+# vertices directly above the ground vertices
+v 2.5 1.6 1.0
+v 2.5 1.6 -1.0
+v -2.5 1.6 -1.0
+v -2.5 1.6 1.0
+
+# top vertices of the two "pointed" walls
+v 2.5 2.6 0.0
+v -2.5 2.6 0.0
+
+# roof vertices
+v 2.8 2.6 0.0
+v -2.8 2.6 0.0
+v -2.8 1.4 1.2
+v 2.8 1.4 1.2
+v 2.8 1.4 -1.2
+v -2.8 1.4 -1.2
+
+
+# side walls (2=wall texture)
+p 4 2
+0.7 0.7 0.7
+0  0.0 0.0
+1  0.0 0.0
+5  0.0 0.0
+4  0.0 0.0
+p 3 2
+0.7 0.7 0.7
+5  0.0 0.0
+8  0.0 0.0
+4  0.0 0.0
+
+p 4 2
+0.7 0.7 0.7
+2  0.0 0.0
+3  0.0 0.0
+7  0.0 0.0
+6  0.0 0.0
+p 3 2
+0.7 0.7 0.7
+7  0.0 0.0
+9  0.0 0.0
+6  0.0 0.0
+
+p 4 2
+0.7 0.7 0.7
+3  0.0 0.0
+0  0.0 0.0
+4  0.0 0.0
+7  0.0 0.0
+
+p 4 2
+0.7 0.7 0.7
+2  0.0 0.0
+6  0.0 0.0
+5  0.0 0.0
+1  0.0 0.0
+
+# roof polygons (1=roof texture)
+
+p 4 1
+0.7 0.7 0.7
+10  0.0 0.0
+11  0.0 0.0
+12  0.0 0.0
+13  0.0 0.0
+
+p 4 1
+0.7 0.7 0.7
+11  0.0 0.0
+10  0.0 0.0
+14  0.0 0.0
+15  0.0 0.0

graphics/ass10/main.c

+/* Computer Graphics, Assignment 10, Texture Mapping
+ *
+ * Filename ........ main.c
+ * Description ..... Creates OpenGL window and draws the scene.
+ * Created by ...... Paul Melis, Robert Belleman, Jurgen Sturm
+ *
+ * Student name ....
+ * Student email ...
+ * Collegekaart ....
+ * Date ............
+ * Comments ........
+ *
+ * (always fill in these fields before submitting!!)
+ */
+
+#if defined(__GNUC__)
+#define glCheckError(s) { GLenum error; if ((error=glGetError())) \
+  fprintf(stderr, "* GL error - %s (%s; %s, line %d): %s\n", \
+  (s), __FUNCTION__,__FILE__,__LINE__, gluErrorString(error)); }
+#else                           /* __GNUC__ */
+#define glCheckError(s) { GLenum error; if ((error=glGetError())) \
+  fprintf(stderr, "* GL error - %s (line %d): %s\n", \
+  (s), __LINE__, gluErrorString(error)); }
+#endif                          /* __GNUC__ */
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <GL/glut.h>
+#include <GL/gl.h>
+#include <GL/glu.h>
+#include <math.h>
+#include "polys.h"
+#include "ppmio.h"
+#include "geometry.h"
+
+struct texture_data
+{
+    const char *filename;
+    int         contains_transparency;
+};
+
+struct texture_data texture_files[] =
+{
+    // File name, transparency flag
+    { "textures/grass.ppm", 0 },
+    { "textures/roof2.ppm", 0 },
+    { "textures/wall3.ppm", 0 },
+    { "textures/wood.ppm", 0 },
+    { "textures/sky.ppm", 0 },
+    { "textures/road.ppm", 0 },
+    { "textures/banana_leaf2.ppm", 1 },
+    { NULL, 0 },
+};
+
+GLuint *texture_names;
+
+polys *polylistHouse = NULL;
+polys *polylistTreeStem = NULL;
+polys *polylistTreeLeafs = NULL;
+polys *polylistGround1 = NULL;
+polys *polylistGround2 = NULL;
+polys *polylistRoad = NULL;
+polys *polylistSkydome = NULL;
+
+int window;
+
+int show_textures = 0;
+int show_polygons_as_lines = 0;
+
+float camDistance = 24.0;
+float camRotZ = -45.0, camAzimuth = 20.0;
+float saved_camRotZ, saved_camAzimuth, saved_camDistance;
+int mouse_mode = 0;
+int mx, my;
+
+GLfloat mat_sun[] = { 1.0, 1.0, 0.5, 1.0 };
+GLfloat mat_no_sun[] = { 0.0, 0.0, 0.0, 1.0 };
+GLfloat mat_ambient[] = { 0.4, 0.4, 0.4, 1.0 };
+GLfloat mat_diffuse[] = { 1, 1, 1, 1.0 };
+GLfloat mat_specular[] = { 1.0, 1.0, 1.0, 1.0 };
+GLfloat mat_no_diffuse[] = { 1.0, 0.0, 0.0, 1.0 };
+GLfloat mat_no_specular[] = { 0.0, 0.0, 0.0, 1.0 };
+GLfloat mat_shininess[] = { 50.0 };
+GLfloat mat_no_shininess[] = { 0.0 };
+
+GLfloat light_position[] = { -8.0, 8.0, -8.0, 0.0 };
+GLfloat light_ambient[] = { 0.4, 0.4, 0.4, 0.0 };
+GLfloat light_diffuse[] = { 0.6, 0.6, 0.6, 0.0 };
+
+// Rough center points of the objects in the scene, used for
+// rotating around them with the camera
+
+#define NUM_OBJECTS 3
+vec3 object_positions[NUM_OBJECTS] =
+{
+    /* house, trees, skydome "center" */
+    {0.0, 0.0, 7.0},
+    {0.0, 0.0, 12.0},
+    {0.0, 10.0, 0.0},
+};
+
+// Determines around which object the camera rotates
+int current_object = 0;
+
+float rand_float(void);
+
+void
+InitializePolygonlists(void)
+{
+    int i;
+
+    // The road
+    polylistRoad = CreatePolylist(10);
+    loadPolygonalObject(polylistRoad, "road.obj", texture_names, 1.0, 0.0, 0.0, 0.0);
+
+    // We load the ground object twice, but use a different translation to put it once
+    // on either side of the road
+    polylistGround1 = CreatePolylist(10);
+    loadPolygonalObject(polylistGround1, "ground.obj", texture_names, 1.0, 0.0, 0.0, -55.0);
+    polylistGround2 = CreatePolylist(10);
+    loadPolygonalObject(polylistGround2, "ground.obj", texture_names, 1.0, 0.0, 0.0, 55.0);
+
+    // The skydome
+    polylistSkydome = CreatePolylist(1000);
+    createHemisphere(polylistSkydome, 100,  0, 0, 0,  1, 1, 1);
+    for (i = 0; i < polylistSkydome->length; i++)
+        polylistSkydome->items[i].texture_id = texture_names[4];
+
+    // Load the house and translate it, to put it next to the road
+    polylistHouse = CreatePolylist(10);
+    loadPolygonalObject(polylistHouse, "house.obj", texture_names, 1.0,
+        object_positions[0].x, object_positions[0].y, object_positions[0].z);
+
+    // A single tree object
+    polylistTreeLeafs = CreatePolylist(10);
+    createSphere(polylistTreeLeafs, 0.7, 0.7, 0.7,  0, 1.7, 0,  0, 1, 0);
+    for (i = 0; i < polylistTreeLeafs->length; i++)
+        polylistTreeLeafs->items[i].texture_id = texture_names[0];
+
+    polylistTreeStem = CreatePolylist(10);
+    createCylinder(polylistTreeStem, 0.075, 1.8,  0, 0, 0,  0.5, 0.3, 0);
+    for (i = 0; i < polylistTreeStem->length; i++)
+        polylistTreeStem->items[i].texture_id = texture_names[3];
+
+    // Done!
+    printf("%d polygons\n",
+           polylistHouse->length + polylistTreeLeafs->length +
+           polylistTreeStem->length + polylistGround1->length + polylistGround2->length);
+}
+
+
+void
+InitGL(void)
+{
+    GLvoid  *image_data;
+    int     i, width, height;
+
+    glMatrixMode(GL_PROJECTION);
+    glLoadIdentity();
+    gluPerspective(45.0f, (GLfloat) 640 / (GLfloat) 480, 0.1f, 1000.0f);
+
+    glClearColor(0.7, 0.7, 1, 1);
+    glDepthFunc(GL_LEQUAL);
+    glEnable(GL_DEPTH_TEST);
+    glCullFace(GL_BACK);
+
+    glEnable(GL_COLOR_MATERIAL);
+    glMaterialfv(GL_FRONT, GL_AMBIENT, mat_ambient);
+    glMaterialfv(GL_FRONT, GL_DIFFUSE, mat_diffuse);
+    glMaterialfv(GL_FRONT, GL_SPECULAR, mat_no_specular);
+    glMaterialfv(GL_FRONT, GL_SHININESS, mat_no_shininess);
+
+    glEnable(GL_LIGHTING);
+    glEnable(GL_LIGHT0);
+    glLightfv(GL_LIGHT0, GL_AMBIENT, light_ambient);
+    glLightfv(GL_LIGHT0, GL_DIFFUSE, light_diffuse);
+
+    glPixelStorei(GL_UNPACK_ALIGNMENT, 1);
+
+    // Using GL_REPLACE will effectively disable lighting.
+    // The default mode, GL_MODULATE, will use the texture color to
+    // module the shaded material color, so the end result still
+    // contains the shading.
+    //glTexEnvf(GL_TEXTURE_ENV, GL_TEXTURE_ENV_MODE, GL_REPLACE);
+
+    // generate texture objects
+
+    for (i = 0; texture_files[i].filename != NULL; i++)
+        ;
+    texture_names = malloc(i * sizeof(GLuint));
+    glGenTextures(i, texture_names);
+    glCheckError("glGenTextures");
+
+    // create textures from image files
+
+    GLint   texture_internal_format;
+    GLenum  texture_format, texture_type;
+
+    for (i = 0; texture_files[i].filename != NULL; i++)
+    {
+        image_data = (GLvoid *) readppm(texture_files[i].filename, &width, &height);
+        if (image_data)
+        {
+            printf("texture %d (OpenGL name = %d): %s, %d x %d\n", i,
+                texture_names[i], texture_files[i].filename, width, height);
+
+            if (texture_files[i].contains_transparency)
+            {
+                // Convert the RGB image to RGBA, by replacing
+                // red pixels with transparent ones
+
+                // Allocate a temporary RGBA image
+                unsigned char *temp_image = (unsigned char*) malloc(width*height*4);
+                unsigned char *input = image_data;
+                unsigned char *output = temp_image;
+
+                // Convert image to RGBA
+                for (int pix = 0; pix < width*height; pix++)
+                {
+                    if (*input == 255 && *(input+1) == 0 && *(input+2) == 0)
+                    {
+                        // Input pixel is pure red, make output pixel
+                        // white color and fully transparent
+                        *output++ = 255;
+                        *output++ = 255;
+                        *output++ = 255;
+                        *output++ = 0;      // alpha
+                    }
+                    else
+                    {
+                        // Copy image color, make fully opaque
+                        *output++ = *input;
+                        *output++ = *(input + 1);
+                        *output++ = *(input + 2);
+                        *output++ = 255;
+                    }
+
+                    input += 3;
+                }
+
+                // Replace the RGB image data with the generated RGBA data
+                free(image_data);
+                image_data = temp_image;
+
+                // This also influences some texture properties
+                texture_internal_format = GL_RGBA8;
+                texture_format = GL_RGBA;
+                texture_type = GL_UNSIGNED_BYTE;
+            }
+            else
+            {
+                texture_internal_format = GL_RGB8;
+                texture_format = GL_RGB;
+                texture_type = GL_UNSIGNED_BYTE;
+            }
+
+            glBindTexture(GL_TEXTURE_2D, texture_names[i]);
+            glCheckError("glBindTexture");
+
+            glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_WRAP_S, GL_CLAMP);
+            glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_WRAP_T, GL_CLAMP);
+            glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MAG_FILTER, GL_LINEAR);
+            glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MIN_FILTER, GL_LINEAR);
+            glCheckError("glTexParameteri");
+
+            glTexImage2D(GL_TEXTURE_2D, 0, texture_internal_format,
+                width, height, 0, texture_format, texture_type, image_data);
+            glCheckError("glTexImage2D");
+
+            // Free the image data, as OpenGL will have made its internal copy by now
+            free(image_data);
+        }
+        else
+        {
+            perror(texture_files[i].filename);
+            exit(0);
+        }
+    }
+}
+
+void
+ReSizeGLScene(int Width, int Height)
+{
+    if (Height == 0)
+        Height = 1;
+
+    glViewport(0, 0, Width, Height);
+
+    glMatrixMode(GL_PROJECTION);
+    glLoadIdentity();
+    gluPerspective(45.0f, (GLfloat) Width / (GLfloat) Height, 0.1f,
+                   1000.0f);
+
+    glMatrixMode(GL_MODELVIEW);
+    glLoadIdentity();
+}
+
+void
+DrawPolylist(polys * list)
+{
+    int i, j;
+
+    for (i = 0; i < list->length; i++)
+    {
+        poly p = list->items[i];
+
+        glColor3f(p.color[0], p.color[1], p.color[2]);
+
+        // Make the correct texture active
+        glBindTexture(GL_TEXTURE_2D, p.texture_id);
+
+        glBegin(GL_POLYGON);
+        for (j = 0; j < p.points; j++)
+        {
+            glNormal3f(p.normal[j].x, p.normal[j].y, p.normal[j].z);
+            glVertex3f(p.pts[j].x, p.pts[j].y, p.pts[j].z);
+        }
+        glEnd();
+    }
+}
+
+void
+SetupCamera(void)
+{
+    float   cx, cy, cz;
+    float   t;
+    float   beta, gamma;
+    vec3    op;
+
+    glMatrixMode(GL_MODELVIEW);
+    glLoadIdentity();
+
+    // Verbose, but straightforward way, of positioning the camera.
+    // Assume the camera's final position is (cx, cy, cz) and Y is up.
+    // Start with c being (camDistance, 0, 0)
+    // First rotate around Z, then around Y.
+    // Now we have c at the given distance from the origin, with specified rotation angles.
+
+    // degrees -> radians
+    beta = camAzimuth / 180.0 * 3.1415926535;
+    gamma = camRotZ / 180.0 * 3.1415926535;
+
+    cx = camDistance;
+    cy = cz = 0.0;
+
+    // Rotate around Z
+    t = cx;
+    cx = cx * cos(beta) + cy * sin(beta);
+    cy = t * sin(beta) + cy * cos(beta);
+
+    // Rotate around Y
+    t = cx;
+    cx = cx * cos(gamma) - cz * sin(gamma);
+    cz = t * sin(gamma) + cz * cos(gamma);
+
+    gluLookAt(cx, cy, cz, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0);
+
+    // Translate view to rotate around the selected object
+    op = object_positions[current_object];
+    glTranslatef(-op.x, -op.y, -op.z);
+}
+
+void
+DrawGLScene(void)
+{
+    float   tx, tz;
+
+    glClear(GL_DEPTH_BUFFER_BIT | GL_COLOR_BUFFER_BIT);
+
+    // set camera based on previous mouse input
+    SetupCamera();
+
+    // Set the light position. We do that after setting up the
+    // camera so that the light position will keep its world position
+    // and rotate along with the scene.
+    glLightfv(GL_LIGHT0, GL_POSITION, light_position);
+
+    // Draw the different objects
+
+    if (show_polygons_as_lines)
+    {
+        glPolygonMode(GL_FRONT_AND_BACK, GL_LINE);
+        glDisable(GL_TEXTURE_2D);
+    }
+    else
+    {
+        glPolygonMode(GL_FRONT_AND_BACK, GL_FILL);
+
+        if (show_textures)
+        {
+            glEnable(GL_ALPHA_TEST);
+            glAlphaFunc(GL_GREATER, 0.5);
+            glEnable(GL_TEXTURE_2D);
+        }
+        else
+            glDisable(GL_TEXTURE_2D);
+    }
+
+    DrawPolylist(polylistHouse);
+    DrawPolylist(polylistGround1);
+    DrawPolylist(polylistGround2);
+    DrawPolylist(polylistRoad);
+
+    // Draw a number of trees around the house.
+    // They are all identical, but differ in location and Y-rotation.
+
+    // Re-seed the random generator, so we always get the same sequence
+    // back from rand_float() below, for different runs of the program.
+    srand(95497452);
+
+    for (int t = 0; t < 12; t++)
+    {
+        glPushMatrix();
+
+        tx = 10 * (rand_float()-0.5) + object_positions[0].x;
+        tz = 3 * rand_float() + 2.0 + object_positions[0].z;
+        glTranslatef(tx, 0, tz);
+
+        glRotatef(rand_float()*360.0, 0, 1, 0);
+        glScalef(1, 1 + (rand_float()-0.5)*0.6, 1);
+
+        DrawPolylist(polylistTreeStem);
+        DrawPolylist(polylistTreeLeafs);
+
+        glPopMatrix();
+    }
+
+    // Draw the skydome with lighting turned off
+
+    glPushAttrib(GL_LIGHTING_BIT);
+    glDisable(GL_LIGHTING);
+    //DrawPolylist(polylistSkydome);
+    glPopAttrib();
+
+    glutSwapBuffers();
+}
+
+// Return a (pseudo-)random floating-point value in the range [0,1]
+float
+rand_float(void)
+{
+    return (float)(1.0 * rand() / RAND_MAX);
+}
+
+void
+keyPressed(unsigned char key, int x, int y)
+{
+    (void)x;
+    (void)y;
+
+    switch (key)
+    {
+        case 'o':
+            current_object = (current_object + 1) % NUM_OBJECTS;
+            glutPostRedisplay();
+            break;
+        case 't':
+            show_textures = 1 - show_textures;
+            glutPostRedisplay();
+            break;
+        case 'l':
+            show_polygons_as_lines = 1;
+            glutPostRedisplay();
+            break;
+        case 'p':
+            show_polygons_as_lines = 0;
+            glutPostRedisplay();
+            break;
+        case 'q':
+            exit(-1);
+    }
+}
+
+void
+specialKeyPressed(int key, int x, int y)
+{
+    (void)x;
+    (void)y;
+    (void)key;
+
+    DrawGLScene();
+}
+
+static void
+mouseFunc(int button, int state, int x, int y)
+{
+    // guard against both left and right buttons being pressed at the same time,
+    // by only responding when a mouse button is pressed while another one
+    // hasn't been pressed yet
+    if (state == GLUT_DOWN && mouse_mode == 0)
+    {
+        if (button == GLUT_LEFT_BUTTON)
+        {
+            mouse_mode = GLUT_LEFT_BUTTON;
+
+            saved_camRotZ = camRotZ;
+            saved_camAzimuth = camAzimuth;
+
+            mx = x;
+            my = y;
+        } else if (button == GLUT_RIGHT_BUTTON)
+        {
+            mouse_mode = GLUT_RIGHT_BUTTON;
+
+            saved_camDistance = camDistance;
+
+            my = y;
+        }
+    } else if (state == GLUT_UP && button == mouse_mode)
+    {
+        // pressed button released
+        mouse_mode = 0;
+    }
+
+}
+
+static void
+motionFunc(int x, int y)
+{
+    int dx, dy;
+
+    if (mouse_mode == GLUT_LEFT_BUTTON)
+    {
+        dx = mx - x;
+        dy = my - y;
+
+        camRotZ = saved_camRotZ - dx * 0.25;
+        camAzimuth = saved_camAzimuth - dy * 0.25;
+
+        if (camAzimuth > 89.99)
+            camAzimuth = 89.99;
+        else if (camAzimuth < -89.99)
+            camAzimuth = -89.99;
+    } else if (mouse_mode == GLUT_RIGHT_BUTTON)
+    {
+        dy = my - y;
+
+        camDistance = saved_camDistance - dy * 0.5;
+        if (camDistance < 0.1)
+            camDistance = 0.1;
+        else if (camDistance > 1000.0)
+            camDistance = 1000.0;
+    }
+}
+
+int
+main(int argc, char **argv)
+{
+    glutInit(&argc, argv);
+
+    glutInitDisplayMode(GLUT_DOUBLE | GLUT_RGBA | GLUT_ALPHA | GLUT_DEPTH);
+    glutInitWindowSize(640, 480);
+    glutInitWindowPosition(0, 0);
+    window = glutCreateWindow("OpenGL Framework");
+
+    glutDisplayFunc(&DrawGLScene);
+    glutIdleFunc(&DrawGLScene);
+    glutReshapeFunc(&ReSizeGLScene);
+    glutSpecialFunc(&specialKeyPressed);
+    glutKeyboardFunc(&keyPressed);
+    glutMouseFunc(&mouseFunc);
+    glutMotionFunc(&motionFunc);
+
+    InitGL();
+    InitializePolygonlists();
+
+    glutMainLoop();
+
+    return 1;
+}

graphics/ass10/polys.c

+/* Computer Graphics
+ *
+ * Filename ........ polys.c
+ * Description ..... Functions to manage lists of polygons
+ * Date ............ 19.08.2008
+ * Created by ...... Jurgen Sturm
+ * Cleaned up by ... Paul Melis
+ *
+ */
+
+#include <stdlib.h>
+#include <stdio.h>
+#include <errno.h>
+#include <string.h>
+#include "polys.h"
+
+/* Create a list of polys, and  immediately reserve space for `n' polys */
+
+polys*
+CreatePolylist(int n)
+{
+    polys *l;
+
+    l = (polys*) malloc(sizeof(polys));
+
+    if (l == NULL)
+    {
+        printf("CreatePolylist(): could not allocate memory (l): %s\n",
+           strerror(errno));
+        exit(-1);
+    }
+
+    if (n > 0)
+    {
+        l->items = (poly*) malloc(n * sizeof(poly));
+
+        if (l->items == NULL)
+        {
+            printf("CreatePolylist(): could not allocate memory (data): %s\n",
+                strerror(errno));
+            exit(-1);
+        }
+
+        l->capacity = n;
+        l->length = 0;
+    }
+    else
+    {
+        l->capacity = 0;
+        l->length = 0;
+        l->items = NULL;
+    }
+
+    return l;
+}
+
+
+void
+DestroyPolylist(polys *list)
+{
+    if (list == NULL)
+        return;
+
+    list->capacity = 0;
+    list->length = 0;
+    free(list->items);
+
+    list = NULL;
+}
+
+void
+ClearPolylist(polys *list)
+{
+    list->length = 0;
+}
+
+/* exit() if no memory available, returns 0 if poly* is NULL, 1 otherwise */
+static int
+GrowPolylist(polys *list, int delta)
+{
+    poly    *newdata;
+
+    if (list == NULL)
+        return 0;
+
+    newdata = (poly*) realloc(list->items, (list->capacity + delta)*sizeof(poly));
+
+    if (newdata == NULL)
+    {
+        printf("GrowPolylist(): could not allocate memory: %s\n",
+            strerror(errno));
+        exit(-1);
+    }
+
+    list->capacity += delta;
+    list->items = newdata;
+
+    return 1;
+}
+
+
+/* exit() if no memory available, returns 0 if poly* is NULL, 1 otherwise.
+ * Decreasing the data segment can, and will, destroy polygons if the new
+ * segment cannot contain all current data.
+ * The data segment cannot be reduced to a size less than 0 (size is clamped
+ * to zero).
+ * The new data segment will be a NULL pointer if the new size is 0. */
+static int
+ShrinkPolylist(polys *list, int delta)
+{
+    int     n;
+    poly    *newdata;
+
+    if (list == NULL)
+        return 0;
+
+    n = list->capacity - delta;
+    if (n < 1)
+    {
+        free(list->items);
+        list->items = NULL;
+        list->capacity = 0;
+        list->length = 0;
+
+        return 1;
+    }
+
+    newdata = (poly*)realloc(list->items, (list->capacity + delta)*sizeof(poly));
+
+    if (newdata == NULL)
+    {
+        printf("ShrinkPolylist(): could not allocate memory: %s\n",
+            strerror(errno));
+        exit(-1);
+    }
+
+    list->capacity -= delta;
+    list->items = newdata;
+
+    /* update list->length if neccesary */
+    if (list->length > list->capacity)
+        list->length = list->capacity;
+
+    return 1;
+}
+
+
+/* poly's are structs of (point) arrays, so they are passed by value instead
+ * of by reference (as do arrays) */
+void
+AddPolyToPolylist(polys *list, poly p)
+{
+    if (list == NULL)
+        return;
+
+    /* polylist is full, so first add some space */
+    if (list->length == list->capacity)
+    {
+        /* grow arbitrary amount */
+        if (GrowPolylist(list, 8) != 1)
+        {
+            printf("AddPolyToList(): failed");
+            exit(-1);
+        }
+    }
+
+    list->items[list->length] = p;
+    list->length++;
+}
+
+/* Append the items of 'to_append' to 'list' */
+void
+AppendPolylist(polys *list, polys *to_append)
+{
+    for (int i = 0; i < to_append->length; i++)
+        AddPolyToPolylist(list, to_append->items[i]);
+}
+
+polys*
+CopyPolylist(polys* list)
+{
+    int     i;
+    polys   *copy;
+
+    if (list == NULL)
+        return NULL;
+
+    copy = CreatePolylist(list->capacity);
+
+    /* not the most efficient, but certainly easiest & least error-prone */
+    for (i = 0; i < list->length; i++)
+        AddPolyToPolylist(copy, list->items[i]);
+
+    return copy;
+}
+
+

graphics/ass10/polys.h

+/* Computer Graphics
+ *
+ * Filename ........ polys.h
+ * Description ..... Functions to manage lists of polygons (header file)
+ * Date ............ 19.08.2008
+ * Created by ...... Jurgen Sturm
+ * Cleaned up by ... Paul Melis
+ */
+
+#ifndef POLYS_H
+#define POLYS_H
+
+#include <GL/gl.h>
+#include "v3math.h"
+
+#define MAX_VERTICES    4
+
+typedef struct poly
+{
+    // Number of vertices
+    int     points;
+
+    // Vertices and their normals
+    vec3    pts[MAX_VERTICES];
+    vec3    normal[MAX_VERTICES];
+
+    // Overall polygon color
+    GLfloat color[4];
+
+    // Texture identifier as set in a .obj file
+    GLuint  texture_id;
+
+    // Texture coordinates per vertex
+    vec3    tcoord[MAX_VERTICES];
+}
+poly;
+
+
+typedef struct polys
+{
+    /* size of the poly array "items" below */
+    int     capacity;
+
+    /* number of poly's stored, also index of the first free position
+     * in the "items" array. If this number is equal to "capacity" above,
+     * then this list is full */
+    int     length;
+
+    /* array of polygons, with length "capacity" of which only the first "length"
+     * items will be in use */
+    poly    *items;
+}
+polys;
+
+
+/* Create an empty list of polys, initially reserving space for 'n' polys */
+polys* CreatePolylist(int n);
+
+/* Destroy a list of polygons */
+void DestroyPolylist(polys *list);
+
+/* poly's are structs of (point) arrays, so they are passed by value instead
+ * of by reference (as do arrays) */
+void AddPolyToPolylist(polys *list, poly p);
+
+/* Append the items of 'to_append' to 'list' */
+void AppendPolylist(polys *list, polys *to_append);
+
+/* Copies the list of polys `list', calls AddPolytoPolylist() for each poly
+ * in the given list. */
+polys* CopyPolylist(polys *list);
+
+/* Clear a list of polygons. Note: this does not de-allocate the item array */
+void ClearPolylist(polys *list);
+
+#endif

graphics/ass10/ppmio.c

+/*
+  By Robert G. Belleman, Section Computational Science,
+  University of Amsterdam, Kruislaan 403, 1098 SJ Amsterdam,
+  the Netherlands.
+
+  Email : robbel@science.uva.nl
+  URL   : http://www.science.uva.nl/~robbel/
+  Tel.  : +31 20 525 7463 (secr)
+
+  $Id: ppmio.c,v 1.1.1.1 2001/10/09 11:43:52 robbel Exp $
+  $Log: ppmio.c,v $
+  Revision 1.1.1.1  2001/10/09 11:43:52  robbel
+  April 1997 CAVE demo developed by Robert G. Belleman in cooperation with
+  ESI France.
+
+
+*/
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+
+typedef unsigned char byte;
+
+/*
+  Read a ppm file. Very untested.
+  One thing this function can't stand is more than one line of comment after
+  the magic number (as xv tends to do).
+
+  Returns a [width][height][3] pointer containing the red, green, blue
+  components of a picture of size width x height. The area for this
+  picture has been allocated using alloc3() and should therefore be freed
+  by a call to free3(). The width and height of the picture are returned
+  in the passed pointers.
+
+  Returns NULL if anything failed.
+*/
+void *readppm(const char *filename, int *width, int *height)
+{
+    char buffer[80];            /* no line should be longer than 70 characters */
+    FILE *fd;
+    int maxcol;
+    byte *area;
+    int j;
+#if defined(DEBUG)
+    int dummy;
+#endif
+
+    if (!filename)
+        return NULL;
+
+    if (!(fd = fopen(filename, "r")))
+    {
+        perror(filename);
+        return NULL;
+    }
+
+    /* read and check magic number */
+    fscanf(fd, "%s\n", buffer);
+    if (strcmp(buffer, "P6"))
+        return NULL;            /* not a PPM raw file */
+
+    /* skip comment line (if any) */
+
+    fgets(buffer, 80, fd);
+    if (buffer[0] == '#')
+        fgets(buffer, 80, fd);
+
+    /* read width, height and maximum color-component value */
+    sscanf(buffer, "%d %d", width, height);
+    fscanf(fd, "%s\n", buffer);
+    maxcol = atoi(buffer);
+
+#if defined(DEBUG)
+    fprintf(stderr, "%s is %d x %d, max %d. Allocating %d bytes.\n",
+            filename, *width, *height, maxcol, (*width) * (*height) * 3);
+#endif
+
+    /* allocate storage, then read data */
+    if (!(area = (byte *) malloc((*width) * (*height) * 3 * sizeof(byte))))
+    {
+        fclose(fd);
+        return NULL;
+    }
+
+    /* From the glTexImage2D man page on the pixels parameter:
+
+    """The  first  element corresponds to the lower left corner of the texture image.  Subsequent elements progress left-to-right through the
+       remaining texels in the lowest row of the texture image, and then in successively higher rows of the texture image.  The final element
+       corresponds to the upper right corner of the texture image."""
+
+    So we read line-by-line here and store the data bottom-top swapped. This way
+    the data can be passed directly to glTexImage2D without getting a bottom-top swap.
+    */
+
+    for (j = *height-1; j >= 0; j--)
+    {
+#if defined(DEBUG)
+        dummy =
+#endif
+            fread(area + j*(*width)*3, sizeof(byte), (*width) * 3, fd);
+
+#if defined(DEBUG)
+        fprintf(stderr, "%d bytes read. First 16 bytes:\n", dummy);
+        for (dummy = 0; dummy < 16; dummy++)
+            fprintf(stderr, "%02x ", area[dummy]);
+        fprintf(stderr, "\n");
+#endif
+    }
+
+    fclose(fd);
+    return area;
+}

graphics/ass10/ppmio.h

+/*
+  By Robert G. Belleman, Section Computational Science,
+  University of Amsterdam, Kruislaan 403, 1098 SJ Amsterdam,
+  the Netherlands.
+
+  Email : robbel@science.uva.nl
+  URL   : http://www.science.uva.nl/~robbel/
+  Tel.  : +31 20 525 7463 (secr)
+
+  $Id: ppmio.h,v 1.1.1.1 2001/10/09 11:43:52 robbel Exp $
+  $Log: ppmio.h,v $
+  Revision 1.1.1.1  2001/10/09 11:43:52  robbel
+  April 1997 CAVE demo developed by Robert G. Belleman in cooperation with
+  ESI France.
+*/
+
+#ifndef _PPMIO_H
+#define _PPMIO_H
+
+void *readppm(const char *filename, int *width, int *height);
+
+#endif
+
+/* _PPMIO_H */

graphics/ass10/road.obj

+v 100 0 5
+v 100 0 -5
+v -100 0 -5
+v -100 0 5
+
+# (5=road texture)
+p 4 5
+0.85 0.85 0.85
+0  0.0 0.0
+1  0.0 0.0
+2  0.0 0.0
+3  0.0 0.0
\ No newline at end of file

graphics/ass10/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

linalg/prac3/freezeColors.m

+function freezeColors(varargin)
+% freezeColors  Lock colors of plot, enabling multiple colormaps per figure. (v2.3)
+%
+%   Problem: There is only one colormap per figure. This function provides
+%       an easy solution when plots using different colomaps are desired 
+%       in the same figure.
+%
+%   freezeColors freezes the colors of graphics objects in the current axis so 
+%       that subsequent changes to the colormap (or caxis) will not change the
+%       colors of these objects. freezeColors works on any graphics object 
+%       with CData in indexed-color mode: surfaces, images, scattergroups, 
+%       bargroups, patches, etc. It works by converting CData to true-color rgb
+%       based on the colormap active at the time freezeColors is called.
+%
+%   The original indexed color data is saved, and can be restored using
+%       unfreezeColors, making the plot once again subject to the colormap and
+%       caxis.
+%
+%
+%   Usage:
+%       freezeColors        applies to all objects in current axis (gca),
+%       freezeColors(axh)   same, but works on axis axh.
+%
+%   Example:
+%       subplot(2,1,1); imagesc(X); colormap hot; freezeColors
+%       subplot(2,1,2); imagesc(Y); colormap hsv; freezeColors etc...
+%
+%       Note: colorbars must also be frozen. Due to Matlab 'improvements' this can
+%				no longer be done with freezeColors. Instead, please
+%				use the function CBFREEZE by Carlos Adrian Vargas Aguilera
+%				that can be downloaded from the MATLAB File Exchange
+%				(http://www.mathworks.com/matlabcentral/fileexchange/24371)
+%
+%       h=colorbar; cbfreeze(h), or simply cbfreeze(colorbar)
+%
+%       For additional examples, see test/test_main.m
+%
+%   Side effect on render mode: freezeColors does not work with the painters
+%       renderer, because Matlab doesn't support rgb color data in
+%       painters mode. If the current renderer is painters, freezeColors
+%       changes it to zbuffer. This may have unexpected effects on other aspects
+%	      of your plots.
+%
+%       See also unfreezeColors, freezeColors_pub.html, cbfreeze.
+%
+%
+%   John Iversen (iversen@nsi.edu) 3/23/05
+%
+
+%   Changes:
+%   JRI (iversen@nsi.edu) 4/19/06   Correctly handles scaled integer cdata
+%   JRI 9/1/06   should now handle all objects with cdata: images, surfaces, 
+%                scatterplots. (v 2.1)
+%   JRI 11/11/06 Preserves NaN colors. Hidden option (v 2.2, not uploaded)
+%   JRI 3/17/07  Preserve caxis after freezing--maintains colorbar scale (v 2.3)
+%   JRI 4/12/07  Check for painters mode as Matlab doesn't support rgb in it.
+%   JRI 4/9/08   Fix preserving caxis for objects within hggroups (e.g. contourf)
+%   JRI 4/7/10   Change documentation for colorbars
+
+% Hidden option for NaN colors:
+%   Missing data are often represented by NaN in the indexed color
+%   data, which renders transparently. This transparency will be preserved
+%   when freezing colors. If instead you wish such gaps to be filled with 
+%   a real color, add 'nancolor',[r g b] to the end of the arguments. E.g. 
+%   freezeColors('nancolor',[r g b]) or freezeColors(axh,'nancolor',[r g b]),
+%   where [r g b] is a color vector. This works on images & pcolor, but not on
+%   surfaces.
+%   Thanks to Fabiano Busdraghi and Jody Klymak for the suggestions. Bugfixes 
+%   attributed in the code.
+
+% Free for all uses, but please retain the following:
+%   Original Author:
+%   John Iversen, 2005-10
+%   john_iversen@post.harvard.edu
+
+appdatacode = 'JRI__freezeColorsData';
+
+[h, nancolor] = checkArgs(varargin);
+
+%gather all children with scaled or indexed CData
+cdatah = getCDataHandles(h);
+
+%current colormap
+cmap = colormap;
+nColors = size(cmap,1);
+cax = caxis;
+
+% convert object color indexes into colormap to true-color data using 
+%  current colormap
+for hh = cdatah',
+    g = get(hh);
+    
+    %preserve parent axis clim
+    parentAx = getParentAxes(hh);
+    originalClim = get(parentAx, 'clim');    
+   
+    %   Note: Special handling of patches: For some reason, setting
+    %   cdata on patches created by bar() yields an error,
+    %   so instead we'll set facevertexcdata instead for patches.
+    if ~strcmp(g.Type,'patch'),
+        cdata = g.CData;
+    else
+        cdata = g.FaceVertexCData; 
+    end
+    
+    %get cdata mapping (most objects (except scattergroup) have it)
+    if isfield(g,'CDataMapping'),
+        scalemode = g.CDataMapping;
+    else
+        scalemode = 'scaled';
+    end
+    
+    %save original indexed data for use with unfreezeColors
+    siz = size(cdata);
+    setappdata(hh, appdatacode, {cdata scalemode});
+
+    %convert cdata to indexes into colormap
+    if strcmp(scalemode,'scaled'),
+        %4/19/06 JRI, Accommodate scaled display of integer cdata:
+        %       in MATLAB, uint * double = uint, so must coerce cdata to double
+        %       Thanks to O Yamashita for pointing this need out
+        idx = ceil( (double(cdata) - cax(1)) / (cax(2)-cax(1)) * nColors);
+    else %direct mapping
+        idx = cdata;
+        %10/8/09 in case direct data is non-int (e.g. image;freezeColors)
+        % (Floor mimics how matlab converts data into colormap index.)
+        % Thanks to D Armyr for the catch
+        idx = floor(idx);
+    end
+    
+    %clamp to [1, nColors]
+    idx(idx<1) = 1;
+    idx(idx>nColors) = nColors;
+
+    %handle nans in idx
+    nanmask = isnan(idx);
+    idx(nanmask)=1; %temporarily replace w/ a valid colormap index
+
+    %make true-color data--using current colormap
+    realcolor = zeros(siz);
+    for i = 1:3,
+        c = cmap(idx,i);
+        c = reshape(c,siz);
+        c(nanmask) = nancolor(i); %restore Nan (or nancolor if specified)
+        realcolor(:,:,i) = c;
+    end
+    
+    %apply new true-color color data
+    
+    %true-color is not supported in painters renderer, so switch out of that
+    if strcmp(get(gcf,'renderer'), 'painters'),
+        set(gcf,'renderer','zbuffer');
+    end
+    
+    %replace original CData with true-color data
+    if ~strcmp(g.Type,'patch'),
+        set(hh,'CData',realcolor);
+    else
+        set(hh,'faceVertexCData',permute(realcolor,[1 3 2]))
+    end
+    
+    %restore clim (so colorbar will show correct limits)
+    if ~isempty(parentAx),
+        set(parentAx,'clim',originalClim)
+    end
+    
+end %loop on indexed-color objects
+
+
+% ============================================================================ %
+% Local functions
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%% getCDataHandles -- get handles of all descendents with indexed CData
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+function hout = getCDataHandles(h)
+% getCDataHandles  Find all objects with indexed CData
+
+%recursively descend object tree, finding objects with indexed CData
+% An exception: don't include children of objects that themselves have CData:
+%   for example, scattergroups are non-standard hggroups, with CData. Changing
+%   such a group's CData automatically changes the CData of its children, 
+%   (as well as the children's handles), so there's no need to act on them.
+
+error(nargchk(1,1,nargin,'struct'))
+
+hout = [];
+if isempty(h),return;end
+
+ch = get(h,'children');
+for hh = ch'
+    g = get(hh);
+    if isfield(g,'CData'),     %does object have CData?
+        %is it indexed/scaled?
+        if ~isempty(g.CData) && isnumeric(g.CData) && size(g.CData,3)==1, 
+            hout = [hout; hh]; %#ok<AGROW> %yes, add to list
+        end
+    else %no CData, see if object has any interesting children
+            hout = [hout; getCDataHandles(hh)]; %#ok<AGROW>
+    end
+end
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%% getParentAxes -- return handle of axes object to which a given object belongs
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+function hAx = getParentAxes(h)
+% getParentAxes  Return enclosing axes of a given object (could be self)
+
+error(nargchk(1,1,nargin,'struct'))
+%object itself may be an axis
+if strcmp(get(h,'type'),'axes'),
+    hAx = h;
+    return
+end
+
+parent = get(h,'parent');
+if (strcmp(get(parent,'type'), 'axes')),
+    hAx = parent;
+else
+    hAx = getParentAxes(parent);
+end
+
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%% checkArgs -- Validate input arguments
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+function [h, nancolor] = checkArgs(args)
+% checkArgs  Validate input arguments to freezeColors
+
+nargs = length(args);
+error(nargchk(0,3,nargs,'struct'))
+
+%grab handle from first argument if we have an odd number of arguments
+if mod(nargs,2),
+    h = args{1};
+    if ~ishandle(h),
+        error('JRI:freezeColors:checkArgs:invalidHandle',...
+            'The first argument must be a valid graphics handle (to an axis)')
+    end
+    % 4/2010 check if object to be frozen is a colorbar
+    if strcmp(get(h,'Tag'),'Colorbar'),
+      if ~exist('cbfreeze.m'),
+        warning('JRI:freezeColors:checkArgs:cannotFreezeColorbar',...
+            ['You seem to be attempting to freeze a colorbar. This no longer'...
+            'works. Please read the help for freezeColors for the solution.'])
+      else
+        cbfreeze(h);
+        return
+      end
+    end
+    args{1} = [];
+    nargs = nargs-1;
+else
+    h = gca;
+end
+
+%set nancolor if that option was specified
+nancolor = [nan nan nan];
+if nargs == 2,
+    if strcmpi(args{end-1},'nancolor'),
+        nancolor = args{end};
+        if ~all(size(nancolor)==[1 3]),
+            error('JRI:freezeColors:checkArgs:badColorArgument',...
+                'nancolor must be [r g b] vector');
+        end
+        nancolor(nancolor>1) = 1; nancolor(nancolor<0) = 0;
+    else
+        error('JRI:freezeColors:checkArgs:unrecognizedOption',...
+            'Unrecognized option (%s). Only ''nancolor'' is valid.',args{end-1})
+    end
+end
+
+

linalg/prac3/leastSquares.m

+function x = leastSquares(A, b)
+	% x = inv(A'*inv(V)*A)*A'*inv(V)*b
+	% x = inv(A' * inv(V) * A) * A' * inv(V) *b
+	x = inv(A' * A) * A' * b
+end

linalg/prac3/planets.m

+% Load the observations
+load observations.mat
+
+% Define the dependent time variable
+Ax1 = [ones(20,1) t_obs t_obs .^2] ;
+Ay1 = Ax1 ;
+Ax2 = [ones(20,1) t_obs t_obs .^2 t_obs .^3] ;
+Ay2 = Ax2 ;
+Ax3 = [ones(20,1) sin(2*pi*t_obs/360) cos(2*pi*t_obs /360)] ;
+Ay3 = Ax3 ;
+
+% Define the independent coordinates variables
+bx1 = x_obs'; 
+by1 = y_obs';
+bx2 = x_obs';
+by2 = y_obs';
+bx3 = x_obs';
+by3 = y_obs';
+
+% Call your least squares function
+cx1 = leastSquares(Ax1, bx1)' ;
+cy1 = leastSquares(Ay1, by1)' ;
+cx2 = leastSquares(Ax2, bx2)' ;
+cy2 = leastSquares(Ay2, by2)' ;
+cx3 = leastSquares(Ax3, bx3)' ;
+cy3 = leastSquares(Ay3, by3)' ;
+
+% Run the test models functions to test your model
+test_models('f1', 'f1', cx1, cy1)
+test_models('f2', 'f2', cx2, cy2)
+test_models('f3', 'f3', cx3, cy3)
+

linalg/prac3/test_models.m

+function test_models(modelx, modely, coef_x, coef_y)
+
+load observations
+
+R = 30 ;
+t = linspace(0, 1, 360) ;
+x = R * cos(2 * pi * t)  ;
+y = R * sin(2 * pi * t) ;
+
+[x_s, y_s, z_s] = sphere ;
+x_s = 5 * x_s ;
+y_s = 5 * y_s ;
+z_s = 5 * z_s ;
+[x_e, y_e, z_e] = sphere ;
+
+tt = 1 : numel(x) ;
+if strcmp(modelx, 'f1')
+    x_est = coef_x(1) + coef_x(2) * tt + coef_x(3) * tt .^ 2 ;
+elseif strcmp(modelx, 'f2')
+    x_est = coef_x(1) + coef_x(2) * tt + coef_x(3) * tt .^ 2 + coef_x(4) * tt .^ 3 ;
+elseif strcmp(modelx, 'f3')
+    x_est = coef_x(1) + coef_x(2) * sin(2*pi*tt/360) + coef_x(3) * cos(2*pi*tt /360) ;
+end
+
+if strcmp(modely, 'f1')
+    y_est = coef_y(1) + coef_y(2) * tt + coef_y(3) * tt .^ 2 ;
+elseif strcmp(modely, 'f2')
+    y_est = coef_y(1) + coef_y(2) * tt + coef_y(3) * tt .^ 2 + coef_y(4) * tt .^ 3 ;
+elseif strcmp(modely, 'f3')
+    y_est = coef_y(1) + coef_y(2) * sin(2*pi*tt/360) + coef_y(3) * cos(2*pi*tt /360) ;
+end
+
+figure
+subplot(2, 2, 3)
+plot(t_obs, x_obs, 'r*')
+hold on
+plot(tt, x_est, 'g-') ;
+plot(tt, x, 'b-')
+legend({'Observations', 'Estimated values', 'True values'}) ;
+title(['LS fitting on x-observations using the ' modelx ' model']) ;
+grid on
+subplot(2, 2, 4)
+plot(t_obs, y_obs, 'r*')
+hold on
+tt = 1 : numel(y_est) ;
+plot(tt, y_est, 'g-') ;
+plot(tt, y, 'b-')
+grid on
+legend({'Observations', 'Estimated values', 'True values'}) ;
+title(['LS fitting on y-observations using the ' modely ' model']) ;
+
+h = subplot(2, 2, [1, 2]) ;
+hold on
+for i = 1 : numel(t)
+    cla(h)
+    plot3(x_obs, y_obs, zeros(numel(x_obs), 1), 'r*', 'MarkerSize', 12) ;
+    hold on
+    surf(x_s, y_s, z_s, 'LineStyle', 'none')  % sphere centered at (3,-2,0)
+    colormap('hot')
+    freezeColors
+    hold on
+    surf(x_e + x(i), y(i) + y_e, z_e, 'LineStyle', 'none')  % sphere centered at (3,-2,0)
+    colormap('winter')
+    plot3(x(1 : i), y(1 : i), zeros(i, 1), 'b-.') ;
+    plot3(x_est(1 : i), y_est(1 : i), zeros(i, 1), 'g-', 'Marker', '+', 'MarkerSize', 12) ;
+    axis([-R, R, -R, R, -5, 5]) ;
+    view([90, 90, 60])
+    set(gca, 'DataAspectRatio', [1, 1, 1]) ;
+    grid on
+%     legend({'Observations', 'Estimated orbit', 'True orbit;'}) ;
+    title('Pandora orbitting around its sun') ;
+    pause(0.003) ;
+%     daspect([20 20 20])
+end

linalg/prac3/unfreezeColors.m

+function unfreezeColors(h)
+% unfreezeColors  Restore colors of a plot to original indexed color. (v2.3)
+%
+%   Useful if you want to apply a new colormap to plots whose
+%       colors were previously frozen with freezeColors.
+%
+%   Usage:
+%       unfreezeColors          unfreezes all objects in current axis, 
+%       unfreezeColors(axh)     same, but works on axis axh. axh can be vector.
+%       unfreezeColors(figh)    same, but for all objects in figure figh.
+%
+%       Has no effect on objects on which freezeColors was not already called.
+%				(Note: if colorbars were frozen using cbfreeze, use cbfreeze('off') to 
+%       unfreeze them. See freezeColors for information on cbfreeze.)
+%
+%
+%   See also freezeColors, freezeColors_pub.html, cbfreeze.
+%
+%
+%   John Iversen (iversen@nsi.edu) 3/23/05
+%
+
+%   Changes:
+%   JRI 9/1/06 now restores any object with frozen CData;
+%              can unfreeze an entire figure at once.
+%   JRI 4/7/10 Change documentation for colorbars
+
+% Free for all uses, but please retain the following:
+%
+%   Original Author:
+%   John Iversen, 2005-10
+%   john_iversen@post.harvard.edu
+
+error(nargchk(0,1,nargin,'struct'))
+
+appdatacode = 'JRI__freezeColorsData';
+
+%default: operate on gca
+if nargin < 1,
+    h = gca;
+end
+
+if ~ishandle(h),
+     error('JRI:unfreezeColors:invalidHandle',...
+            'The argument must be a valid graphics handle to a figure or axis')
+end
+
+%if h is a figure, loop on its axes
+if strcmp(get(h,'type'),'figure'),
+    h = get(h,'children');
+end
+
+for h1 = h', %loop on axes
+
+    %process all children, acting only on those with saved CData
+    %   ( in appdata JRI__freezeColorsData)
+    ch = findobj(h1);
+    
+    for hh = ch',
+        
+        %some object handles may be invalidated when their parent changes 
+        %   (e.g. restoring colors of a scattergroup unfortunately changes
+        %   the handles of all its children). So, first check to make sure
+        %   it's a valid handle
+        if ishandle(hh)
+            if isappdata(hh,appdatacode),
+                ad = getappdata(hh,appdatacode);
+                %get oroginal cdata
+                %patches have to be handled separately (see note in freezeColors)
+                if ~strcmp(get(hh,'type'),'patch'),
+                    cdata = get(hh,'CData');
+                else
+                    cdata = get(hh,'faceVertexCData');
+                    cdata = permute(cdata,[1 3 2]);
+                end
+                indexed = ad{1};
+                scalemode = ad{2};
+                
+                %size consistency check
+                if all(size(indexed) == size(cdata(:,:,1))),
+                    %ok, restore indexed cdata
+                    if ~strcmp(get(hh,'type'),'patch'),
+                        set(hh,'CData',indexed);
+                    else
+                        set(hh,'faceVertexCData',indexed);
+                    end
+                    %restore cdatamapping, if needed
+                    g = get(hh);
+                    if isfield(g,'CDataMapping'),
+                        set(hh,'CDataMapping',scalemode);
+                    end
+                    %clear appdata
+                    rmappdata(hh,appdatacode)
+                else
+                    warning('JRI:unfreezeColors:internalCdataInconsistency',...
+                        ['Could not restore indexed data: it is the wrong size. ' ...
+                        'Were the axis contents changed since the call to freezeColors?'])
+                end
+                
+            end %test if has our appdata
+        end %test ishandle
+
+    end %loop on children
+
+end %loop on axes
+