diff --git a/examples/cpp/lib/svl/Basics.h b/examples/cpp/lib/svl/Basics.h
new file mode 100644
index 0000000000000000000000000000000000000000..c59852c4e3502300c84ecb97470ab7488d0c6cff
--- /dev/null
+++ b/examples/cpp/lib/svl/Basics.h
@@ -0,0 +1,162 @@
+/*
+ File: Basics.h
+
+ Function: Basic definitions for all files. Contains type
+ definitions, assertion and debugging facilities, and
+ miscellaneous useful template functions.
+
+ This is a cut-down version for SVL.
+
+ Author(s): Andrew Willmott
+
+ Copyright: (c) 1995-2001, Andrew Willmott
+
+ Notes: This header is affected by the following defines:
+
+ VL_CHECKING - Include code for assertions,
+ range errors and warnings.
+ VL_FLOAT - Use floats for real numbers. (Doubles
+ are the default.)
+ VL_NO_BOOL - There is no bool type.
+ VL_NO_TF - true and false are not predefined.
+*/
+
+#ifndef __Basics__
+#define __Basics__
+
+#include "SVLConfig.h"
+#include <iostream>
+#include <cmath>
+
+
+// --- Basic types -------------------------------------------------------------
+
+typedef void Void;
+typedef float Float;
+typedef double Double;
+typedef char Char;
+typedef int Short;
+typedef int Int;
+typedef long Long;
+typedef unsigned char Byte;
+typedef unsigned int UInt;
+
+#ifndef VL_FLOAT
+typedef Double Real;
+#else
+typedef Float Real;
+#endif
+
+#define SELF (*this) // A syntactic convenience.
+
+
+// --- Boolean type ------------------------------------------------------------
+
+// X11 #defines 'Bool' -- typical.
+
+#ifdef Bool
+#undef Bool
+#endif
+
+#ifndef VL_NO_BOOL
+// if the compiler implements the bool type...
+typedef bool Bool;
+#else
+// if not, make up our own.
+class Bool
+{
+public:
+
+ Bool() : val(0) {};
+ Bool(Int b) : val(b) {};
+
+ operator Int() { return val; };
+
+private:
+ Int val;
+};
+#ifdef VL_NO_TF
+enum {false, true};
+#endif
+#endif
+
+
+// --- Assertions and Range checking -------------------------------------------
+
+#define _Error(e) _Assert(false, e, __FILE__, __LINE__)
+#define _Warning(w) _Expect(false, w, __FILE__, __LINE__)
+
+#if defined(VL_CHECKING)
+#define Assert(b, e) _Assert(b, e, __FILE__, __LINE__)
+ // Assert that b is true. e is an error message to be printed if b
+ // is false.
+#define Expect(b, w) _Expect(b, w, __FILE__, __LINE__)
+ // Prints warning w if b is false
+#define CheckRange(i, l, u, r) _CheckRange(i, l, u, r, __FILE__, __LINE__)
+ // Checks whether i is in the range [lowerBound, upperBound).
+#else
+#define Assert(b, e)
+#define Expect(b, w)
+#define CheckRange(a, l, u, r)
+#endif
+
+Void _Assert(Int condition, const Char *errorMessage, const Char *file, Int line);
+Void _Expect(Int condition, const Char *warningMessage, const Char *file, Int line);
+Void _CheckRange(Int i, Int lowerBound, Int upperBound, const Char *rangeMessage,
+ const Char *file, Int line);
+
+
+// --- Inlines -----------------------------------------------------------------
+
+template<class Value>
+ inline Value Min(Value x, Value y)
+ {
+ if (x <= y)
+ return(x);
+ else
+ return(y);
+ };
+
+template<class Value>
+ inline Value Max(Value x, Value y)
+ {
+ if (x >= y)
+ return(x);
+ else
+ return(y);
+ };
+
+template<class Value>
+ inline Void Swap(Value &x, Value &y)
+ {
+ Value t;
+
+ t = x;
+ x = y;
+ y = t;
+ };
+
+template<class Value>
+ inline Value Mix(Value x, Value y, Real s)
+ {
+ return(x + (y - x) * s);
+ };
+
+template<class Value>
+ inline Value Clip(Value x, Value min, Value max)
+ {
+ if (x < min)
+ return(min);
+ else if (x > max)
+ return(max);
+ else
+ return(x);
+ };
+
+template<class Value>
+ inline Value sqr(Value x)
+ {
+ return(x * x);
+ };
+
+#endif
diff --git a/examples/cpp/lib/svl/Constants.h b/examples/cpp/lib/svl/Constants.h
new file mode 100644
index 0000000000000000000000000000000000000000..bc3a5b8875fb5a65b0946e1b6894c157159b4471
--- /dev/null
+++ b/examples/cpp/lib/svl/Constants.h
@@ -0,0 +1,43 @@
+/*
+ File: Constants.h
+
+ Function: Contains various constants for VL.
+
+ Author: Andrew Willmott
+
+ Copyright: (c) 1999-2001, Andrew Willmott
+*/
+
+#ifndef __VLConstants__
+#define __VLConstants__
+
+#include <cmath>
+
+
+// --- Mathematical constants -------------------------------------------------
+
+
+#ifdef M_PI
+const Real vl_pi = M_PI;
+const Real vl_halfPi = M_PI_2;
+#elif defined(_PI)
+const Real vl_pi = _PI;
+const Real vl_halfPi = vl_pi / 2.0;
+#else
+const Real vl_pi = 3.14159265358979323846;
+const Real vl_halfPi = vl_pi / 2.0;
+#endif
+
+#ifdef HUGE_VAL
+const Double vl_inf = HUGE_VAL;
+#endif
+
+enum ZeroOrOne { vl_zero = 0, vl_0 = 0, vl_one = 1, vl_I = 1, vl_1 = 1 };
+enum Block { vl_Z = 0, vl_B = 1, vl_block = 1 };
+enum Axis { vl_x, vl_y, vl_z, vl_w };
+typedef Axis vl_axis;
+
+const UInt VL_REF_FLAG = UInt(1) << (sizeof(UInt) * 8 - 1);
+const UInt VL_REF_MASK = (~VL_REF_FLAG);
+
+#endif
diff --git a/examples/cpp/lib/svl/LICENSE b/examples/cpp/lib/svl/LICENSE
new file mode 100644
index 0000000000000000000000000000000000000000..7c1e717dab0c6ea10c7aa55a848865312bd6fd87
--- /dev/null
+++ b/examples/cpp/lib/svl/LICENSE
@@ -0,0 +1,23 @@
+
+Except where otherwise noted, Copyright (c) 2002 Andrew Willmott.
+All rights reserved.
+
+Redistribution and use in source and binary forms, with or without
+modification, are permitted provided that the following conditions
+are met:
+
+1. Redistributions of source code must retain the above copyright
+ notice, this list of conditions and the following disclaimer.
+
+2. Redistributions in binary form have no conditions on them.
+
+THIS SOFTWARE IS PROVIDED ``AS IS'' AND ANY EXPRESSED OR IMPLIED
+WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF
+MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN
+NO EVENT SHALL THE AUTHORS BE LIABLE FOR ANY DIRECT, INDIRECT,
+INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT
+NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF
+USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON
+ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
+(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF
+THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
diff --git a/examples/cpp/lib/svl/Mat.h b/examples/cpp/lib/svl/Mat.h
new file mode 100644
index 0000000000000000000000000000000000000000..8162ed892c3f9caa148e771ee9bc400bb54621a0
--- /dev/null
+++ b/examples/cpp/lib/svl/Mat.h
@@ -0,0 +1,213 @@
+/*
+ File: Mat.h
+
+ Function: Defines a generic resizeable matrix.
+
+ Author(s): Andrew Willmott
+
+ Copyright: (c) 1995-2001, Andrew Willmott
+ */
+
+#ifndef __Mat__
+#define __Mat__
+
+#include "svl/Vec.h"
+
+class Mat2;
+class Mat3;
+class Mat4;
+
+
+// --- Mat Class --------------------------------------------------------------
+
+class Mat
+{
+public:
+
+ // Constructors
+
+ Mat(); // Null matrix
+ Mat(Int rows, Int cols); // Uninitialised matrix
+ Mat(Int rows, Int cols, Double elt0 ...); // Mat(2, 2, 1.0, 2.0, 3.0, 4.0)
+ Mat(Int nrows, Int ncols, Real *ndata); // Create reference matrix
+ Mat(const Mat &m); // Copy constructor
+ Mat(const Mat2 &m); // reference to a Mat2
+ Mat(const Mat3 &m); // reference to a Mat3
+ Mat(const Mat4 &m); // reference to a Mat4
+ Mat(Int rows, Int cols, ZeroOrOne k); // I * k
+ Mat(Int rows, Int cols, Block k); // block matrix (m[i][j] = k)
+
+ ~Mat();
+
+ // Accessor methods
+
+ Int Rows() const { return(rows & VL_REF_MASK); };
+ Int Cols() const { return(cols); };
+
+ Vec operator [] (Int i); // Indexing by row
+ Vec operator [] (Int i) const; // Indexing by row
+
+ Real &Elt(Int i, Int j); // Indexing by elt
+ Real Elt(Int i, Int j) const;
+
+ Void SetSize(Int nrows, Int ncols);
+ Void SetSize(const Mat &m);
+ Bool IsSquare() const
+ { return((rows & VL_REF_MASK) == cols); };
+
+ Real *Ref() const; // Return pointer to data
+
+ // Assignment operators
+
+ Mat &operator = (const Mat &m); // Assignment of a matrix
+ Mat &operator = (ZeroOrOne k); // Set to k * I...
+ Mat &operator = (Block k); // Set to a block matrix...
+ Mat &operator = (const Mat2 &m);
+ Mat &operator = (const Mat3 &m);
+ Mat &operator = (const Mat4 &m);
+
+ // In-Place Operators
+
+ Mat &operator += (const Mat &m);
+ Mat &operator -= (const Mat &m);
+ Mat &operator *= (const Mat &m);
+ Mat &operator *= (Real s);
+ Mat &operator /= (Real s);
+
+ // Matrix initialisers
+
+ Void MakeZero();
+ Void MakeDiag(Real k);
+ Void MakeDiag();
+ Void MakeBlock(Real k);
+ Void MakeBlock();
+
+ Mat &Clamp(Real fuzz);
+ Mat &Clamp();
+
+ // Private...
+
+protected:
+
+ Real *data;
+ UInt rows;
+ UInt cols;
+
+ Bool IsRef() { return((rows & VL_REF_FLAG) != 0); };
+};
+
+
+// --- Mat Comparison Operators -----------------------------------------------
+
+Bool operator == (const Mat &m, const Mat &n);
+Bool operator != (const Mat &m, const Mat &n);
+
+
+// --- Mat Arithmetic Operators -----------------------------------------------
+
+Mat operator + (const Mat &m, const Mat &n);
+Mat operator - (const Mat &m, const Mat &n);
+Mat operator - (const Mat &m);
+Mat operator * (const Mat &m, const Mat &n);
+Mat operator * (const Mat &m, Real s);
+inline Mat operator * (Real s, const Mat &m);
+Mat operator / (const Mat &m, Real s);
+
+Vec operator * (const Mat &m, const Vec &v);
+Vec operator * (const Vec &v, const Mat &m);
+
+Mat trans(const Mat &m); // Transpose
+Real trace(const Mat &m); // Trace
+Mat inv(const Mat &m, Real *determinant = 0, Real pEps = 1e-20);
+ // Inverse
+Mat oprod(const Vec &a, const Vec &b); // Outer product
+
+Mat clamped(const Mat &m, Real fuzz);
+Mat clamped(const Mat &m);
+
+
+// --- Mat Input & Output -----------------------------------------------------
+
+std::ostream &operator << (std::ostream &s, const Mat &m);
+std::istream &operator >> (std::istream &s, Mat &m);
+
+
+// --- Mat Inlines ------------------------------------------------------------
+
+inline Mat::Mat() : data(0), rows(0), cols(0)
+{
+}
+
+inline Mat::Mat(Int rows, Int cols) : rows(rows), cols(cols)
+{
+ Assert(rows > 0 && cols > 0, "(Mat) illegal matrix size");
+
+ data = new Real[rows * cols];
+}
+
+inline Mat::Mat(Int nrows, Int ncols, Real *ndata) :
+ data(ndata), rows(nrows | VL_REF_FLAG), cols(ncols)
+{
+}
+
+inline Vec Mat::operator [] (Int i)
+{
+ CheckRange(i, 0, Rows(), "(Mat::[i]) i index out of range");
+
+ return(Vec(cols, data + i * cols));
+}
+
+inline Vec Mat::operator [] (Int i) const
+{
+ CheckRange(i, 0, Rows(), "(Mat::[i]) i index out of range");
+
+ return(Vec(cols, data + i * cols));
+}
+
+inline Real &Mat::Elt(Int i, Int j)
+{
+ CheckRange(i, 0, Rows(), "(Mat::e(i,j)) i index out of range");
+ CheckRange(j, 0, Cols(), "(Mat::e(i,j)) j index out of range");
+
+ return(data[i * cols + j]);
+}
+
+inline Real Mat::Elt(Int i, Int j) const
+{
+ CheckRange(i, 0, Rows(), "(Mat::e(i,j)) i index out of range");
+ CheckRange(j, 0, Cols(), "(Mat::e(i,j)) j index out of range");
+
+ return(data[i * cols + j]);
+}
+
+inline Real *Mat::Ref() const
+{
+ return(data);
+}
+
+inline Mat operator * (Real s, const Mat &m)
+{
+ return(m * s);
+}
+
+inline Mat &Mat::operator = (ZeroOrOne k)
+{
+ MakeDiag(k);
+
+ return(SELF);
+}
+
+inline Mat &Mat::operator = (Block k)
+{
+ MakeBlock((ZeroOrOne) k);
+
+ return(SELF);
+}
+
+inline Mat::~Mat()
+{
+ if (!IsRef())
+ delete[] data;
+}
+
+#endif
diff --git a/examples/cpp/lib/svl/Mat2.h b/examples/cpp/lib/svl/Mat2.h
new file mode 100644
index 0000000000000000000000000000000000000000..7f0a6fa6ef902bb44eecb2351cd0f1962b5cc000
--- /dev/null
+++ b/examples/cpp/lib/svl/Mat2.h
@@ -0,0 +1,379 @@
+/*
+ File: Mat2.h
+
+ Function: Defines a 2 x 2 matrix.
+
+ Author(s): Andrew Willmott
+
+ Copyright: (c) 1995-2001, Andrew Willmott
+ */
+
+#ifndef __Mat2__
+#define __Mat2__
+
+#include "svl/Vec2.h"
+
+
+// --- Mat2 Class -------------------------------------------------------------
+
+class Mat2
+{
+public:
+
+ // Constructors
+
+ Mat2();
+ Mat2(Real a, Real b, Real c, Real d); // Create from rows
+ Mat2(const Mat2 &m); // Copy constructor
+ Mat2(ZeroOrOne k);
+ Mat2(Block k);
+
+ // Accessor functions
+
+ Int Rows() const { return(2); };
+ Int Cols() const { return(2); };
+
+ Vec2 &operator [] (Int i);
+ const Vec2 &operator [] (Int i) const;
+
+ Real *Ref() const; // Return pointer to data
+
+ // Assignment operators
+
+ Mat2 &operator = (const Mat2 &m);
+ Mat2 &operator = (ZeroOrOne k);
+ Mat2 &operator = (Block k);
+ Mat2 &operator += (const Mat2 &m);
+ Mat2 &operator -= (const Mat2 &m);
+ Mat2 &operator *= (const Mat2 &m);
+ Mat2 &operator *= (Real s);
+ Mat2 &operator /= (Real s);
+
+ // Comparison operators
+
+ Bool operator == (const Mat2 &m) const; // M == N?
+ Bool operator != (const Mat2 &m) const; // M != N?
+
+ // Arithmetic operators
+
+ Mat2 operator + (const Mat2 &m) const; // M + N
+ Mat2 operator - (const Mat2 &m) const; // M - N
+ Mat2 operator - () const; // -M
+ Mat2 operator * (const Mat2 &m) const; // M * N
+ Mat2 operator * (Real s) const; // M * s
+ Mat2 operator / (Real s) const; // M / s
+
+ // Initialisers
+
+ Void MakeZero(); // Zero matrix
+ Void MakeDiag(Real k = vl_one); // I
+ Void MakeBlock(Real k = vl_one); // all elts=k
+
+ // Vector Transformations
+
+ Mat2& MakeRot(Real theta);
+ Mat2& MakeScale(const Vec2 &s);
+
+ // Private...
+
+protected:
+
+ Vec2 row[2]; // Rows of the matrix
+};
+
+
+// --- Matrix operators -------------------------------------------------------
+
+inline Vec2 &operator *= (Vec2 &v, const Mat2 &m); // v *= m
+inline Vec2 operator * (const Mat2 &m, const Vec2 &v); // m * v
+inline Vec2 operator * (const Vec2 &v, const Mat2 &m); // v * m
+inline Mat2 operator * (Real s, const Mat2 &m); // s * m
+
+inline Mat2 trans(const Mat2 &m); // Transpose
+inline Real trace(const Mat2 &m); // Trace
+inline Mat2 adj(const Mat2 &m); // Adjoint
+Real det(const Mat2 &m); // Determinant
+Mat2 inv(const Mat2 &m); // Inverse
+Mat2 oprod(const Vec2 &a, const Vec2 &b);
+ // Outer product
+
+// The xform functions help avoid dependence on whether row or column
+// vectors are used to represent points and vectors.
+inline Vec2 xform(const Mat2 &m, const Vec2 &v); // Transform of v by m
+inline Mat2 xform(const Mat2 &m, const Mat2 &n); // xform v -> m(n(v))
+
+std::ostream &operator << (std::ostream &s, const Mat2 &m);
+std::istream &operator >> (std::istream &s, Mat2 &m);
+
+
+// --- Inlines ----------------------------------------------------------------
+
+inline Vec2 &Mat2::operator [] (Int i)
+{
+ CheckRange(i, 0, 2, "(Mat2::[i]) index out of range");
+ return(row[i]);
+}
+
+inline const Vec2 &Mat2::operator [] (Int i) const
+{
+ CheckRange(i, 0, 2, "(Mat2::[i]) index out of range");
+ return(row[i]);
+}
+
+inline Real *Mat2::Ref() const
+{
+ return((Real*) row);
+}
+
+inline Mat2::Mat2()
+{
+}
+
+inline Mat2::Mat2(Real a, Real b, Real c, Real d)
+{
+ row[0][0] = a; row[0][1] = b;
+ row[1][0] = c; row[1][1] = d;
+}
+
+inline Mat2::Mat2(const Mat2 &m)
+{
+ row[0] = m[0];
+ row[1] = m[1];
+}
+
+
+inline Void Mat2::MakeZero()
+{
+ row[0][0] = vl_zero; row[0][1] = vl_zero;
+ row[1][0] = vl_zero; row[1][1] = vl_zero;
+}
+
+inline Void Mat2::MakeDiag(Real k)
+{
+ row[0][0] = k; row[0][1] = vl_zero;
+ row[1][0] = vl_zero; row[1][1] = k;
+}
+
+inline Void Mat2::MakeBlock(Real k)
+{
+ row[0][0] = k; row[0][1] = k;
+ row[1][0] = k; row[1][1] = k;
+}
+
+inline Mat2::Mat2(ZeroOrOne k)
+{
+ MakeDiag(k);
+}
+
+inline Mat2::Mat2(Block k)
+{
+ MakeBlock((ZeroOrOne) k);
+}
+
+inline Mat2 &Mat2::operator = (ZeroOrOne k)
+{
+ MakeDiag(k);
+
+ return(SELF);
+}
+
+inline Mat2 &Mat2::operator = (Block k)
+{
+ MakeBlock((ZeroOrOne) k);
+
+ return(SELF);
+}
+
+inline Mat2 &Mat2::operator = (const Mat2 &m)
+{
+ row[0] = m[0];
+ row[1] = m[1];
+
+ return(SELF);
+}
+
+inline Mat2 &Mat2::operator += (const Mat2 &m)
+{
+ row[0] += m[0];
+ row[1] += m[1];
+
+ return(SELF);
+}
+
+inline Mat2 &Mat2::operator -= (const Mat2 &m)
+{
+ row[0] -= m[0];
+ row[1] -= m[1];
+
+ return(SELF);
+}
+
+inline Mat2 &Mat2::operator *= (const Mat2 &m)
+{
+ SELF = SELF * m;
+
+ return(SELF);
+}
+
+inline Mat2 &Mat2::operator *= (Real s)
+{
+ row[0] *= s;
+ row[1] *= s;
+
+ return(SELF);
+}
+
+inline Mat2 &Mat2::operator /= (Real s)
+{
+ row[0] /= s;
+ row[1] /= s;
+
+ return(SELF);
+}
+
+
+inline Mat2 Mat2::operator + (const Mat2 &m) const
+{
+ Mat2 result;
+
+ result[0] = row[0] + m[0];
+ result[1] = row[1] + m[1];
+
+ return(result);
+}
+
+inline Mat2 Mat2::operator - (const Mat2 &m) const
+{
+ Mat2 result;
+
+ result[0] = row[0] - m[0];
+ result[1] = row[1] - m[1];
+
+ return(result);
+}
+
+inline Mat2 Mat2::operator - () const
+{
+ Mat2 result;
+
+ result[0] = -row[0];
+ result[1] = -row[1];
+
+ return(result);
+}
+
+inline Mat2 Mat2::operator * (const Mat2 &m) const
+{
+#define N(x,y) row[x][y]
+#define M(x,y) m.row[x][y]
+#define R(x,y) result[x][y]
+
+ Mat2 result;
+
+ R(0,0) = N(0,0) * M(0,0) + N(0,1) * M(1,0);
+ R(0,1) = N(0,0) * M(0,1) + N(0,1) * M(1,1);
+ R(1,0) = N(1,0) * M(0,0) + N(1,1) * M(1,0);
+ R(1,1) = N(1,0) * M(0,1) + N(1,1) * M(1,1);
+
+ return(result);
+
+#undef N
+#undef M
+#undef R
+}
+
+inline Mat2 Mat2::operator * (Real s) const
+{
+ Mat2 result;
+
+ result[0] = row[0] * s;
+ result[1] = row[1] * s;
+
+ return(result);
+}
+
+inline Mat2 Mat2::operator / (Real s) const
+{
+ Mat2 result;
+
+ result[0] = row[0] / s;
+ result[1] = row[1] / s;
+
+ return(result);
+}
+
+inline Mat2 operator * (Real s, const Mat2 &m)
+{
+ return(m * s);
+}
+
+inline Vec2 operator * (const Mat2 &m, const Vec2 &v)
+{
+ Vec2 result;
+
+ result[0] = m[0][0] * v[0] + m[0][1] * v[1];
+ result[1] = m[1][0] * v[0] + m[1][1] * v[1];
+
+ return(result);
+}
+
+inline Vec2 operator * (const Vec2 &v, const Mat2 &m)
+{
+ Vec2 result;
+
+ result[0] = v[0] * m[0][0] + v[1] * m[1][0];
+ result[1] = v[0] * m[0][1] + v[1] * m[1][1];
+
+ return(result);
+}
+
+inline Vec2 &operator *= (Vec2 &v, const Mat2 &m)
+{
+ Real t;
+
+ t = v[0] * m[0][0] + v[1] * m[1][0];
+ v[1] = v[0] * m[0][1] + v[1] * m[1][1];
+ v[0] = t;
+
+ return(v);
+}
+
+
+inline Mat2 trans(const Mat2 &m)
+{
+ Mat2 result;
+
+ result[0][0] = m[0][0]; result[0][1] = m[1][0];
+ result[1][0] = m[0][1]; result[1][1] = m[1][1];
+
+ return(result);
+}
+
+inline Real trace(const Mat2 &m)
+{
+ return(m[0][0] + m[1][1]);
+}
+
+inline Mat2 adj(const Mat2 &m)
+{
+ Mat2 result;
+
+ result[0] = cross(m[1]);
+ result[1] = -cross(m[0]);
+
+ return(result);
+}
+
+#ifdef VL_ROW_ORIENT
+inline Vec2 xform(const Mat2 &m, const Vec2 &v)
+{ return(v * m); }
+inline Mat2 xform(const Mat2 &m, const Mat2 &n)
+{ return(n * m); }
+#else
+inline Vec2 xform(const Mat2 &m, const Vec2 &v)
+{ return(m * v); }
+inline Mat2 xform(const Mat2 &m, const Mat2 &n)
+{ return(m * n); }
+#endif
+
+#endif
diff --git a/examples/cpp/lib/svl/Mat3.h b/examples/cpp/lib/svl/Mat3.h
new file mode 100644
index 0000000000000000000000000000000000000000..a528c7aa4bb300cf6c4062a3b50f06e9fc054120
--- /dev/null
+++ b/examples/cpp/lib/svl/Mat3.h
@@ -0,0 +1,224 @@
+/*
+ File: Mat3.h
+
+ Function: Defines a 3 x 3 matrix.
+
+ Author(s): Andrew Willmott
+
+ Copyright: (c) 1995-2001, Andrew Willmott
+*/
+
+#ifndef __Mat3__
+#define __Mat3__
+
+#include "svl/Vec3.h"
+
+
+// --- Mat3 Class -------------------------------------------------------------
+
+
+class Vec4;
+
+class Mat3
+{
+public:
+
+ // Constructors
+
+ Mat3();
+ Mat3(Real a, Real b, Real c,
+ Real d, Real e, Real f,
+ Real g, Real h, Real i);
+ Mat3(const Mat3 &m);
+ Mat3(ZeroOrOne k);
+ Mat3(Block k);
+
+ // Accessor functions
+
+ Int Rows() const { return(3); };
+ Int Cols() const { return(3); };
+
+ Vec3 &operator [] (Int i);
+ const Vec3 &operator [] (Int i) const;
+
+ Real *Ref() const; // Return pointer to data
+
+ // Assignment operators
+
+ Mat3 &operator = (const Mat3 &m);
+ Mat3 &operator = (ZeroOrOne k);
+ Mat3 &operator = (Block k);
+ Mat3 &operator += (const Mat3 &m);
+ Mat3 &operator -= (const Mat3 &m);
+ Mat3 &operator *= (const Mat3 &m);
+ Mat3 &operator *= (Real s);
+ Mat3 &operator /= (Real s);
+
+ // Comparison operators
+
+ Bool operator == (const Mat3 &m) const; // M == N?
+ Bool operator != (const Mat3 &m) const; // M != N?
+
+ // Arithmetic operators
+
+ Mat3 operator + (const Mat3 &m) const; // M + N
+ Mat3 operator - (const Mat3 &m) const; // M - N
+ Mat3 operator - () const; // -M
+ Mat3 operator * (const Mat3 &m) const; // M * N
+ Mat3 operator * (Real s) const; // M * s
+ Mat3 operator / (Real s) const; // M / s
+
+ // Initialisers
+
+ Void MakeZero(); // Zero matrix
+ Void MakeDiag(Real k = vl_one); // I
+ Void MakeBlock(Real k = vl_one); // all elts = k
+
+ // Vector Transforms
+
+ Mat3& MakeRot(const Vec3 &axis, Real theta);
+ Mat3& MakeRot(const Vec4 &q); // Rotate by quaternion
+ Mat3& MakeScale(const Vec3 &s);
+
+ // Homogeneous Transforms
+
+ Mat3& MakeHRot(Real theta); // Rotate by theta rads
+ Mat3& MakeHScale(const Vec2 &s); // Scale by s
+ Mat3& MakeHTrans(const Vec2 &t); // Translation by t
+
+ // Private...
+
+protected:
+
+ Vec3 row[3];
+};
+
+
+// --- Matrix operators -------------------------------------------------------
+
+inline Vec3 &operator *= (Vec3 &v, const Mat3 &m); // v *= m
+inline Vec3 operator * (const Mat3 &m, const Vec3 &v); // m * v
+inline Vec3 operator * (const Vec3 &v, const Mat3 &m); // v * m
+inline Mat3 operator * (const Real s, const Mat3 &m); // s * m
+
+Mat3 trans(const Mat3 &m); // Transpose
+Real trace(const Mat3 &m); // Trace
+Mat3 adj(const Mat3 &m); // Adjoint
+Real det(const Mat3 &m); // Determinant
+Mat3 inv(const Mat3 &m); // Inverse
+Mat3 oprod(const Vec3 &a, const Vec3 &b); // Outer product
+
+// The xform functions help avoid dependence on whether row or column
+// vectors are used to represent points and vectors.
+inline Vec3 xform(const Mat3 &m, const Vec3 &v); // Transform of v by m
+inline Vec2 xform(const Mat3 &m, const Vec2 &v); // Hom. xform of v by m
+inline Mat3 xform(const Mat3 &m, const Mat3 &n); // Xform v -> m(n(v))
+
+std::ostream &operator << (std::ostream &s, const Mat3 &m);
+std::istream &operator >> (std::istream &s, Mat3 &m);
+
+
+// --- Inlines ----------------------------------------------------------------
+
+inline Mat3::Mat3()
+{
+}
+
+inline Vec3 &Mat3::operator [] (Int i)
+{
+ CheckRange(i, 0, 3, "(Mat3::[i]) index out of range");
+ return(row[i]);
+}
+
+inline const Vec3 &Mat3::operator [] (Int i) const
+{
+ CheckRange(i, 0, 3, "(Mat3::[i]) index out of range");
+ return(row[i]);
+}
+
+inline Real *Mat3::Ref() const
+{
+ return((Real *) row);
+}
+
+inline Mat3::Mat3(ZeroOrOne k)
+{
+ MakeDiag(k);
+}
+
+inline Mat3::Mat3(Block k)
+{
+ MakeBlock((ZeroOrOne) k);
+}
+
+inline Mat3 &Mat3::operator = (ZeroOrOne k)
+{
+ MakeDiag(k);
+
+ return(SELF);
+}
+
+inline Mat3 &Mat3::operator = (Block k)
+{
+ MakeBlock((ZeroOrOne) k);
+
+ return(SELF);
+}
+
+inline Mat3 operator * (const Real s, const Mat3 &m)
+{
+ return(m * s);
+}
+
+inline Vec3 operator * (const Mat3 &m, const Vec3 &v)
+{
+ Vec3 result;
+
+ result[0] = v[0] * m[0][0] + v[1] * m[0][1] + v[2] * m[0][2];
+ result[1] = v[0] * m[1][0] + v[1] * m[1][1] + v[2] * m[1][2];
+ result[2] = v[0] * m[2][0] + v[1] * m[2][1] + v[2] * m[2][2];
+
+ return(result);
+}
+
+inline Vec3 operator * (const Vec3 &v, const Mat3 &m)
+{
+ Vec3 result;
+
+ result[0] = v[0] * m[0][0] + v[1] * m[1][0] + v[2] * m[2][0];
+ result[1] = v[0] * m[0][1] + v[1] * m[1][1] + v[2] * m[2][1];
+ result[2] = v[0] * m[0][2] + v[1] * m[1][2] + v[2] * m[2][2];
+
+ return(result);
+}
+
+inline Vec3 &operator *= (Vec3 &v, const Mat3 &m)
+{
+ Real t0, t1;
+
+ t0 = v[0] * m[0][0] + v[1] * m[1][0] + v[2] * m[2][0];
+ t1 = v[0] * m[0][1] + v[1] * m[1][1] + v[2] * m[2][1];
+ v[2] = v[0] * m[0][2] + v[1] * m[1][2] + v[2] * m[2][2];
+ v[0] = t0;
+ v[1] = t1;
+
+ return(v);
+}
+
+#ifdef VL_ROW_ORIENT
+inline Vec2 xform(const Mat3 &m, const Vec2 &v)
+{ return(proj(Vec3(v, 1.0) * m)); }
+inline Vec3 xform(const Mat3 &m, const Vec3 &v)
+{ return(v * m); }
+inline Mat3 xform(const Mat3 &m, const Mat3 &n)
+{ return(n * m); }
+#else
+inline Vec2 xform(const Mat3 &m, const Vec2 &v)
+{ return(proj(m * Vec3(v, 1.0))); }
+inline Vec3 xform(const Mat3 &m, const Vec3 &v)
+{ return(m * v); }
+inline Mat3 xform(const Mat3 &m, const Mat3 &n)
+{ return(m * n); }
+#endif
+
+#endif
diff --git a/examples/cpp/lib/svl/Mat4.h b/examples/cpp/lib/svl/Mat4.h
new file mode 100644
index 0000000000000000000000000000000000000000..d4d1dd6d7b7f1ef5ca3f2fbadf96ab54c8c96995
--- /dev/null
+++ b/examples/cpp/lib/svl/Mat4.h
@@ -0,0 +1,189 @@
+/*
+ File: Mat4.h
+
+ Function: Defines a 4 x 4 matrix.
+
+ Author(s): Andrew Willmott
+
+ Copyright: (c) 1995-2001, Andrew Willmott
+ */
+
+#ifndef __Mat4__
+#define __Mat4__
+
+#include "svl/Vec3.h"
+#include "svl/Vec4.h"
+
+
+// --- Mat4 Class -------------------------------------------------------------
+
+class Mat4
+{
+public:
+
+ // Constructors
+
+ Mat4();
+ Mat4(Real a, Real b, Real c, Real d,
+ Real e, Real f, Real g, Real h,
+ Real i, Real j, Real k, Real l,
+ Real m, Real n, Real o, Real p);
+ Mat4(const Mat4 &m);
+ Mat4(ZeroOrOne k);
+ Mat4(Block k);
+
+ // Accessor functions
+
+ Int Rows() const { return(4); };
+ Int Cols() const { return(4); };
+
+ Vec4 &operator [] (Int i);
+ const Vec4 &operator [] (Int i) const;
+
+ Real *Ref() const;
+
+ // Assignment operators
+
+ Mat4 &operator = (const Mat4 &m);
+ Mat4 &operator = (ZeroOrOne k);
+ Mat4 &operator = (Block k);
+ Mat4 &operator += (const Mat4 &m);
+ Mat4 &operator -= (const Mat4 &m);
+ Mat4 &operator *= (const Mat4 &m);
+ Mat4 &operator *= (Real s);
+ Mat4 &operator /= (Real s);
+
+ // Comparison operators
+
+ Bool operator == (const Mat4 &m) const; // M == N?
+ Bool operator != (const Mat4 &m) const; // M != N?
+
+ // Arithmetic operators
+
+ Mat4 operator + (const Mat4 &m) const; // M + N
+ Mat4 operator - (const Mat4 &m) const; // M - N
+ Mat4 operator - () const; // -M
+ Mat4 operator * (const Mat4 &m) const; // M * N
+ Mat4 operator * (Real s) const; // M * s
+ Mat4 operator / (Real s) const; // M / s
+
+ // Initialisers
+
+ Void MakeZero(); // Zero matrix
+ Void MakeDiag(Real k = vl_one); // I
+ Void MakeBlock(Real k = vl_one); // all elts = k
+
+ // Homogeneous Transforms
+
+ Mat4& MakeHRot(const Vec3 &axis, Real theta);
+ // Rotate by theta radians about axis
+ Mat4& MakeHRot(const Vec4 &q); // Rotate by quaternion
+ Mat4& MakeHScale(const Vec3 &s); // Scale by components of s
+
+ Mat4& MakeHTrans(const Vec3 &t); // Translation by t
+
+ Mat4& Transpose(); // transpose in place
+ Mat4& AddShift(const Vec3 &t); // Concatenate shift
+
+ // Private...
+
+protected:
+
+ Vec4 row[4];
+};
+
+
+// --- Matrix operators -------------------------------------------------------
+
+Vec4 operator * (const Mat4 &m, const Vec4 &v); // m * v
+Vec4 operator * (const Vec4 &v, const Mat4 &m); // v * m
+Vec4 &operator *= (Vec4 &a, const Mat4 &m); // v *= m
+inline Mat4 operator * (Real s, const Mat4 &m); // s * m
+
+Mat4 trans(const Mat4 &m); // Transpose
+Real trace(const Mat4 &m); // Trace
+Mat4 adj(const Mat4 &m); // Adjoint
+Real det(const Mat4 &m); // Determinant
+Mat4 inv(const Mat4 &m); // Inverse
+Mat4 oprod(const Vec4 &a, const Vec4 &b);
+ // Outer product
+
+// The xform functions help avoid dependence on whether row or column
+// vectors are used to represent points and vectors.
+inline Vec4 xform(const Mat4 &m, const Vec4 &v); // Transform of v by m
+inline Vec3 xform(const Mat4 &m, const Vec3 &v); // Hom. xform of v by m
+inline Mat4 xform(const Mat4 &m, const Mat4 &n); // Xform v -> m(n(v))
+
+std::ostream &operator << (std::ostream &s, const Mat4 &m);
+std::istream &operator >> (std::istream &s, Mat4 &m);
+
+
+// --- Inlines ----------------------------------------------------------------
+
+inline Mat4::Mat4()
+{
+}
+
+inline Vec4 &Mat4::operator [] (Int i)
+{
+ CheckRange(i, 0, 4, "(Mat4::[i]) index out of range");
+ return(row[i]);
+}
+
+inline const Vec4 &Mat4::operator [] (Int i) const
+{
+ CheckRange(i, 0, 4, "(Mat4::[i]) index out of range");
+ return(row[i]);
+}
+
+inline Real *Mat4::Ref() const
+{
+ return((Real *) row);
+}
+
+inline Mat4::Mat4(ZeroOrOne k)
+{
+ MakeDiag(k);
+}
+
+inline Mat4::Mat4(Block k)
+{
+ MakeBlock((ZeroOrOne) k);
+}
+
+inline Mat4 &Mat4::operator = (ZeroOrOne k)
+{
+ MakeDiag(k);
+
+ return(SELF);
+}
+
+inline Mat4 &Mat4::operator = (Block k)
+{
+ MakeBlock((ZeroOrOne) k);
+
+ return(SELF);
+}
+
+inline Mat4 operator * (Real s, const Mat4 &m)
+{
+ return(m * s);
+}
+
+#ifdef VL_ROW_ORIENT
+inline Vec3 xform(const Mat4 &m, const Vec3 &v)
+{ return(proj(Vec4(v, 1.0) * m)); }
+inline Vec4 xform(const Mat4 &m, const Vec4 &v)
+{ return(v * m); }
+inline Mat4 xform(const Mat4 &m, const Mat4 &n)
+{ return(n * m); }
+#else
+inline Vec3 xform(const Mat4 &m, const Vec3 &v)
+{ return(proj(m * Vec4(v, 1.0))); }
+inline Vec4 xform(const Mat4 &m, const Vec4 &v)
+{ return(m * v); }
+inline Mat4 xform(const Mat4 &m, const Mat4 &n)
+{ return(m * n); }
+#endif
+
+#endif
diff --git a/examples/cpp/lib/svl/README b/examples/cpp/lib/svl/README
new file mode 100644
index 0000000000000000000000000000000000000000..b9435663a6065d9682a7d7da01028d97c7f79ce5
--- /dev/null
+++ b/examples/cpp/lib/svl/README
@@ -0,0 +1,94 @@
+Contents
+
+ The SVL vector & matrix package, version 1.5.
+
+ SVL provides 2-, 3- and 4-vector and matrix types, as well as
+ arbitrarily-sized vectors and matrices, and various useful functions
+ and arithmetic operators.
+
+ SVL is free for commercial and non-commercial use; see the LICENSE file
+ for redistribution conditions. (These apply only to the source code.
+ Binaries may be freely redistributed, no strings attached.)
+
+Author
+
+ Andrew Willmott, ajw+svl@cs.cmu.edu
+ Please send me any suggestions/comments/bug fixes you have.
+
+History
+
+ Briefly: SVL, and its more complex cousin, VL, have their roots in
+ an attempt replace the various (inevitably heavily-macroized)
+ vector/matrix packages used by the members of the CMU graphics group
+ with a single C++ one, so that people could share code more easily.
+ SVL makes heavy use of inlines to achieve performance equal to that
+ of the C macros and functions it replaced.
+
+To Make
+
+ Running 'make' with no arguments will list possible architectures to
+ compile for. The makefiles should work with most makes that have the
+ include directive, but try using gmake if something goes wrong.
+
+ If your architecture is not listed, copy and then modify the config file
+ in the makefiles directory that's closest to what you need.
+
+ To install, set the DEST variable appropriately and 'make install'.
+ DEST defaults to /usr/local.
+
+ E.g., to install into /usr on a linux machine:
+ make DEST=/usr linux install
+ from this directory.
+
+To Use
+
+ #include svl/SVL.h, and link with -lsvl. To use in debug mode, add the
+ switch "-DVL_DEBUG" to your compile line, and link with -lsvl.dbg. See
+ the documentation for more information.
+
+Documentation
+
+ Run your favourite web browser with file:doc/svl.html
+
+Changes
+
+ 1.4 Fixed problem with SVLgl.h -- had some wrappers for OpenGL calls that
+ don't actually exist(!)
+ Fixed bad memory problem with Mat::SetSize
+ Makefile for OSX, split sunos/solaris
+
+ 1.3.2 Fixed minor MSVC warnings
+ Added windows packaging
+ 1.3.1 Sunday 11/03/01
+ Reordered initializers, eliminated signed/unsigned comparisons
+ Fixed make clean problem
+ Irix CC is warning that 'Mat2::*' is being called before its inline
+ definition, but appears to be on crack as the inline '*' is
+ defined before the xform() that calls it.
+ Fixed afs-lib dependency for afs build
+ 1.3 Sunday 11/03/01
+ Split off from VL 1.3. All VL-dependent code and type
+ parameterization macros removed to make code more understandable.
+ Fix for the C++ "standard" no longer allowing stack temporaries to
+ be passes as non-const refs. (Fixing this for VL requires more
+ significant changes, so it's proved easier to fork SVL at this
+ point for the sake of simplicity.)
+ Renamed files, general cleanups.
+ Added ">>" operators for Vec and Mat.
+ Disallowed resize of a vector reference. (This caused more bugs
+ than it was worth.)
+ Removed all asserts checking that the return value of 'new' is
+ not zero. The C++ "standard" has been changed so that new
+ instead throws an exception in this situation.
+ Added option to use memcpy for general vector copies: VL_USE_MEMCPY.
+ It seems that nowadays, on Intel processors at least, using
+ memcpy is faster than loops for most cases. See note in Vec.cc.
+ Changed the packaging to match the current vl/gcl scheme.
+ Updated afs-setup rule to ensure SVLConfig.h winds up in the
+ architecture-dependent lib directory.
+ 1.2 Release to match VL 1.2.
+ 1.1.3 Fixed some windows compile issues under VC++.
+ 1.1.2 Added comments about the Real, Vec, etc. types to VL.h
+ 1.1.1 Fixed problem with += operators and gcc under linux.
+ 1.1 Added oprod operator, irix n32 CC support, changed makefiles & afs
+ support, eliminated most g++ warnings.
\ No newline at end of file
diff --git a/examples/cpp/lib/svl/SVL.h b/examples/cpp/lib/svl/SVL.h
new file mode 100644
index 0000000000000000000000000000000000000000..5f2f6435d3582b82eeb38ce5512b0ec270bbc228
--- /dev/null
+++ b/examples/cpp/lib/svl/SVL.h
@@ -0,0 +1,43 @@
+/*
+ File: SVL.h
+
+ Function: Master header for a simple version of the VL library.
+ The various classes are named Vec2, Mat3, Vec, etc.
+ Link with -lsvl, or define the symbol VL_DEBUG and
+ link with -lsvl.dbg for the debugging version.
+
+ Author(s): Andrew Willmott
+
+ Copyright: (c) 1995-2001, Andrew Willmott
+ */
+
+
+#ifndef __SVL__
+#define __SVL__
+
+#define SVL_VERSION "1.5"
+#define SVL_VER_NUM 10500
+
+#ifdef VL_DEBUG
+#define VL_CHECKING
+#endif
+
+#include <iostream>
+
+#include "svl/Basics.h"
+#include "svl/Constants.h"
+#include "svl/Utils.h"
+
+#include "svl/Vec2.h"
+#include "svl/Vec3.h"
+#include "svl/Vec4.h"
+#include "svl/Vec.h"
+
+#include "svl/Mat2.h"
+#include "svl/Mat3.h"
+#include "svl/Mat4.h"
+#include "svl/Mat.h"
+
+#include "svl/Transform.h"
+
+#endif
diff --git a/examples/cpp/lib/svl/Transform.h b/examples/cpp/lib/svl/Transform.h
new file mode 100644
index 0000000000000000000000000000000000000000..d40f329f12c3ef4819e8ea360eef9bd8fa776fa3
--- /dev/null
+++ b/examples/cpp/lib/svl/Transform.h
@@ -0,0 +1,41 @@
+/*
+ File: Transform.h
+
+ Function: Provides transformation constructors.
+
+ Author(s): Andrew Willmott
+
+ Copyright: (c) 1995-2001, Andrew Willmott
+*/
+
+#ifndef __SVL_TRANSFORM__
+#define __SVL_TRANSFORM__
+
+inline Mat2 Rot2(Real theta)
+ { Mat2 result; result.MakeRot(theta); return(result); }
+inline Mat2 Scale2(const Vec2 &s)
+ { Mat2 result; result.MakeScale(s); return(result); }
+
+inline Mat3 Rot3(const Vec3 &axis, Real theta)
+ { Mat3 result; result.MakeRot(axis, theta); return(result); }
+inline Mat3 Rot3(const Vec4 &q)
+ { Mat3 result; result.MakeRot(q); return(result); }
+inline Mat3 Scale3(const Vec3 &s)
+ { Mat3 result; result.MakeScale(s); return(result); }
+inline Mat3 HRot3(Real theta)
+ { Mat3 result; result.MakeHRot(theta); return(result); }
+inline Mat3 HScale3(const Vec2 &s)
+ { Mat3 result; result.MakeHScale(s); return(result); }
+inline Mat3 HTrans3(const Vec2 &t)
+ { Mat3 result; result.MakeHTrans(t); return(result); }
+
+inline Mat4 HRot4(const Vec3 &axis, Real theta)
+ { Mat4 result; result.MakeHRot(axis, theta); return(result); }
+inline Mat4 HRot4(const Vec4 &q)
+ { Mat4 result; result.MakeHRot(q); return(result); }
+inline Mat4 HScale4(const Vec3 &s)
+ { Mat4 result; result.MakeHScale(s); return(result); }
+inline Mat4 HTrans4(const Vec3 &t)
+ { Mat4 result; result.MakeHTrans(t); return(result); }
+
+#endif
diff --git a/examples/cpp/lib/svl/Utils.h b/examples/cpp/lib/svl/Utils.h
new file mode 100644
index 0000000000000000000000000000000000000000..803332817b15e60df0227ffc22b45a10c914f944
--- /dev/null
+++ b/examples/cpp/lib/svl/Utils.h
@@ -0,0 +1,85 @@
+/*
+ File: Utils.h
+
+ Function: Various math definitions for VL
+
+ Author(s): Andrew Willmott
+
+ Copyright: (c) 1995-2001, Andrew Willmott
+ */
+
+#ifndef __VL_MATH__
+#define __VL_MATH__
+
+#include <cstdlib>
+
+// --- Inlines ----------------------------------------------------------------
+
+// additions to arithmetic functions
+
+#ifdef VL_HAS_IEEEFP
+#include <ieeefp.h>
+#define vl_is_finite(X) finite(X)
+#elif defined (__GNUC__) && defined(__USE_MISC)
+#define vl_is_finite(X) finite(X)
+#else
+#define vl_is_finite(X) (1)
+#endif
+
+#ifdef VL_HAS_DRAND
+inline Double vl_rand()
+{ return(drand48()); }
+#else
+#ifndef RAND_MAX
+// we're on something totally sucky, like SunOS
+#define RAND_MAX (Double(1 << 30) * 4.0 - 1.0)
+#endif
+inline Double vl_rand()
+{ return(rand() / (RAND_MAX + 1.0)); }
+#endif
+
+#ifndef __CMATH__
+// GNU's complex.h defines its own abs(double)
+#ifdef VL_HAS_ABSF
+inline Float abs(Float x)
+{ return (fabsf(x)); }
+#endif
+inline Double abs(Double x)
+{ return (fabs(x)); }
+#endif
+#ifdef VL_HAS_ABSF
+inline Float len(Float x)
+{ return (fabsf(x)); }
+#endif
+inline Double len(Double x)
+{ return (fabs(x)); }
+
+inline Float sqrlen(Float r)
+{ return(sqr(r)); }
+inline Double sqrlen(Double r)
+{ return(sqr(r)); }
+
+inline Float mix(Float a, Float b, Float s)
+{ return((1.0f - s) * a + s * b); }
+inline Double mix(Double a, Double b, Double s)
+{ return((1.0 - s) * a + s * b); }
+
+inline Double sign(Double d)
+{
+ if (d < 0)
+ return(-1.0);
+ else
+ return(1.0);
+}
+
+// useful routines
+
+inline Bool IsPowerOfTwo(Int a)
+{ return((a & -a) == a); };
+
+inline Void SetReal(Float &a, Double b)
+{ a = Float(b); }
+inline Void SetReal(Double &a, Double b)
+{ a = b; }
+
+#endif
diff --git a/examples/cpp/lib/svl/Vec.h b/examples/cpp/lib/svl/Vec.h
new file mode 100644
index 0000000000000000000000000000000000000000..faef974e0fc389014e32d20a48d7d08f470411c2
--- /dev/null
+++ b/examples/cpp/lib/svl/Vec.h
@@ -0,0 +1,237 @@
+/*
+ File: Vec.h
+
+ Function: Defines a generic resizeable vector.
+
+ Author(s): Andrew Willmott
+
+ Copyright: (c) 1995-2001, Andrew Willmott
+ */
+
+#ifndef __Vec__
+#define __Vec__
+
+class Vec2;
+class Vec3;
+class Vec4;
+
+
+// --- Vec Class --------------------------------------------------------------
+
+class Vec
+{
+public:
+
+ // Constructors
+
+ Vec(); // Null vector
+ explicit Vec(Int n); // n-element vector
+ Vec(Int n, double elt0, ...); // e.g. Vec(3, 1.1, 2.0, 3.4)
+ Vec(Int n, Real *data); // Vector reference
+ Vec(const Vec &v); // Copy constructor
+ Vec(const Vec2 &v); // reference to a Vec2
+ Vec(const Vec3 &v); // reference to a Vec3
+ Vec(const Vec4 &v); // reference to a Vec4
+ Vec(Int n, ZeroOrOne); // Zero or all-ones vector
+ Vec(Int n, Axis a); // Unit vector
+ ~Vec(); // Destructor
+
+ // Accessor functions
+
+ Int Elts() const;
+
+ Real &operator [] (Int i);
+ Real operator [] (Int i) const;
+
+ Void SetSize(Int n); // Resize the vector
+ Real *Ref() const; // Return pointer to data
+
+ // Assignment operators
+
+ Vec &operator = (const Vec &v); // v = a etc.
+ Vec &operator = (ZeroOrOne k);
+ Vec &operator = (Axis a);
+ Vec &operator = (const Vec2 &v);
+ Vec &operator = (const Vec3 &v);
+ Vec &operator = (const Vec4 &v);
+
+ // In-Place operators
+
+ Vec &operator += (const Vec &v);
+ Vec &operator -= (const Vec &v);
+ Vec &operator *= (const Vec &v);
+ Vec &operator *= (Real s);
+ Vec &operator /= (const Vec &v);
+ Vec &operator /= (Real s);
+
+ // Vector initialisers
+
+ Vec &MakeZero();
+ Vec &MakeUnit(Int i, Real k = vl_one);
+ Vec &MakeBlock(Real k = vl_one);
+
+ Vec &Normalise(); // Normalise vector
+ Vec &Clamp(Real fuzz);
+ Vec &Clamp();
+
+ Bool IsRef() const { return((elts & VL_REF_FLAG) != 0); };
+
+ // Private...
+
+protected:
+
+ Real *data;
+ UInt elts;
+};
+
+
+// --- Vec Comparison Operators -----------------------------------------------
+
+Bool operator == (const Vec &a, const Vec &b);
+Bool operator != (const Vec &a, const Vec &b);
+
+
+// --- Vec Arithmetic Operators -----------------------------------------------
+
+Vec operator + (const Vec &a, const Vec &b);
+Vec operator - (const Vec &a, const Vec &b);
+Vec operator - (const Vec &v);
+Vec operator * (const Vec &a, const Vec &b);
+Vec operator * (const Vec &v, Real s);
+Vec operator / (const Vec &a, const Vec &b);
+Vec operator / (const Vec &v, Real s);
+Vec operator * (Real s, const Vec &v);
+
+Real dot(const Vec &a, const Vec &b);// v . a
+inline Real len(const Vec &v); // || v ||
+inline Real sqrlen(const Vec &v); // v . v
+inline Vec norm(const Vec &v); // v / || v ||
+inline Void normalise(Vec &v); // v = norm(v)
+
+Vec clamped(const Vec &v, Real fuzz);
+Vec clamped(const Vec &v);
+
+
+// --- Vec Input & Output -----------------------------------------------------
+
+std::ostream &operator << (std::ostream &s, const Vec &v);
+std::istream &operator >> (std::istream &s, Vec &v);
+
+
+// --- Sub-vector functions ---------------------------------------------------
+
+inline Vec sub(const Vec &v, Int start, Int length);
+inline Vec first(const Vec &v, Int length);
+inline Vec last(const Vec &v, Int length);
+
+
+// --- Vec inlines ------------------------------------------------------------
+
+inline Vec::Vec() : data(0), elts(0)
+{
+}
+
+inline Vec::Vec(Int n) : elts(n)
+{
+ Assert(n > 0,"(Vec) illegal vector size");
+
+ data = new Real[n];
+}
+
+inline Vec::Vec(Int n, Real *data) : data(data), elts(n | VL_REF_FLAG)
+{
+}
+
+inline Int Vec::Elts() const
+{
+ return(elts & VL_REF_MASK);
+}
+
+inline Real &Vec::operator [] (Int i)
+{
+ CheckRange(i, 0, Elts(), "Vec::[i]");
+
+ return(data[i]);
+}
+
+inline Real Vec::operator [] (Int i) const
+{
+ CheckRange(i, 0, Elts(), "Vec::[i]");
+
+ return(data[i]);
+}
+
+inline Real *Vec::Ref() const
+{
+ return(data);
+}
+
+inline Vec &Vec::operator = (ZeroOrOne k)
+{
+ MakeBlock(k);
+
+ return(SELF);
+}
+
+inline Vec &Vec::operator = (Axis a)
+{
+ MakeUnit(a);
+
+ return(SELF);
+}
+
+inline Real len(const Vec &v)
+{
+ return(sqrt(dot(v, v)));
+}
+
+inline Real sqrlen(const Vec &v)
+{
+ return(dot(v, v));
+}
+
+inline Vec norm(const Vec &v)
+{
+ Assert(sqrlen(v) > 0.0, "normalising length-zero vector");
+ return(v / len(v));
+}
+
+inline Void normalise(Vec &v)
+{
+ v /= len(v);
+}
+
+inline Vec sub(const Vec &v, Int start, Int length)
+{
+ Assert(start >= 0 && length > 0 && start + length <= v.Elts(),
+ "(sub(Vec)) illegal subset of vector");
+
+ return(Vec(length, v.Ref() + start));
+}
+
+inline Vec first(const Vec &v, Int length)
+{
+ Assert(length > 0 && length <= v.Elts(),
+ "(first(Vec)) illegal subset of vector");
+
+ return(Vec(length, v.Ref()));
+}
+
+inline Vec last(const Vec &v, Int length)
+{
+ Assert(length > 0 && length <= v.Elts(),
+ "(last(Vec)) illegal subset of vector");
+
+ return(Vec(length, v.Ref() + v.Elts() - length));
+}
+
+inline Vec &Vec::Normalise()
+{
+ Assert(sqrlen(SELF) > 0.0, "normalising length-zero vector");
+ SELF /= len(SELF);
+ return(SELF);
+}
+
+
+#endif
+
diff --git a/examples/cpp/lib/svl/Vec2.h b/examples/cpp/lib/svl/Vec2.h
new file mode 100644
index 0000000000000000000000000000000000000000..10eac08a3f150cab612e031526e12e6df8205a6d
--- /dev/null
+++ b/examples/cpp/lib/svl/Vec2.h
@@ -0,0 +1,364 @@
+/*
+ File: Vec2.h
+
+ Function: Defines a length-2 vector.
+
+ Author(s): Andrew Willmott
+
+ Copyright: (c) 1995-2001, Andrew Willmott
+ */
+
+#ifndef __Vec2__
+#define __Vec2__
+
+
+// --- Vec2 Class -------------------------------------------------------------
+
+
+class Vec2
+{
+public:
+
+ // Constructors
+
+ Vec2();
+ Vec2(Real x, Real y); // (x, y)
+ Vec2(const Vec2 &v); // Copy constructor
+ Vec2(ZeroOrOne k); // v[i] = vl_zero
+ Vec2(Axis k); // v[k] = 1
+
+ // Accessor functions
+
+ Real &operator [] (Int i);
+ const Real &operator [] (Int i) const;
+
+ Int Elts() const { return(2); };
+ Real *Ref() const; // Return ptr to data
+
+ // Assignment operators
+
+ Vec2 &operator = (const Vec2 &a);
+ Vec2 &operator = (ZeroOrOne k);
+ Vec2 &operator = (Axis k);
+
+ Vec2 &operator += (const Vec2 &a);
+ Vec2 &operator -= (const Vec2 &a);
+ Vec2 &operator *= (const Vec2 &a);
+ Vec2 &operator *= (Real s);
+ Vec2 &operator /= (const Vec2 &a);
+ Vec2 &operator /= (Real s);
+
+ // Comparison operators
+
+ Bool operator == (const Vec2 &a) const; // v == a?
+ Bool operator != (const Vec2 &a) const; // v != a?
+
+ // Arithmetic operators
+
+ Vec2 operator + (const Vec2 &a) const; // v + a
+ Vec2 operator - (const Vec2 &a) const; // v - a
+ Vec2 operator - () const; // -v
+ Vec2 operator * (const Vec2 &a) const; // v * a (vx * ax, ...)
+ Vec2 operator * (Real s) const; // v * s
+ Vec2 operator / (const Vec2 &a) const; // v / a (vx / ax, ...)
+ Vec2 operator / (Real s) const; // v / s
+
+ // Initialisers
+
+ Vec2 &MakeZero(); // Zero vector
+ Vec2 &MakeUnit(Int i, Real k = vl_one); // I[i]
+ Vec2 &MakeBlock(Real k = vl_one); // All-k vector
+
+ Vec2 &Normalise(); // normalise vector
+
+ // Private...
+
+protected:
+
+ Real elt[2];
+};
+
+
+// --- Vec operators ----------------------------------------------------------
+
+inline Vec2 operator * (Real s, const Vec2 &v); // s * v
+inline Real dot(const Vec2 &a, const Vec2 &b); // v . a
+inline Real len(const Vec2 &v); // || v ||
+inline Real sqrlen(const Vec2 &v); // v . v
+inline Vec2 norm(const Vec2 &v); // v / || v ||
+inline Void normalise(Vec2 &v); // v = norm(v)
+inline Vec2 cross(const Vec2 &v); // cross prod.
+
+std::ostream &operator << (std::ostream &s, const Vec2 &v);
+std::istream &operator >> (std::istream &s, Vec2 &v);
+
+
+// --- Inlines ----------------------------------------------------------------
+
+inline Real &Vec2::operator [] (Int i)
+{
+ CheckRange(i, 0, 2, "(Vec2::[i]) index out of range");
+ return(elt[i]);
+}
+
+inline const Real &Vec2::operator [] (Int i) const
+{
+ CheckRange(i, 0, 2, "(Vec2::[i]) index out of range");
+ return(elt[i]);
+}
+
+inline Vec2::Vec2()
+{
+}
+
+inline Vec2::Vec2(Real x, Real y)
+{
+ elt[0] = x;
+ elt[1] = y;
+}
+
+inline Vec2::Vec2(const Vec2 &v)
+{
+ elt[0] = v[0];
+ elt[1] = v[1];
+}
+
+inline Real *Vec2::Ref() const
+{
+ return((Real *) elt);
+}
+
+inline Vec2 &Vec2::operator = (const Vec2 &v)
+{
+ elt[0] = v[0];
+ elt[1] = v[1];
+
+ return(SELF);
+}
+
+inline Vec2 &Vec2::operator += (const Vec2 &v)
+{
+ elt[0] += v[0];
+ elt[1] += v[1];
+
+ return(SELF);
+}
+
+inline Vec2 &Vec2::operator -= (const Vec2 &v)
+{
+ elt[0] -= v[0];
+ elt[1] -= v[1];
+
+ return(SELF);
+}
+
+inline Vec2 &Vec2::operator *= (const Vec2 &v)
+{
+ elt[0] *= v[0];
+ elt[1] *= v[1];
+
+ return(SELF);
+}
+
+inline Vec2 &Vec2::operator *= (Real s)
+{
+ elt[0] *= s;
+ elt[1] *= s;
+
+ return(SELF);
+}
+
+inline Vec2 &Vec2::operator /= (const Vec2 &v)
+{
+ elt[0] /= v[0];
+ elt[1] /= v[1];
+
+ return(SELF);
+}
+
+inline Vec2 &Vec2::operator /= (Real s)
+{
+ elt[0] /= s;
+ elt[1] /= s;
+
+ return(SELF);
+}
+
+inline Vec2 Vec2::operator + (const Vec2 &a) const
+{
+ Vec2 result;
+
+ result[0] = elt[0] + a[0];
+ result[1] = elt[1] + a[1];
+
+ return(result);
+}
+
+inline Vec2 Vec2::operator - (const Vec2 &a) const
+{
+ Vec2 result;
+
+ result[0] = elt[0] - a[0];
+ result[1] = elt[1] - a[1];
+
+ return(result);
+}
+
+inline Vec2 Vec2::operator - () const
+{
+ Vec2 result;
+
+ result[0] = -elt[0];
+ result[1] = -elt[1];
+
+ return(result);
+}
+
+inline Vec2 Vec2::operator * (const Vec2 &a) const
+{
+ Vec2 result;
+
+ result[0] = elt[0] * a[0];
+ result[1] = elt[1] * a[1];
+
+ return(result);
+}
+
+inline Vec2 Vec2::operator * (Real s) const
+{
+ Vec2 result;
+
+ result[0] = elt[0] * s;
+ result[1] = elt[1] * s;
+
+ return(result);
+}
+
+inline Vec2 operator * (Real s, const Vec2 &v)
+{
+ return(v * s);
+}
+
+inline Vec2 Vec2::operator / (const Vec2 &a) const
+{
+ Vec2 result;
+
+ result[0] = elt[0] / a[0];
+ result[1] = elt[1] / a[1];
+
+ return(result);
+}
+
+inline Vec2 Vec2::operator / (Real s) const
+{
+ Vec2 result;
+
+ result[0] = elt[0] / s;
+ result[1] = elt[1] / s;
+
+ return(result);
+}
+
+inline Real dot(const Vec2 &a, const Vec2 &b)
+{
+ return(a[0] * b[0] + a[1] * b[1]);
+}
+
+inline Vec2 cross(const Vec2 &a)
+{
+ Vec2 result;
+
+ result[0] = a[1];
+ result[1] = -a[0];
+
+ return(result);
+}
+
+inline Real len(const Vec2 &v)
+{
+ return(sqrt(dot(v, v)));
+}
+
+inline Real sqrlen(const Vec2 &v)
+{
+ return(dot(v, v));
+}
+
+inline Vec2 norm(const Vec2 &v)
+{
+ Assert(sqrlen(v) > 0.0, "normalising length-zero vector");
+ return(v / len(v));
+}
+
+inline Void normalise(Vec2 &v)
+{
+ v /= len(v);
+}
+
+inline Vec2 &Vec2::MakeUnit(Int i, Real k)
+{
+ if (i == 0)
+ { elt[0] = k; elt[1] = vl_zero; }
+ else if (i == 1)
+ { elt[0] = vl_zero; elt[1] = k; }
+ else
+ _Error("(Vec2::Unit) illegal unit vector");
+ return(SELF);
+}
+
+inline Vec2 &Vec2::MakeZero()
+{
+ elt[0] = vl_zero; elt[1] = vl_zero;
+ return(SELF);
+}
+
+inline Vec2 &Vec2::MakeBlock(Real k)
+{
+ elt[0] = k; elt[1] = k;
+ return(SELF);
+}
+
+inline Vec2 &Vec2::Normalise()
+{
+ Assert(sqrlen(SELF) > 0.0, "normalising length-zero vector");
+ SELF /= len(SELF);
+ return(SELF);
+}
+
+
+inline Vec2::Vec2(ZeroOrOne k)
+{
+ elt[0] = k;
+ elt[1] = k;
+}
+
+inline Vec2::Vec2(Axis k)
+{
+ MakeUnit(k, vl_one);
+}
+
+inline Vec2 &Vec2::operator = (ZeroOrOne k)
+{
+ elt[0] = k; elt[1] = k;
+
+ return(SELF);
+}
+
+inline Vec2 &Vec2::operator = (Axis k)
+{
+ MakeUnit(k, vl_1);
+
+ return(SELF);
+}
+
+inline Bool Vec2::operator == (const Vec2 &a) const
+{
+ return(elt[0] == a[0] && elt[1] == a[1]);
+}
+
+inline Bool Vec2::operator != (const Vec2 &a) const
+{
+ return(elt[0] != a[0] || elt[1] != a[1]);
+}
+
+#endif
diff --git a/examples/cpp/lib/svl/Vec3.h b/examples/cpp/lib/svl/Vec3.h
new file mode 100644
index 0000000000000000000000000000000000000000..7a9b64f0477c272b952c2ee21752292bb52f2778
--- /dev/null
+++ b/examples/cpp/lib/svl/Vec3.h
@@ -0,0 +1,414 @@
+/*
+ File: Vec3.h
+
+ Function: Defines a length-3 vector.
+
+ Author(s): Andrew Willmott
+
+ Copyright: (c) 1995-2001, Andrew Willmott
+ */
+
+#ifndef __Vec3__
+#define __Vec3__
+
+#include "svl/Vec2.h"
+
+
+// --- Vec3 Class -------------------------------------------------------------
+
+class Vec3
+{
+public:
+
+ // Constructors
+
+ Vec3();
+ Vec3(Real x, Real y, Real z); // [x, y, z]
+ Vec3(const Vec3 &v); // Copy constructor
+ Vec3(const Vec2 &v, Real w); // Hom. 2D vector
+ Vec3(ZeroOrOne k);
+ Vec3(Axis a);
+
+ // Accessor functions
+
+ Int Elts() const { return(3); };
+
+ Real &operator [] (Int i);
+ const Real &operator [] (Int i) const;
+
+ Real *Ref() const; // Return pointer to data
+
+ // Assignment operators
+
+ Vec3 &operator = (const Vec3 &a);
+ Vec3 &operator = (ZeroOrOne k);
+ Vec3 &operator += (const Vec3 &a);
+ Vec3 &operator -= (const Vec3 &a);
+ Vec3 &operator *= (const Vec3 &a);
+ Vec3 &operator *= (Real s);
+ Vec3 &operator /= (const Vec3 &a);
+ Vec3 &operator /= (Real s);
+
+ // Comparison operators
+
+ Bool operator == (const Vec3 &a) const; // v == a?
+ Bool operator != (const Vec3 &a) const; // v != a?
+ Bool operator < (const Vec3 &a) const; // v < a?
+ Bool operator >= (const Vec3 &a) const; // v >= a?
+
+ // Arithmetic operators
+
+ Vec3 operator + (const Vec3 &a) const; // v + a
+ Vec3 operator - (const Vec3 &a) const; // v - a
+ Vec3 operator - () const; // -v
+ Vec3 operator * (const Vec3 &a) const; // v * a (vx * ax, ...)
+ Vec3 operator * (Real s) const; // v * s
+ Vec3 operator / (const Vec3 &a) const; // v / a (vx / ax, ...)
+ Vec3 operator / (Real s) const; // v / s
+
+ // Initialisers
+
+ Vec3 &MakeZero(); // Zero vector
+ Vec3 &MakeUnit(Int i, Real k = vl_one); // I[i]
+ Vec3 &MakeBlock(Real k = vl_one); // All-k vector
+
+ Vec3 &Normalise(); // normalise vector
+
+ // Private...
+
+protected:
+
+ Real elt[3];
+};
+
+
+// --- Vec operators ----------------------------------------------------------
+
+inline Vec3 operator * (Real s, const Vec3 &v); // s * v
+inline Real dot(const Vec3 &a, const Vec3 &b); // v . a
+inline Real len(const Vec3 &v); // || v ||
+inline Real sqrlen(const Vec3 &v); // v . v
+inline Vec3 norm(const Vec3 &v); // v / || v ||
+inline Void normalise(Vec3 &v); // v = norm(v)
+inline Vec3 cross(const Vec3 &a, const Vec3 &b);// a x b
+inline Vec2 proj(const Vec3 &v); // hom. projection
+
+std::ostream &operator << (std::ostream &s, const Vec3 &v);
+std::istream &operator >> (std::istream &s, Vec3 &v);
+
+
+// --- Inlines ----------------------------------------------------------------
+
+inline Real &Vec3::operator [] (Int i)
+{
+ CheckRange(i, 0, 3, "(Vec3::[i]) index out of range");
+ return(elt[i]);
+}
+
+inline const Real &Vec3::operator [] (Int i) const
+{
+ CheckRange(i, 0, 3, "(Vec3::[i]) index out of range");
+ return(elt[i]);
+}
+
+inline Vec3::Vec3()
+{
+}
+
+inline Vec3::Vec3(Real x, Real y, Real z)
+{
+ elt[0] = x;
+ elt[1] = y;
+ elt[2] = z;
+}
+
+inline Vec3::Vec3(const Vec3 &v)
+{
+ elt[0] = v[0];
+ elt[1] = v[1];
+ elt[2] = v[2];
+}
+
+inline Vec3::Vec3(const Vec2 &v, Real w)
+{
+ elt[0] = v[0];
+ elt[1] = v[1];
+ elt[2] = w;
+}
+
+inline Real *Vec3::Ref() const
+{
+ return((Real *) elt);
+}
+
+inline Vec3 &Vec3::operator = (const Vec3 &v)
+{
+ elt[0] = v[0];
+ elt[1] = v[1];
+ elt[2] = v[2];
+
+ return(SELF);
+}
+
+inline Vec3 &Vec3::operator += (const Vec3 &v)
+{
+ elt[0] += v[0];
+ elt[1] += v[1];
+ elt[2] += v[2];
+
+ return(SELF);
+}
+
+inline Vec3 &Vec3::operator -= (const Vec3 &v)
+{
+ elt[0] -= v[0];
+ elt[1] -= v[1];
+ elt[2] -= v[2];
+
+ return(SELF);
+}
+
+inline Vec3 &Vec3::operator *= (const Vec3 &a)
+{
+ elt[0] *= a[0];
+ elt[1] *= a[1];
+ elt[2] *= a[2];
+
+ return(SELF);
+}
+
+inline Vec3 &Vec3::operator *= (Real s)
+{
+ elt[0] *= s;
+ elt[1] *= s;
+ elt[2] *= s;
+
+ return(SELF);
+}
+
+inline Vec3 &Vec3::operator /= (const Vec3 &a)
+{
+ elt[0] /= a[0];
+ elt[1] /= a[1];
+ elt[2] /= a[2];
+
+ return(SELF);
+}
+
+inline Vec3 &Vec3::operator /= (Real s)
+{
+ elt[0] /= s;
+ elt[1] /= s;
+ elt[2] /= s;
+
+ return(SELF);
+}
+
+inline Vec3 Vec3::operator + (const Vec3 &a) const
+{
+ Vec3 result;
+
+ result[0] = elt[0] + a[0];
+ result[1] = elt[1] + a[1];
+ result[2] = elt[2] + a[2];
+
+ return(result);
+}
+
+inline Vec3 Vec3::operator - (const Vec3 &a) const
+{
+ Vec3 result;
+
+ result[0] = elt[0] - a[0];
+ result[1] = elt[1] - a[1];
+ result[2] = elt[2] - a[2];
+
+ return(result);
+}
+
+inline Vec3 Vec3::operator - () const
+{
+ Vec3 result;
+
+ result[0] = -elt[0];
+ result[1] = -elt[1];
+ result[2] = -elt[2];
+
+ return(result);
+}
+
+inline Vec3 Vec3::operator * (const Vec3 &a) const
+{
+ Vec3 result;
+
+ result[0] = elt[0] * a[0];
+ result[1] = elt[1] * a[1];
+ result[2] = elt[2] * a[2];
+
+ return(result);
+}
+
+inline Vec3 Vec3::operator * (Real s) const
+{
+ Vec3 result;
+
+ result[0] = elt[0] * s;
+ result[1] = elt[1] * s;
+ result[2] = elt[2] * s;
+
+ return(result);
+}
+
+inline Vec3 Vec3::operator / (const Vec3 &a) const
+{
+ Vec3 result;
+
+ result[0] = elt[0] / a[0];
+ result[1] = elt[1] / a[1];
+ result[2] = elt[2] / a[2];
+
+ return(result);
+}
+
+inline Vec3 Vec3::operator / (Real s) const
+{
+ Vec3 result;
+
+ result[0] = elt[0] / s;
+ result[1] = elt[1] / s;
+ result[2] = elt[2] / s;
+
+ return(result);
+}
+
+inline Vec3 operator * (Real s, const Vec3 &v)
+{
+ return(v * s);
+}
+
+inline Vec3 &Vec3::MakeUnit(Int n, Real k)
+{
+ if (n == 0)
+ { elt[0] = k; elt[1] = vl_zero; elt[2] = vl_zero; }
+ else if (n == 1)
+ { elt[0] = vl_zero; elt[1] = k; elt[2] = vl_zero; }
+ else if (n == 2)
+ { elt[0] = vl_zero; elt[1] = vl_zero; elt[2] = k; }
+ else
+ _Error("(Vec3::Unit) illegal unit vector");
+
+ return(SELF);
+}
+
+inline Vec3 &Vec3::MakeZero()
+{
+ elt[0] = vl_zero; elt[1] = vl_zero; elt[2] = vl_zero;
+
+ return(SELF);
+}
+
+inline Vec3 &Vec3::MakeBlock(Real k)
+{
+ elt[0] = k; elt[1] = k; elt[2] = k;
+
+ return(SELF);
+}
+
+inline Vec3 &Vec3::Normalise()
+{
+ Assert(sqrlen(SELF) > 0.0, "normalising length-zero vector");
+ SELF /= len(SELF);
+
+ return(SELF);
+}
+
+
+inline Vec3::Vec3(ZeroOrOne k)
+{
+ elt[0] = k; elt[1] = k; elt[2] = k;
+}
+
+inline Vec3 &Vec3::operator = (ZeroOrOne k)
+{
+ elt[0] = k; elt[1] = k; elt[2] = k;
+
+ return(SELF);
+}
+
+inline Vec3::Vec3(Axis a)
+{
+ MakeUnit(a, vl_one);
+}
+
+
+inline Bool Vec3::operator == (const Vec3 &a) const
+{
+ return(elt[0] == a[0] && elt[1] == a[1] && elt[2] == a[2]);
+}
+
+inline Bool Vec3::operator != (const Vec3 &a) const
+{
+ return(elt[0] != a[0] || elt[1] != a[1] || elt[2] != a[2]);
+}
+
+inline Bool Vec3::operator < (const Vec3 &a) const
+{
+ return(elt[0] < a[0] && elt[1] < a[1] && elt[2] < a[2]);
+}
+
+inline Bool Vec3::operator >= (const Vec3 &a) const
+{
+ return(elt[0] >= a[0] && elt[1] >= a[1] && elt[2] >= a[2]);
+}
+
+
+inline Real dot(const Vec3 &a, const Vec3 &b)
+{
+ return(a[0] * b[0] + a[1] * b[1] + a[2] * b[2]);
+}
+
+inline Real len(const Vec3 &v)
+{
+ return(sqrt(dot(v, v)));
+}
+
+inline Real sqrlen(const Vec3 &v)
+{
+ return(dot(v, v));
+}
+
+inline Vec3 norm(const Vec3 &v)
+{
+ Assert(sqrlen(v) > 0.0, "normalising length-zero vector");
+ return(v / len(v));
+}
+
+inline Void normalise(Vec3 &v)
+{
+ v /= len(v);
+}
+
+inline Vec3 cross(const Vec3 &a, const Vec3 &b)
+{
+ Vec3 result;
+
+ result[0] = a[1] * b[2] - a[2] * b[1];
+ result[1] = a[2] * b[0] - a[0] * b[2];
+ result[2] = a[0] * b[1] - a[1] * b[0];
+
+ return(result);
+}
+
+inline Vec2 proj(const Vec3 &v)
+{
+ Vec2 result;
+
+ Assert(v[2] != 0, "(Vec3/proj) last elt. is zero");
+
+ result[0] = v[0] / v[2];
+ result[1] = v[1] / v[2];
+
+ return(result);
+}
+
+#endif
diff --git a/examples/cpp/lib/svl/Vec4.h b/examples/cpp/lib/svl/Vec4.h
new file mode 100644
index 0000000000000000000000000000000000000000..7586d675e543aeae2e6c663063d2efc3d08350f5
--- /dev/null
+++ b/examples/cpp/lib/svl/Vec4.h
@@ -0,0 +1,379 @@
+/*
+ File: Vec4.h
+
+ Function: Defines a length-4 vector.
+
+ Author(s): Andrew Willmott
+
+ Copyright: (c) 1995-2001, Andrew Willmott
+ */
+
+#ifndef __Vec4__
+#define __Vec4__
+
+#include "svl/Vec3.h"
+
+
+// --- Vec4 Class -------------------------------------------------------------
+
+class Vec4
+{
+public:
+
+ // Constructors
+
+ Vec4();
+ Vec4(Real x, Real y, Real z, Real w); // [x, y, z, w]
+ Vec4(const Vec4 &v); // Copy constructor
+ Vec4(const Vec3 &v, Real w); // Hom. 3D vector
+ Vec4(ZeroOrOne k);
+ Vec4(Axis k);
+
+ // Accessor functions
+
+ Int Elts() const { return(4); };
+
+ Real &operator [] (Int i);
+ const Real &operator [] (Int i) const;
+
+ Real *Ref() const; // Return pointer to data
+
+ // Assignment operators
+
+ Vec4 &operator = (const Vec4 &a);
+ Vec4 &operator = (ZeroOrOne k);
+ Vec4 &operator = (Axis k);
+ Vec4 &operator += (const Vec4 &a);
+ Vec4 &operator -= (const Vec4 &a);
+ Vec4 &operator *= (const Vec4 &a);
+ Vec4 &operator *= (Real s);
+ Vec4 &operator /= (const Vec4 &a);
+ Vec4 &operator /= (Real s);
+
+ // Comparison operators
+
+ Bool operator == (const Vec4 &a) const; // v == a ?
+ Bool operator != (const Vec4 &a) const; // v != a ?
+
+ // Arithmetic operators
+
+ Vec4 operator + (const Vec4 &a) const; // v + a
+ Vec4 operator - (const Vec4 &a) const; // v - a
+ Vec4 operator - () const; // -v
+ Vec4 operator * (const Vec4 &a) const; // v * a (vx * ax, ...)
+ Vec4 operator * (Real s) const; // v * s
+ Vec4 operator / (const Vec4 &a) const; // v / a (vx / ax, ...)
+ Vec4 operator / (Real s) const; // v / s
+
+
+ // Initialisers
+
+ Vec4 &MakeZero(); // Zero vector
+ Vec4 &MakeUnit(Int i, Real k = vl_one); // kI[i]
+ Vec4 &MakeBlock(Real k = vl_one); // All-k vector
+
+ Vec4 &Normalise(); // normalise vector
+
+ // Private...
+
+protected:
+
+ Real elt[4];
+};
+
+
+// --- Vec operators ----------------------------------------------------------
+
+inline Vec4 operator * (Real s, const Vec4 &v); // Left mult. by s
+inline Real dot(const Vec4 &a, const Vec4 &b); // v . a
+inline Real len(const Vec4 &v); // || v ||
+inline Real sqrlen(const Vec4 &v); // v . v
+inline Vec4 norm(const Vec4 &v); // v / || v ||
+inline Void normalise(Vec4 &v); // v = norm(v)
+Vec4 cross(const Vec4 &a, const Vec4 &b, const Vec4 &c);
+ // 4D cross prod.
+Vec3 proj(const Vec4 &v); // hom. projection
+
+std::ostream &operator << (std::ostream &s, const Vec4 &v);
+std::istream &operator >> (std::istream &s, Vec4 &v);
+
+
+// --- Inlines ----------------------------------------------------------------
+
+inline Real &Vec4::operator [] (Int i)
+{
+ CheckRange(i, 0, 4, "(Vec4::[i]) index out of range");
+ return(elt[i]);
+}
+
+inline const Real &Vec4::operator [] (Int i) const
+{
+ CheckRange(i, 0, 4, "(Vec4::[i]) index out of range");
+ return(elt[i]);
+}
+
+
+inline Vec4::Vec4()
+{
+}
+
+inline Vec4::Vec4(Real x, Real y, Real z, Real w)
+{
+ elt[0] = x;
+ elt[1] = y;
+ elt[2] = z;
+ elt[3] = w;
+}
+
+inline Vec4::Vec4(const Vec4 &v)
+{
+ elt[0] = v[0];
+ elt[1] = v[1];
+ elt[2] = v[2];
+ elt[3] = v[3];
+}
+
+inline Vec4::Vec4(const Vec3 &v, Real w)
+{
+ elt[0] = v[0];
+ elt[1] = v[1];
+ elt[2] = v[2];
+ elt[3] = w;
+}
+
+inline Real *Vec4::Ref() const
+{
+ return((Real *) elt);
+}
+
+inline Vec4 &Vec4::operator = (const Vec4 &v)
+{
+ elt[0] = v[0];
+ elt[1] = v[1];
+ elt[2] = v[2];
+ elt[3] = v[3];
+
+ return(SELF);
+}
+
+inline Vec4 &Vec4::operator += (const Vec4 &v)
+{
+ elt[0] += v[0];
+ elt[1] += v[1];
+ elt[2] += v[2];
+ elt[3] += v[3];
+
+ return(SELF);
+}
+
+inline Vec4 &Vec4::operator -= (const Vec4 &v)
+{
+ elt[0] -= v[0];
+ elt[1] -= v[1];
+ elt[2] -= v[2];
+ elt[3] -= v[3];
+
+ return(SELF);
+}
+
+inline Vec4 &Vec4::operator *= (const Vec4 &v)
+{
+ elt[0] *= v[0];
+ elt[1] *= v[1];
+ elt[2] *= v[2];
+ elt[3] *= v[3];
+
+ return(SELF);
+}
+
+inline Vec4 &Vec4::operator *= (Real s)
+{
+ elt[0] *= s;
+ elt[1] *= s;
+ elt[2] *= s;
+ elt[3] *= s;
+
+ return(SELF);
+}
+
+inline Vec4 &Vec4::operator /= (const Vec4 &v)
+{
+ elt[0] /= v[0];
+ elt[1] /= v[1];
+ elt[2] /= v[2];
+ elt[3] /= v[3];
+
+ return(SELF);
+}
+
+inline Vec4 &Vec4::operator /= (Real s)
+{
+ elt[0] /= s;
+ elt[1] /= s;
+ elt[2] /= s;
+ elt[3] /= s;
+
+ return(SELF);
+}
+
+
+inline Vec4 Vec4::operator + (const Vec4 &a) const
+{
+ Vec4 result;
+
+ result[0] = elt[0] + a[0];
+ result[1] = elt[1] + a[1];
+ result[2] = elt[2] + a[2];
+ result[3] = elt[3] + a[3];
+
+ return(result);
+}
+
+inline Vec4 Vec4::operator - (const Vec4 &a) const
+{
+ Vec4 result;
+
+ result[0] = elt[0] - a[0];
+ result[1] = elt[1] - a[1];
+ result[2] = elt[2] - a[2];
+ result[3] = elt[3] - a[3];
+
+ return(result);
+}
+
+inline Vec4 Vec4::operator - () const
+{
+ Vec4 result;
+
+ result[0] = -elt[0];
+ result[1] = -elt[1];
+ result[2] = -elt[2];
+ result[3] = -elt[3];
+
+ return(result);
+}
+
+inline Vec4 Vec4::operator * (const Vec4 &a) const
+{
+ Vec4 result;
+
+ result[0] = elt[0] * a[0];
+ result[1] = elt[1] * a[1];
+ result[2] = elt[2] * a[2];
+ result[3] = elt[3] * a[3];
+
+ return(result);
+}
+
+inline Vec4 Vec4::operator * (Real s) const
+{
+ Vec4 result;
+
+ result[0] = elt[0] * s;
+ result[1] = elt[1] * s;
+ result[2] = elt[2] * s;
+ result[3] = elt[3] * s;
+
+ return(result);
+}
+
+inline Vec4 Vec4::operator / (const Vec4 &a) const
+{
+ Vec4 result;
+
+ result[0] = elt[0] / a[0];
+ result[1] = elt[1] / a[1];
+ result[2] = elt[2] / a[2];
+ result[3] = elt[3] / a[3];
+
+ return(result);
+}
+
+inline Vec4 Vec4::operator / (Real s) const
+{
+ Vec4 result;
+
+ result[0] = elt[0] / s;
+ result[1] = elt[1] / s;
+ result[2] = elt[2] / s;
+ result[3] = elt[3] / s;
+
+ return(result);
+}
+
+inline Vec4 operator * (Real s, const Vec4 &v)
+{
+ return(v * s);
+}
+
+inline Vec4 &Vec4::MakeZero()
+{
+ elt[0] = vl_zero; elt[1] = vl_zero; elt[2] = vl_zero; elt[3] = vl_zero;
+ return(SELF);
+}
+
+inline Vec4 &Vec4::MakeBlock(Real k)
+{
+ elt[0] = k; elt[1] = k; elt[2] = k; elt[3] = k;
+ return(SELF);
+}
+
+inline Vec4 &Vec4::Normalise()
+{
+ Assert(sqrlen(SELF) > 0.0, "normalising length-zero vector");
+ SELF /= len(SELF);
+ return(SELF);
+}
+
+inline Vec4::Vec4(ZeroOrOne k)
+{
+ MakeBlock(k);
+}
+
+inline Vec4::Vec4(Axis k)
+{
+ MakeUnit(k, vl_1);
+}
+
+inline Vec4 &Vec4::operator = (ZeroOrOne k)
+{
+ MakeBlock(k);
+
+ return(SELF);
+}
+
+inline Vec4 &Vec4::operator = (Axis k)
+{
+ MakeUnit(k, vl_1);
+
+ return(SELF);
+}
+
+
+inline Real dot(const Vec4 &a, const Vec4 &b)
+{
+ return(a[0] * b[0] + a[1] * b[1] + a[2] * b[2] + a[3] * b[3]);
+}
+
+inline Real len(const Vec4 &v)
+{
+ return(sqrt(dot(v, v)));
+}
+
+inline Real sqrlen(const Vec4 &v)
+{
+ return(dot(v, v));
+}
+
+inline Vec4 norm(const Vec4 &v)
+{
+ Assert(sqrlen(v) > 0.0, "normalising length-zero vector");
+ return(v / len(v));
+}
+
+inline Void normalise(Vec4 &v)
+{
+ v /= len(v);
+}
+
+#endif