This section of the archives stores flipcode's complete Developer Toolbox collection, featuring a variety of mini-articles and source code contributions from our readers.

 

  Smallest Enclosing Spheres
  Submitted by



Bounding spheres have a reputation of being a bad bounding volume because they are said to have a bad fit. They only seem interesting for the very rough first culling test or collision test because of their cheap intersection tests. They do not have a bad fit because of their geometric qualities, but rather because we use the wrong algorithms for constructing them. All methods used in 3D game engines are bad approximation of what could be the smallest sphere. Most algorithms that compute the smallest enclosing sphere use a refinement algorithm that is hard to implement and slow so we don't bother using them... In this COTD I will try to introduce you to an algorithm that computes the smallest enclosing sphere in an analytical way which is easy to implement and has very good performance. The algorithm was developed by Prof. Dr. Emo Welzl in 1991. You can download his paper here: Smallest enclosing disks (balls and ellipsoids). The first implementation was done by Dr. Bernd Gärtner, and as far as I know it's hasn't been implemented by many other people. He has also written a paper about his work: Fast and Robust Smallest Enclosing Balls, and you can download his code in the CGAL project here: Computational Geometry Algorithms Library. So why don't we use his implementation? Well, he has written a very general version of the algorithm, and in my humble opinion, it's a little too general to use it in a 3D game engine. It allows any dimension so it is over-complicated and in no way specialized for 3D. I can't give a thorough explanation of the algorithm, but I'll give you the basics because I think it's important you know a bit about an algorithm before you use it. I strongly suggest to read both papers, but be aware that they are relatively difficult. First of all, you need to understand that a sphere can be defined by four points on it's boundary. If we give less than four points, there's an infinite number of solutions, but only one of them is the smallest sphere and they can all be computed easily. Since this means a sphere can only be defined by at most four points, we only need to find those points of the object we want to enclose. This will be our starting point, so we first need functions that give the smallest sphere for 0 to 4 points. I'll only give one example to show how this works. Suppose we get two points and we are supposed to construct the smallest sphere through them. It should be clear that the midpoint of the line segment connecting the two points is the sphere's center, and half the length of the segment is it's radius. Sphere's through three and four points are a bit harder to compute but the code at the end of this COTD gives your the necessary formula's. I've implemented all these functions as constructors of the Sphere class, with 0 to 4 points as arguments. Now we have these functions we need to find the points that will lie on the smallest sphere's boundary, construct the sphere trough them and we're done. However, this is not a trivial task, and you need to know one more thing before I show you Welzl's algorithm. We can construct the smallest enclosing sphere by adding one point after another, so eventually all points will be added and we have the smallest sphere enclosing our object. If we already have the smallest sphere for a certain set of points, and we add another point, two things can happen: it can either be inside our outside the current sphere. If it's inside, we can immediately go to the next point because we already have the smallest enclosing sphere. If it's not inside, we have to construct a new sphere. It can be proven that the current point will be one the new sphere's boundary, and this point will also be important for the next steps. All this suggests we can use a recursive algorithm...

Sphere Sphere::MiniBall(Point P[], int p, Point B[], int b)
{
	if(p == 0 || b == 4)
	{
		switch(b)
		{
		case 0:
			return Sphere();
		case 1:
			return Sphere(B[0]);
		case 2:
			return Sphere(B[0], B[1]);
		case 3:
			return Sphere(B[0], B[1], B[2]);
		case 4:
			return Sphere(B[0], B[1], B[2], B[3]);
		}
	}

Sphere MB = MiniBall(P, p - 1, B, b);

if(MB.d(P[p - 1]) 0) // Signed distance to sphere { for(int i = p - 2; i = 0; i--) { Point T = P[i]; P[i] = P[i + 1]; P[i + 1] = T; }

MB = MiniBall(P + 1, p - 1, B, b + 1); }

return MB; }

This is the most basic form of Welzl's algorithm. It takes two arguments: a set of points P to be enclosed and a set of points B that must lie on the sphere's boundary. The sets don't have common points. If we want to construct the smallest enclosing sphere of an object, we pass all vertices in P, and pass an empty set for B because we don't know these points yet. To understand this recursive algorithm, we only have to understand one step. If you try to follow the recursive calls you'll definitely loose track and won't understand it. So let's start with the first if() statement. Here we check for the trivial situations. If P is empty or B has the maximum amount of points that can lie on the sphere's boundary, we have to compute the smallest enclosing sphere of B. This is done with the Sphere constructors I've talked about. After this we can return because we have finished the task of constructing the smallest sphere around P with B on it's boundary. If we don't have the trivial case, we construct the smallest sphere without the last point of P. Don't worry about how this recursive call works, just assume it will return the correct result. Now that we have this sphere we can check if the last point of P is inside of it or not. If it's inside we return because the task has been finished, if not, we know the last point of P will be a point of B. So we move the point to B and again we do a recursive call. There's a little trick here. Because we know the last point of P is an important point, we can move it to the front of the list so 'other' recursive calls will only have to add the last point of P, and Welzl proves this also speeds up the algorithm. Welzl first calls the function and then moves the point, but for this implementation it was easier to do it the other way around because I use one array for P and B. This implementation is not very efficient, but at least it compiles and works. We can make a lot of improvements that will make it faster and more reliable. First of all, we don't need to copy the points because we don't change anything, so we can use pointers. Secondly, moving the last point to the front is done here by swapping each element, which can be inefficient with huge point sets. Gärtner uses a linked list to move points in O(1) time. If we have large point sets, the amount of recursive calls can make the stack overflow quite easily, so we need to restructure the algorithm to reduce recursive calls. Gärtner does this by using a loop that adds the points from the first to the last instead of just adding the last and doing a new recursive call. A last optimization that also makes the algorithm more robust is to add the 'farthest' points first, so a lot of points will already be inside the sphere and it 'grows' very fast. However, for small point sets like used in 3D engines this is slower because it takes some time to find the farthest point. Check Gärtner's paper for details about his 'pivoting' method. So what do we win with all this? Well, if you're using bounding sphere's for visibility culling against the viewing frustum, you need the best 'plane intersection chance'. What I mean by this is that if the bounding volume is intersected by a plane, there should be a big chance the enclosed object is intersected too. The worst case is when a thin rod has to be enclosed. Since we now have an algorithm that will compute the smallest sphere possible, we know it will go through the rod's endpoints. A little bit of math shows the intersection chance is 2 / Pi = 0.64! In my opinion this is extremely good for the worst case. A cube has an intersection chance of sqrt(3 / 4) = 0.87. Combine this with the cheap intersection test of a sphere with a plane, and you'll cull a lot of geometry almost for free. Finally, here's the implementation I'm currently using: MiniBall.zip (100 kB). I have made some optimizations, but not all. For instance, the move-to-front is still done by swapping each pointer, this is because it's not a bottleneck of the algorithm. I also provide a DOS test program. It constructs 10,000 points and compares a classical algorithm with Welzl's algorithm. The approximation algorithm still constructs pretty small spheres, but that's only because I use random points in the unit cube. With the average 3D model the difference between the two is very noticeable. Best regards, Nicolas "Nick" Capens


Currently browsing [MiniBall.zip] (54,486 bytes) - [BoundingVolume.h] - (419 bytes)

#ifndef BOUNDINGVOLUME_H
#define BOUNDINGVOLUME_H

namespace MiniBall { class ViewFrustum; class Plane;

class BoundingVolume { public: virtual bool visible(const ViewFrustum &viewFrustum) const; virtual bool intersect(const Plane &plane) const; virtual bool frontSide(const Plane &plane) const = 0; virtual bool backSide(const Plane &plane) const; }; }

#endif // BOUNDINGVOLUME_H

Currently browsing [MiniBall.zip] (54,486 bytes) - [Matrix.h] - (7,140 bytes)

#ifndef MATRIX_H
#define MATRIX_H

namespace MiniBall { class Vector; class Point; class Quaternion;

class Matrix { private: // Row major order float m[4][4];

public: Matrix(); Matrix(const int i); Matrix(const float m[16]); Matrix(const float m[4][4]); Matrix(float m11, float m12, float m13, float m21, float m22, float m23, float m31, float m32, float m33); Matrix(float m11, float m12, float m13, float m14, float m21, float m22, float m23, float m24, float m31, float m32, float m33, float m34, float m41, float m42, float m43, float m44); Matrix(const Vector &v1, const Vector &v2, const Vector &v3); // Column vectors Matrix &operator=(const Matrix &N);

static Matrix diag(float m11, float m22, float m33, float m44);

operator float*();

Matrix operator+() const; Matrix operator-() const;

Matrix operator!() const; // Inverse Matrix operator~() const; // Transpose Matrix &operator+=(const Matrix &N); Matrix &operator-=(const Matrix &N); Matrix &operator*=(float s); Matrix &operator*=(const Matrix &N); Matrix &operator/=(float s);

float *operator[](int i); // Access element [row, col], starting with [0][0] const float *operator[](int i) const;

float &operator()(int i, int j); // Access element (row, col), starting with (1, 1) const float &operator()(int i, int j) const;

friend bool operator==(const Matrix &M, const Matrix &N); friend bool operator!=(const Matrix &M, const Matrix &N);

friend Matrix operator+(const Matrix &M, const Matrix &N); friend Matrix operator-(const Matrix &M, const Matrix &N); friend Matrix operator*(float s, const Matrix &M); friend Matrix operator*(const Matrix &M, const Matrix &N); friend Matrix operator/(const Matrix &M, float s); friend Matrix operator^(const Matrix &M, const int n);

Vector operator*(const Vector& v) const; Point operator*(const Point& P) const; Quaternion operator*(const Quaternion& Q) const;

friend Vector operator*(const Vector &v, const Matrix &M); friend Point operator*(const Point &P, const Matrix &M); friend Quaternion operator*(const Quaternion &Q, const Matrix &M);

friend Vector &operator*=(Vector &v, const Matrix &M); friend Point &operator*=(Point &P, const Matrix &M); friend Quaternion &operator*=(Quaternion &Q, const Matrix &M);

static float det(const Matrix &M); static float det(float m11); static float det(float m11, float m12, float m21, float m22); static float det(float m11, float m12, float m13, float m21, float m22, float m23, float m31, float m32, float m33); static float det(float m11, float m12, float m13, float m14, float m21, float m22, float m23, float m24, float m31, float m32, float m33, float m34, float m41, float m42, float m43, float m44); static float det(const Vector &v1, const Vector &v2, const Vector &v3);

static float tr(const Matrix &M);

Matrix &orthogonalise(); // Gram-Schmidt orthogonalisation of 3x3 submatrix static Matrix eulerRotate(float x, float y, float z); static Matrix eulerRotate(const Vector& v);

static Matrix translate(float x, float y, float z); static Matrix translate(const Vector& v); }; }

#include "Vector.h"

#include <assert.h>

namespace MiniBall { inline Matrix::Matrix() { }

inline Matrix::Matrix(const int i) { const float s = (float)i;

Matrix &M = *this;

M(1, 1) = s; M(1, 2) = 0; M(1, 3) = 0; M(1, 4) = 0; M(2, 1) = 0; M(2, 2) = s; M(2, 3) = 0; M(2, 4) = 0; M(3, 1) = 0; M(3, 2) = 0; M(3, 3) = s; M(3, 4) = 0; M(4, 1) = 0; M(4, 2) = 0; M(4, 3) = 0; M(4, 4) = s; }

inline Matrix::Matrix(const float m[16]) { Matrix &M = *this;

M(1, 1) = m[0]; M(1, 2) = m[1]; M(1, 3) = m[2]; M(1, 4) = m[3]; M(2, 1) = m[4]; M(2, 2) = m[5]; M(2, 3) = m[6]; M(2, 4) = m[7]; M(3, 1) = m[8]; M(3, 2) = m[8]; M(3, 3) = m[10]; M(3, 4) = m[11]; M(4, 1) = m[12]; M(4, 2) = m[13]; M(4, 3) = m[14]; M(4, 4) = m[15]; }

inline Matrix::Matrix(const float m[4][4]) { Matrix &M = *this;

M[0][0] = m[0][0]; M[0][1] = m[0][1]; M[0][2] = m[0][2]; M[0][3] = m[0][3]; M[1][0] = m[1][0]; M[1][1] = m[1][1]; M[1][2] = m[1][2]; M[1][3] = m[1][3]; M[2][0] = m[2][0]; M[2][1] = m[2][1]; M[2][2] = m[2][2]; M[2][3] = m[2][3]; M[3][0] = m[3][0]; M[3][1] = m[3][1]; M[3][2] = m[3][2]; M[3][3] = m[3][3]; }

inline Matrix::Matrix(float m11, float m12, float m13, float m21, float m22, float m23, float m31, float m32, float m33) { Matrix &M = *this;

M(1, 1) = m11; M(1, 2) = m12; M(1, 3) = m13; M(1, 4) = 0; M(2, 1) = m21; M(2, 2) = m22; M(2, 3) = m23; M(2, 4) = 0; M(3, 1) = m31; M(3, 2) = m32; M(3, 3) = m33; M(3, 4) = 0; M(4, 1) = 0; M(4, 2) = 0; M(4, 3) = 0; M(4, 4) = 1; }

inline Matrix::Matrix(float m11, float m12, float m13, float m14, float m21, float m22, float m23, float m24, float m31, float m32, float m33, float m34, float m41, float m42, float m43, float m44) { Matrix &M = *this;

M(1, 1) = m11; M(1, 2) = m12; M(1, 3) = m13; M(1, 4) = m14; M(2, 1) = m21; M(2, 2) = m22; M(2, 3) = m23; M(2, 4) = m24; M(3, 1) = m31; M(3, 2) = m32; M(3, 3) = m33; M(3, 4) = m34; M(4, 1) = m41; M(4, 2) = m42; M(4, 3) = m43; M(4, 4) = m44; }

inline Matrix::Matrix(const Vector &v1, const Vector &v2, const Vector &v3) { Matrix &M = *this;

M(1, 1) = v1.x; M(1, 2) = v2.x; M(1, 3) = v3.x; M(1, 4) = 0; M(2, 1) = v1.y; M(2, 2) = v2.y; M(2, 3) = v3.y; M(2, 4) = 0; M(3, 1) = v1.z; M(3, 2) = v2.z; M(3, 3) = v3.z; M(3, 4) = 0; M(4, 1) = 0; M(4, 2) = 0; M(4, 3) = 0; M(4, 4) = 1; }

inline Matrix &Matrix::operator=(const Matrix &N) { Matrix &M = *this;

M(1, 1) = N(1, 1); M(1, 2) = N(1, 2); M(1, 3) = N(1, 3); M(1, 4) = N(1, 4); M(2, 1) = N(2, 1); M(2, 2) = N(2, 2); M(2, 3) = N(2, 3); M(2, 4) = N(2, 4); M(3, 1) = N(3, 1); M(3, 2) = N(3, 2); M(3, 3) = N(3, 3); M(3, 4) = N(3, 4); M(4, 1) = N(4, 1); M(4, 2) = N(4, 2); M(4, 3) = N(4, 3); M(4, 4) = N(4, 4);

return M; }

inline float *Matrix::operator[](int i) { return m[i]; }

inline const float *Matrix::operator[](int i) const { return m[i]; }

inline float &Matrix::operator()(int i, int j) { return m[i - 1][j - 1]; }

inline const float &Matrix::operator()(int i, int j) const { return m[i - 1][j - 1]; }

inline Matrix operator^(const Matrix &M, const int n) { switch(n) { case -1: return !M; case 0: return 1; case 1: return M; case 2: return M*M; case 3: return M*M*M; default: assert(n >= -1 && n <= 3); // FIXME: Compile time assert? return 0; } } }

#endif // MATRIX_H

Currently browsing [MiniBall.zip] (54,486 bytes) - [Plane.h] - (2,219 bytes)

#ifndef PLANE_H
#define PLANE_H

#include "Vector.h"

namespace MiniBall { class Matrix; class Point;

class Plane { public: union { struct { float A; float B; float C; }; struct { Vector n; }; };

float D; // Distance to origin along normal Plane(); Plane(const Plane &p); Plane(const Vector &n, float D); // Normal and distance to origin Plane(const Vector &n, const Point &P); // Normal and point on plane Plane(const Point &P0, const Point &P1, const Point &P2); // Through three points Plane(float A, float B, float C, float D); // Plane equation Plane &operator=(const Plane &p);

Plane &operator+(); Plane operator-() const; // Flip normal Plane &operator*=(const Matrix &A); // Transform plane by matrix (post-multiply) friend Plane operator*(const Plane &p, const Matrix &A); // Transform plane by matrix (post-multiply) friend Plane operator*(const Matrix &A, const Plane &p); // Transform plane by matrix (pre-multiply) friend float operator^(const Plane &p1, const Plane &p2); // Angle between planes float d(const Point &P) const; // Oriented distance between point and plane static float d(const Point &P, const Plane &p); // Oriented distance between point and plane static float d(const Plane &p, const Point &P); // Oriented distance between plane and point Plane &normalise(); // Normalise the Plane equation }; }

#include "Point.h"

namespace MiniBall { inline Plane::Plane() { }

inline Plane::Plane(const Plane &p) { n = p.n; D = p.D; }

inline Plane::Plane(const Vector &p_n, float p_D) { n = p_n; D = p_D; }

inline Plane::Plane(const Vector &p_n, const Point &P) { n = p_n; D = -n * P; }

inline Plane::Plane(const Point &P0, const Point &P1, const Point &P2) { n = (P1 - P0) % (P2 - P0); D = -n * P0; }

inline Plane::Plane(float p_A, float p_B, float p_C, float p_D) { A = p_A; B = p_B; C = p_C; D = p_D; }

inline Plane &Plane::operator=(const Plane &p) { n = p.n; D = p.D;

return *this; } }

#endif // PLANE_H

Currently browsing [MiniBall.zip] (54,486 bytes) - [Point.h] - (2,218 bytes)

#ifndef POINT_H
#define POINT_H

#include "Tuple.h"

namespace MiniBall { class Vector;

class Point : public Tuple<3> { public: Point(); Point(const int i); Point(const Point &P); Point(const Vector &v); Point(float Px, float Py, float Pz);

Point &operator=(const Point &P);

float getX() const; void setX(float x); __declspec(property(get = getX, put = setX)) float x;

float getY() const; void setY(float y); __declspec(property(get = getY, put = setY)) float y;

float getZ() const; void setZ(float z); __declspec(property(get = getZ, put = setZ)) float z;

Point &operator+=(const Vector &v); Point &operator-=(const Vector &v);

friend Point operator+(const Point &P, const Vector &v); friend Point operator-(const Point &P, const Vector &v);

friend Vector operator-(const Point &P, const Point &Q);

float d(const Point &P) const; // Distance between two points float d2(const Point &P) const; // Squared distance between two points static float d(const Point &P, const Point &Q); // Distance between two points static float d2(const Point &P, const Point &Q); // Squared distance between two points }; }

#include "Vector.h"

namespace MiniBall { inline Point::Point() { }

inline Point::Point(const int i) { const float s = (float)i;

x = s; y = s; z = s; }

inline Point::Point(const Point &P) { x = P.x; y = P.y; z = P.z; }

inline Point::Point(const Vector &v) { x = v.x; y = v.y; z = v.z; }

inline Point::Point(float P_x, float P_y, float P_z) { x = P_x; y = P_y; z = P_z; }

inline Point &Point::operator=(const Point &P) { x = P.x; y = P.y; z = P.z;

return *this; }

inline float Point::getX() const { return element[0]; }

inline void Point::setX(float x_) { element[0] = x_; }

inline float Point::getY() const { return element[1]; }

inline void Point::setY(float y_) { element[1] = y_; }

inline float Point::getZ() const { return element[2]; }

inline void Point::setZ(float z_) { element[2] = z_; } }

#endif // POINT_H

Currently browsing [MiniBall.zip] (54,486 bytes) - [Quaternion.h] - (4,193 bytes)

#ifndef QUATERNION_H
#define QUATERNION_H

#include "Tuple.h"

namespace MiniBall { class Point; class Vector; class Matrix;

class Quaternion : public Tuple<4> { public: Quaternion(); Quaternion(float Q_w); Quaternion(const Quaternion &Q); Quaternion(const Quaternion &Q, const Quaternion &R, float t); // Slerp Quaternion(const Vector &Q_v, float Q_w = 0.0f); Quaternion(const Point &P, float Q_w = 1.0f); Quaternion(const Matrix &A); // Must be a rotation Matrix Quaternion(float Q_x, float Q_y, float Q_z, float Q_w);

Quaternion &operator=(const Quaternion &Q);

float getX() const; void setX(float x); __declspec(property(get = getX, put = setX)) float x;

float getY() const; void setY(float y); __declspec(property(get = getY, put = setY)) float y;

float getZ() const; void setZ(float z); __declspec(property(get = getZ, put = setZ)) float z;

float getW() const; void setW(float w); __declspec(property(get = getW, put = setW)) float w;

operator float() const; operator Vector() const; operator Matrix() const;

Quaternion operator+() const; Quaternion operator-() const;

Quaternion operator!() const; // Inverse Quaternion operator~() const; // Conjugate Quaternion &operator+=(const Quaternion &V); Quaternion &operator-=(const Quaternion &V); Quaternion &operator*=(float s); Quaternion &operator*=(const Quaternion &Q); Quaternion &operator/=(float s);

friend bool operator==(const Quaternion &Q, const Quaternion &R); friend bool operator!=(const Quaternion &Q, const Quaternion &R);

friend Quaternion operator+(const Quaternion &Q, const Quaternion &R); friend Quaternion operator-(const Quaternion &Q, const Quaternion &R); friend Quaternion operator*(const Quaternion &Q, const Quaternion &R); friend Quaternion operator*(float s, const Quaternion &Q); friend Quaternion operator*(const Quaternion &Q, float s); friend Quaternion operator/(const Quaternion &Q, float s);

friend float N(const Quaternion &Q); // Norm friend float N2(const Quaternion &Q); // Squared norm friend Quaternion slerp(const Quaternion &Q, const Quaternion &R, float t); }; }

#include "Matrix.h" #include "Point.h" #include "Vector.h"

#include <math.h>

namespace MiniBall { inline Quaternion::Quaternion() { }

inline Quaternion::Quaternion(const Quaternion &Q) { x = Q.x; y = Q.y; z = Q.z; w = Q.w; }

inline Quaternion::Quaternion(const Quaternion &Q, const Quaternion &R, float t) { float a = acosf(Q.x * R.x + Q.y * R.y + Q.z * R.z + Q.w * R.w); // Angle float s = sinf(a);

*this = (Q * sinf(a * (1 - t)) + R * sinf(a * t)) / s; }

inline Quaternion::Quaternion(float Q_w) { x = 0; y = 0; z = 0; w = Q_w; }

inline Quaternion::Quaternion(const Vector &Q_v, float Q_w) { x = Q_v.x; y = Q_v.y; z = Q_v.z; w = Q_w; }

inline Quaternion::Quaternion(const Point &P, float Q_w) { x = P.x; y = P.y; z = P.z; w = Q_w; }

inline Quaternion::Quaternion(const Matrix &A) { w = 0.5F * sqrtf(Matrix::tr(A)); x = (A(3, 2) - A(2, 3)) / (4 * w); y = (A(1, 3) - A(3, 1)) / (4 * w); z = (A(2, 1) - A(1, 2)) / (4 * w); }

inline Quaternion::Quaternion(float Q_x, float Q_y, float Q_z, float Q_w) { x = Q_x; y = Q_y; z = Q_z; w = Q_w; }

inline Quaternion &Quaternion::operator=(const Quaternion &Q) { x = Q.x; y = Q.y; z = Q.z; w = Q.w;

return *this; }

inline float Quaternion::getX() const { return element[0]; }

inline void Quaternion::setX(float x_) { element[0] = x_; }

inline float Quaternion::getY() const { return element[1]; }

inline void Quaternion::setY(float y_) { element[1] = y_; }

inline float Quaternion::getZ() const { return element[2]; }

inline void Quaternion::setZ(float z_) { element[2] = z_; }

inline float Quaternion::getW() const { return element[3]; }

inline void Quaternion::setW(float w_) { element[3] = w_; } }

#endif // QUATERNION_H

Currently browsing [MiniBall.zip] (54,486 bytes) - [Sphere.h] - (1,599 bytes)

#ifndef SPHERE_H
#define SPHERE_H

#include "Point.h"

#include <math.h>

namespace MiniBall { class Vector; class Matrix;

class Sphere { public: Point center; float radius;

Sphere(); Sphere(const Sphere &X); Sphere(const Point &O); // Point-Sphere Sphere(const Point &O, float R); // Center and radius (not squared) Sphere(const Point &O, const Point &A); // Sphere through two points Sphere(const Point &O, const Point &A, const Point &B); // Sphere through three points Sphere(const Point &O, const Point &A, const Point &B, const Point &C); // Sphere through four points Sphere &operator=(const Sphere &S);

float d(const Point &P) const; // Distance from p to boundary of the Sphere float d2(const Point &P) const; // Square distance from p to boundary of the Sphere static float d(const Sphere &S, const Point &P); // Distance from p to boundary of the Sphere static float d(const Point &P, const Sphere &S); // Distance from p to boundary of the Sphere static float d2(const Sphere &S, const Point &P); // Square distance from p to boundary of the Sphere static float d2(const Point &P, const Sphere &S); // Square distance from p to boundary of the Sphere static Sphere miniBall(Point P[], unsigned int p); // Smallest enclosing sphere static Sphere smallBall(Point P[], unsigned int p); // Enclosing sphere approximation private: static Sphere recurseMini(Point *P[], unsigned int p, unsigned int b = 0); };

extern const float radiusEpsilon; }

#endif // SPHERE_H

Currently browsing [MiniBall.zip] (54,486 bytes) - [Timer.h] - (216 bytes)

#ifndef TIMER_H
#define TIMER_H

namespace MiniBall { class Timer { public: Timer();

static double seconds(); static double ticks(); static void pause(double p); }; }

#endif // TIMER_H

Currently browsing [MiniBall.zip] (54,486 bytes) - [Tuple.h] - (1,375 bytes)

#ifndef TUPLE_H
#define TUPLE_H

namespace MiniBall { template<unsigned int n, class Type = float> class Tuple { protected: Type element[n];

public: Tuple(); Tuple(Type e[n]);

Tuple &operator=(const Tuple &T);

float &operator[](int i); float &operator()(int i);

const float &operator[](int i) const; const float &operator()(int i) const; };

template<unsigned int n, class Type> inline Tuple<n, Type>::Tuple() { }

template<unsigned int n, class Type> inline Tuple<n, Type>::Tuple(Type e[n]) { for(int i = 0; i < n; i++) element[i] = e[i]; }

template<unsigned int n, class Type> inline Tuple<n, Type> &Tuple<n, Type>::operator=(const Tuple<n, Type> &T) { for(int i = 0; i < n; i++) element[i] = T.element[i];

return *this; }

template<unsigned int n, class Type> inline float &Tuple<n, Type>::operator[](int i) { return element[i]; }

template<unsigned int n, class Type> inline float &Tuple<n, Type>::operator()(int i) { return element[(i - 1) & 3]; }

template<unsigned int n, class Type> inline const float &Tuple<n, Type>::operator[](int i) const { return element[i]; }

template<unsigned int n, class Type> inline const float &Tuple<n, Type>::operator()(int i) const { return element[(i - 1) & 3]; } }

#endif // TUPLE_H

Currently browsing [MiniBall.zip] (54,486 bytes) - [Vector.h] - (2,921 bytes)

#ifndef VECTOR_H
#define VECTOR_H

#include "Tuple.h"

namespace MiniBall { class Point;

class Vector : public Tuple<3> { public: Vector(); Vector(const int i); Vector(const Vector &v); Vector(const Point &p); Vector(float v_x, float v_y, float v_z);

Vector &operator=(const Vector &v);

float getX() const; void setX(float x); __declspec(property(get = getX, put = setX)) float x;

float getY() const; void setY(float y); __declspec(property(get = getY, put = setY)) float y;

float getZ() const; void setZ(float z); __declspec(property(get = getZ, put = setZ)) float z;

Vector operator+() const; Vector operator-() const;

Vector &operator+=(const Vector &v); Vector &operator-=(const Vector &v); Vector &operator*=(float s); Vector &operator/=(float s);

friend bool operator==(const Vector &u, const Vector &v); friend bool operator!=(const Vector &u, const Vector &v);

friend Vector operator+(const Vector &u, const Vector &v); friend Vector operator-(const Vector &u, const Vector &v); friend float operator*(const Vector &u, const Vector &v); // Dot product friend Vector operator*(float s, const Vector &v); friend Vector operator*(const Vector &v, float s); friend Vector operator/(const Vector &v, float s); friend float operator^(const Vector &v, const int n); // Vector power friend float operator^(const Vector &u, const Vector &v); // Angle between vectors friend Vector operator%(const Vector &u, const Vector &v); // Cross product static float N(const Vector &v); // Norm static float N2(const Vector &v); // Squared norm }; }

#include "Point.h"

#include <assert.h>

namespace MiniBall { inline Vector::Vector() { }

inline Vector::Vector(const int i) { const float s = (float)i;

x = s; y = s; z = s; }

inline Vector::Vector(const Vector &v) { x = v.x; y = v.y; z = v.z; }

inline Vector::Vector(const Point &P) { x = P.x; y = P.y; z = P.z; }

inline Vector::Vector(float v_x, float v_y, float v_z) { x = v_x; y = v_y; z = v_z; }

inline Vector &Vector::operator=(const Vector &v) { x = v.x; y = v.y; z = v.z;

return *this; }

inline float Vector::getX() const { return element[0]; }

inline void Vector::setX(float x_) { element[0] = x_; }

inline float Vector::getY() const { return element[1]; }

inline void Vector::setY(float y_) { element[1] = y_; }

inline float Vector::getZ() const { return element[2]; }

inline void Vector::setZ(float z_) { element[2] = z_; }

inline float operator^(const Vector &v, const int n) { switch(n) { case 2: return v*v; default: assert(n == 2); // FIXME: Compile time assert? return 1; } } }

#endif // VECTOR_H

Currently browsing [MiniBall.zip] (54,486 bytes) - [ViewFrustum.h] - (540 bytes)

#ifndef VIEWFRUSTUM_H
#define VIEWFRUSTUM_H

#include "Plane.h"

namespace MiniBall { class Matrix;

class ViewFrustum { public: Plane nearPlane; Plane farPlane; Plane leftPlane; Plane rightPlane; Plane topPlane; Plane bottomPlane;

ViewFrustum(); ViewFrustum(const ViewFrustum &viewFrustum);

~ViewFrustum();

ViewFrustum &operator=(const ViewFrustum &viewFrustum);

friend ViewFrustum operator*(const Matrix &matrix, const ViewFrustum &viewFrustum); }; }

#endif // VIEWFRUSTUM_H

Currently browsing [MiniBall.zip] (54,486 bytes) - [BoundingSphere.cpp] - (1,738 bytes)

#include "BoundingSphere.h"

#include "ViewFrustum.h" #include "Plane.h"

namespace MiniBall { BoundingSphere::BoundingSphere() { }

BoundingSphere::BoundingSphere(const BoundingSphere &boundingSphere) { center = boundingSphere.center; radius = boundingSphere.radius; }

BoundingSphere::BoundingSphere(const Sphere &sphere) { center = sphere.center; radius = sphere.radius; }

BoundingSphere &BoundingSphere::operator=(const BoundingSphere &boundingSphere) { center = boundingSphere.center; radius = boundingSphere.radius;

return *this; }

BoundingSphere &BoundingSphere::operator=(const Sphere &sphere) { center = sphere.center; radius = sphere.radius;

return *this; }

BoundingSphere::~BoundingSphere() { }

bool BoundingSphere::visible(const ViewFrustum &viewFrustum) const { // NOTE: Conservative test, not exact if(Plane::d(viewFrustum.nearPlane, center) < -radius) return false;

if(Plane::d(viewFrustum.farPlane, center) < -radius) return false;

if(Plane::d(viewFrustum.leftPlane, center) < -radius) return false;

if(Plane::d(viewFrustum.rightPlane, center) < -radius) return false;

if(Plane::d(viewFrustum.topPlane, center) < -radius) return false;

if(Plane::d(viewFrustum.bottomPlane, center) < -radius) return false;

return true; }

bool BoundingSphere::intersect(const Plane &plane) const { return fabs(Plane::d(plane, center)) < radius; }

bool BoundingSphere::frontSide(const Plane &plane) const { return Plane::d(plane, center) > radius; }

bool BoundingSphere::backSide(const Plane &plane) const { return Plane::d(plane, center) < -radius; } }

Currently browsing [MiniBall.zip] (54,486 bytes) - [BoundingVolume.cpp] - (715 bytes)

#include "BoundingVolume.h"

#include "ViewFrustum.h" #include "Plane.h"

namespace MiniBall { bool BoundingVolume::visible(const ViewFrustum &viewFrustum) const { // NOTE: Conservative test, not exact return !backSide(viewFrustum.nearPlane) && !backSide(viewFrustum.farPlane) && !backSide(viewFrustum.leftPlane) && !backSide(viewFrustum.rightPlane) && !backSide(viewFrustum.topPlane) && !backSide(viewFrustum.bottomPlane); }

bool BoundingVolume::intersect(const Plane &plane) const { return !frontSide(plane) && !backSide(plane); }

bool BoundingVolume::backSide(const Plane &plane) const { return frontSide(-plane); } }

Currently browsing [MiniBall.zip] (54,486 bytes) - [Main.cpp] - (2,149 bytes)

#include "BoundingSphere.h"
#include "Point.h"
#include "Timer.h"

#include <stdlib.h> #include <stdio.h> #include <time.h> #include <conio.h>

float rand(float min, float max) { return min + (max - min) * rand() / (float)RAND_MAX; }

int main() { using namespace MiniBall;

const float epsilon = radiusEpsilon + 1e-5f; // Used for testing correctness of algorithms BoundingSphere BS; unsigned int p = 0;

do { printf("Enter number of points: ");

if(!scanf("%d", &p)) { printf("Invalid input.\n\n"); break; }

printf("\nAllocating memory.\n");

Point *P = new Point[p];

if(!P) { printf("Memory could not be allocated.\n\n"); break; }

printf("\nConstructing random points: ");

srand((unsigned)time(0));

for(unsigned int i = 0; i < 10; i++) { for(unsigned int j = i * p / 10; j < (i + 1) * p / 10; j++) { P[j].x = rand(-1, 1); P[j].y = rand(-1, 1); P[j].z = rand(-1, 1); }

printf("."); }

printf("\n\nBounding sphere approximation:\n\n");

double time0 = Timer::seconds();

BS = Sphere::smallBall(P, p);

double time1 = Timer::seconds();

printf(" Construction time: %f\n", time1 - time0); printf(" Center = (%f, %f, %f)\n", BS.center.x, BS.center.y, BS.center.z); printf(" Radius = %f\n\n", BS.radius);

for(i = 0; i < p; i++) { if(Sphere::d(BS, P[i]) >= -epsilon) { printf(" Boundary point distance check: %.3f\n", Sphere::d(BS, P[i])); } }

printf("\nSmallest enclosing sphere:\n\n");

time0 = Timer::seconds();

BS = Sphere::miniBall(P, p);

time1 = Timer::seconds();

printf(" Construction time: %f\n", time1 - time0); printf(" Center = (%f, %f, %f)\n", BS.center.x, BS.center.y, BS.center.z); printf(" Radius = %f\n\n", BS.radius);

for(i = 0; i < p; i++) { if(Sphere::d(BS, P[i]) >= -epsilon) { printf(" Boundary point distance check: %.3f\n", Sphere::d(BS, P[i])); } }

printf("\n");

if(P) { delete[] P; P = 0; } } while(true);

return 0; }

Currently browsing [MiniBall.zip] (54,486 bytes) - [Matrix.cpp] - (14,586 bytes)

#include "Matrix.h"

#include "Point.h" #include "Quaternion.h"

#include <math.h>

namespace MiniBall { Matrix Matrix::diag(float m11, float m22, float m33, float m44) { return Matrix(m11, 0, 0, 0, 0, m22, 0, 0, 0, 0, m33, 0, 0, 0, 0, m44); }

Matrix::operator float*() { return &(*this)(1, 1); }

Matrix Matrix::operator+() const { return *this; }

Matrix Matrix::operator-() const { const Matrix &M = *this;

return Matrix(-M(1, 1), -M(1, 2), -M(1, 3), -M(1, 4), -M(2, 1), -M(2, 2), -M(2, 3), -M(2, 4), -M(3, 1), -M(3, 2), -M(3, 3), -M(3, 4), -M(4, 1), -M(4, 2), -M(4, 3), -M(4, 4)); }

Matrix Matrix::operator!() const { const Matrix &M = *this; Matrix I;

float M3344 = M(3, 3) * M(4, 4) - M(4, 3) * M(3, 4); float M2344 = M(2, 3) * M(4, 4) - M(4, 3) * M(2, 4); float M2334 = M(2, 3) * M(3, 4) - M(3, 3) * M(2, 4); float M3244 = M(3, 2) * M(4, 4) - M(4, 2) * M(3, 4); float M2244 = M(2, 2) * M(4, 4) - M(4, 2) * M(2, 4); float M2234 = M(2, 2) * M(3, 4) - M(3, 2) * M(2, 4); float M3243 = M(3, 2) * M(4, 3) - M(4, 2) * M(3, 3); float M2243 = M(2, 2) * M(4, 3) - M(4, 2) * M(2, 3); float M2233 = M(2, 2) * M(3, 3) - M(3, 2) * M(2, 3); float M1344 = M(1, 3) * M(4, 4) - M(4, 3) * M(1, 4); float M1334 = M(1, 3) * M(3, 4) - M(3, 3) * M(1, 4); float M1244 = M(1, 2) * M(4, 4) - M(4, 2) * M(1, 4); float M1234 = M(1, 2) * M(3, 4) - M(3, 2) * M(1, 4); float M1243 = M(1, 2) * M(4, 3) - M(4, 2) * M(1, 3); float M1233 = M(1, 2) * M(3, 3) - M(3, 2) * M(1, 3); float M1324 = M(1, 3) * M(2, 4) - M(2, 3) * M(1, 4); float M1224 = M(1, 2) * M(2, 4) - M(2, 2) * M(1, 4); float M1223 = M(1, 2) * M(2, 3) - M(2, 2) * M(1, 3);

// Adjoint Matrix I(1, 1) = M(2, 2) * M3344 - M(3, 2) * M2344 + M(4, 2) * M2334; I(2, 1) = -M(2, 1) * M3344 + M(3, 1) * M2344 - M(4, 1) * M2334; I(3, 1) = M(2, 1) * M3244 - M(3, 1) * M2244 + M(4, 1) * M2234; I(4, 1) = -M(2, 1) * M3243 + M(3, 1) * M2243 - M(4, 1) * M2233;

I(1, 2) = -M(1, 2) * M3344 + M(3, 2) * M1344 - M(4, 2) * M1334; I(2, 2) = M(1, 1) * M3344 - M(3, 1) * M1344 + M(4, 1) * M1334; I(3, 2) = -M(1, 1) * M3244 + M(3, 1) * M1244 - M(4, 1) * M1234; I(4, 2) = M(1, 1) * M3243 - M(3, 1) * M1243 + M(4, 1) * M1233;

I(1, 3) = M(1, 2) * M2344 - M(2, 2) * M1344 + M(4, 2) * M1324; I(2, 3) = -M(1, 1) * M2344 + M(2, 1) * M1344 - M(4, 1) * M1324; I(3, 3) = M(1, 1) * M2244 - M(2, 1) * M1244 + M(4, 1) * M1224; I(4, 3) = -M(1, 1) * M2243 + M(2, 1) * M1243 - M(4, 1) * M1223;

I(1, 4) = -M(1, 2) * M2334 + M(2, 2) * M1334 - M(3, 2) * M1324; I(2, 4) = M(1, 1) * M2334 - M(2, 1) * M1334 + M(3, 1) * M1324; I(3, 4) = -M(1, 1) * M2234 + M(2, 1) * M1234 - M(3, 1) * M1224; I(4, 4) = M(1, 1) * M2233 - M(2, 1) * M1233 + M(3, 1) * M1223;

// Division by determinant I /= M(1, 1) * I(1, 1) + M(2, 1) * I(1, 2) + M(3, 1) * I(1, 3) + M(4, 1) * I(1, 4);

return I; }

Matrix Matrix::operator~() const { const Matrix &M = *this;

return Matrix(M(1, 1), M(2, 1), M(3, 1), M(4, 1), M(1, 2), M(2, 2), M(3, 2), M(4, 2), M(1, 3), M(2, 3), M(3, 3), M(4, 3), M(1, 4), M(2, 4), M(3, 4), M(4, 4)); }

Matrix &Matrix::operator+=(const Matrix &N) { Matrix &M = *this;

M(1, 1) += N(1, 1); M(1, 2) += N(1, 2); M(1, 3) += N(1, 3); M(1, 4) += N(1, 4); M(2, 1) += N(2, 1); M(2, 2) += N(2, 2); M(2, 3) += N(2, 3); M(2, 4) += N(2, 4); M(3, 1) += N(3, 1); M(3, 2) += N(3, 2); M(3, 3) += N(3, 3); M(3, 4) += N(3, 4); M(4, 1) += N(4, 1); M(4, 2) += N(4, 2); M(4, 3) += N(4, 3); M(4, 4) += N(4, 4);

return M; }

Matrix &Matrix::operator-=(const Matrix &N) { Matrix &M = *this;

M(1, 1) -= N(1, 1); M(1, 2) -= N(1, 2); M(1, 3) -= N(1, 3); M(1, 4) -= N(1, 4); M(2, 1) -= N(2, 1); M(2, 2) -= N(2, 2); M(2, 3) -= N(2, 3); M(2, 4) -= N(2, 4); M(3, 1) -= N(3, 1); M(3, 2) -= N(3, 2); M(3, 3) -= N(3, 3); M(3, 4) -= N(3, 4); M(4, 1) -= N(4, 1); M(4, 2) -= N(4, 2); M(4, 3) -= N(4, 3); M(4, 4) -= N(4, 4);

return M; }

Matrix &Matrix::operator*=(float s) { Matrix &M = *this;

M(1, 1) *= s; M(1, 2) *= s; M(1, 3) *= s; M(1, 4) *= s; M(2, 1) *= s; M(2, 2) *= s; M(2, 3) *= s; M(2, 4) *= s; M(3, 1) *= s; M(3, 2) *= s; M(3, 3) *= s; M(3, 4) *= s; M(4, 1) *= s; M(4, 2) *= s; M(4, 3) *= s; M(4, 4) *= s;

return M; }

Matrix &Matrix::operator*=(const Matrix &M) { return *this = *this * M; }

Matrix &Matrix::operator/=(float s) { float r = 1.0f / s;

return *this *= r; }

bool operator==(const Matrix &M, const Matrix &N) { if(M(1, 1) == N(1, 1) && M(1, 2) == N(1, 2) && M(1, 3) == N(1, 3) && M(1, 4) == N(1, 4) && M(2, 1) == N(2, 1) && M(2, 2) == N(2, 2) && M(2, 3) == N(2, 3) && M(2, 4) == N(2, 4) && M(3, 1) == N(3, 1) && M(3, 2) == N(3, 2) && M(3, 3) == N(3, 3) && M(3, 4) == N(3, 4) && M(4, 1) == N(4, 1) && M(4, 2) == N(4, 2) && M(4, 3) == N(4, 3) && M(4, 4) == N(4, 4)) return true; else return false; }

bool operator!=(const Matrix &M, const Matrix &N) { if(M(1, 1) != N(1, 1) || M(1, 2) != N(1, 2) || M(1, 3) != N(1, 3) || M(1, 4) != N(1, 4) || M(2, 1) != N(2, 1) || M(2, 2) != N(2, 2) || M(2, 3) != N(2, 3) || M(2, 4) != N(2, 4) || M(3, 1) != N(3, 1) || M(3, 2) != N(3, 2) || M(3, 3) != N(3, 3) || M(3, 4) != N(3, 4) || M(4, 1) != N(4, 1) || M(4, 2) != N(4, 2) || M(4, 3) != N(4, 3) || M(4, 4) != N(4, 4)) return true; else return false; }

Matrix operator+(const Matrix &M, const Matrix &N) { return Matrix(M(1, 1) + N(1, 1), M(1, 2) + N(1, 2), M(1, 3) + N(1, 3), M(1, 4) + N(1, 4), M(2, 1) + N(2, 1), M(2, 2) + N(2, 2), M(2, 3) + N(2, 3), M(2, 4) + N(2, 4), M(3, 1) + N(3, 1), M(3, 2) + N(3, 2), M(3, 3) + N(3, 3), M(3, 4) + N(3, 4), M(4, 1) + N(4, 1), M(4, 2) + N(4, 2), M(4, 3) + N(4, 3), M(4, 4) + N(4, 4)); }

Matrix operator-(const Matrix &M, const Matrix &N) { return Matrix(M(1, 1) - N(1, 1), M(1, 2) - N(1, 2), M(1, 3) - N(1, 3), M(1, 4) - N(1, 4), M(2, 1) - N(2, 1), M(2, 2) - N(2, 2), M(2, 3) - N(2, 3), M(2, 4) - N(2, 4), M(3, 1) - N(3, 1), M(3, 2) - N(3, 2), M(3, 3) - N(3, 3), M(3, 4) - N(3, 4), M(4, 1) - N(4, 1), M(4, 2) - N(4, 2), M(4, 3) - N(4, 3), M(4, 4) - N(4, 4)); }

Matrix operator*(float s, const Matrix &M) { return Matrix(s * M(1, 1), s * M(1, 2), s * M(1, 3), s * M(1, 4), s * M(2, 1), s * M(2, 2), s * M(2, 3), s * M(2, 4), s * M(3, 1), s * M(3, 2), s * M(3, 3), s * M(3, 4), s * M(4, 1), s * M(4, 2), s * M(4, 3), s * M(4, 4)); }

Matrix operator*(const Matrix &M, float s) { return Matrix(M(1, 1) * s, M(1, 2) * s, M(1, 3) * s, M(1, 4) * s, M(2, 1) * s, M(2, 2) * s, M(2, 3) * s, M(2, 4) * s, M(3, 1) * s, M(3, 2) * s, M(3, 3) * s, M(3, 4) * s, M(4, 1) * s, M(4, 2) * s, M(4, 3) * s, M(4, 4) * s); }

Matrix operator*(const Matrix &M, const Matrix &N) { return Matrix(M(1, 1) * N(1, 1) + M(1, 2) * N(2, 1) + M(1, 3) * N(3, 1) + M(1, 4) * N(4, 1), M(1, 1) * N(1, 2) + M(1, 2) * N(2, 2) + M(1, 3) * N(3, 2) + M(1, 4) * N(4, 2), M(1, 1) * N(1, 3) + M(1, 2) * N(2, 3) + M(1, 3) * N(3, 3) + M(1, 4) * N(4, 3), M(1, 1) * N(1, 4) + M(1, 2) * N(2, 4) + M(1, 3) * N(3, 4) + M(1, 4) * N(4, 4), M(2, 1) * N(1, 1) + M(2, 2) * N(2, 1) + M(2, 3) * N(3, 1) + M(2, 4) * N(4, 1), M(2, 1) * N(1, 2) + M(2, 2) * N(2, 2) + M(2, 3) * N(3, 2) + M(2, 4) * N(4, 2), M(2, 1) * N(1, 3) + M(2, 2) * N(2, 3) + M(2, 3) * N(3, 3) + M(2, 4) * N(4, 3), M(2, 1) * N(1, 4) + M(2, 2) * N(2, 4) + M(2, 3) * N(3, 4) + M(2, 4) * N(4, 4), M(3, 1) * N(1, 1) + M(3, 2) * N(2, 1) + M(3, 3) * N(3, 1) + M(3, 4) * N(4, 1), M(3, 1) * N(1, 2) + M(3, 2) * N(2, 2) + M(3, 3) * N(3, 2) + M(3, 4) * N(4, 2), M(3, 1) * N(1, 3) + M(3, 2) * N(2, 3) + M(3, 3) * N(3, 3) + M(3, 4) * N(4, 3), M(3, 1) * N(1, 4) + M(3, 2) * N(2, 4) + M(3, 3) * N(3, 4) + M(3, 4) * N(4, 4), M(4, 1) * N(1, 1) + M(4, 2) * N(2, 1) + M(4, 3) * N(3, 1) + M(4, 4) * N(4, 1), M(4, 1) * N(1, 2) + M(4, 2) * N(2, 2) + M(4, 3) * N(3, 2) + M(4, 4) * N(4, 2), M(4, 1) * N(1, 3) + M(4, 2) * N(2, 3) + M(4, 3) * N(3, 3) + M(4, 4) * N(4, 3), M(4, 1) * N(1, 4) + M(4, 2) * N(2, 4) + M(4, 3) * N(3, 4) + M(4, 4) * N(4, 4)); }

Matrix operator/(const Matrix &M, float s) { float r = 1.0f / s;

return M * s; }

Vector Matrix::operator*(const Vector &v) const { const Matrix &M = *this;

return Vector(M(1, 1) * v.x + M(1, 2) * v.y + M(1, 3) * v.z, M(2, 1) * v.x + M(2, 2) * v.y + M(2, 3) * v.z, M(3, 1) * v.x + M(3, 2) * v.y + M(3, 3) * v.z); }

Point Matrix::operator*(const Point &P) const { const Matrix &M = *this;

return Point(M(1, 1) * P.x + M(1, 2) * P.y + M(1, 3) * P.z + M(1, 4), M(2, 1) * P.x + M(2, 2) * P.y + M(2, 3) * P.z + M(2, 4), M(3, 1) * P.x + M(3, 2) * P.y + M(3, 3) * P.z + M(3, 4)); }

Quaternion Matrix::operator*(const Quaternion &Q) const { const Matrix &M = *this;

return Quaternion(M(1, 1) * Q.x + M(1, 2) * Q.y + M(1, 3) * Q.z + M(1, 4) * Q.w, M(2, 1) * Q.x + M(2, 2) * Q.y + M(2, 3) * Q.z + M(2, 4) * Q.w, M(3, 1) * Q.x + M(3, 2) * Q.y + M(3, 3) * Q.z + M(3, 4) * Q.w, M(4, 1) * Q.x + M(4, 2) * Q.y + M(4, 3) * Q.z + M(4, 4) * Q.w); }

Vector operator*(const Vector &v, const Matrix &M) { return Vector(v.x * M(1, 1) + v.y * M(2, 1) + v.z * M(3, 1) + M(4, 1), v.x * M(1, 2) + v.y * M(2, 2) + v.z * M(3, 2) + M(4, 2), v.x * M(1, 3) + v.y * M(2, 3) + v.z * M(3, 3) + M(4, 3)); }

Point operator*(const Point &P, const Matrix &M) { return Point(P.x * M(1, 1) + P.y * M(2, 1) + P.z * M(3, 1), P.x * M(1, 2) + P.y * M(2, 2) + P.z * M(3, 2), P.x * M(1, 3) + P.y * M(2, 3) + P.z * M(3, 3)); }

Quaternion operator*(const Quaternion &Q, const Matrix &M) { return Quaternion(Q.x * M(1, 1) + Q.y * M(2, 1) + Q.z + M(3, 1) + Q.w * M(4, 1), Q.x * M(1, 2) + Q.y * M(2, 2) + Q.z + M(3, 2) + Q.w * M(4, 2), Q.x * M(1, 3) + Q.y * M(2, 3) + Q.z + M(3, 3) + Q.w * M(4, 3), Q.x * M(1, 4) + Q.y * M(2, 4) + Q.z + M(3, 4) + Q.w * M(4, 4)); }

Vector &operator*=(Vector &v, const Matrix &M) { return v = v * M; }

Point &operator*=(Point &P, const Matrix &M) { return P = P * M; }

Quaternion &operator*=(Quaternion &Q, const Matrix &M) { return Q = Q * M; }

float Matrix::det(const Matrix &M) { float M3344 = M(3, 3) * M(4, 4) - M(4, 3) * M(3, 4); float M2344 = M(2, 3) * M(4, 4) - M(4, 3) * M(2, 4); float M2334 = M(2, 3) * M(3, 4) - M(3, 3) * M(2, 4); float M1344 = M(1, 3) * M(4, 4) - M(4, 3) * M(1, 4); float M1334 = M(1, 3) * M(3, 4) - M(3, 3) * M(1, 4); float M1324 = M(1, 3) * M(2, 4) - M(2, 3) * M(1, 4);

return M(1, 1) * (M(2, 2) * M3344 - M(3, 2) * M2344 + M(4, 2) * M2334) - M(2, 1) * (M(1, 2) * M3344 - M(3, 2) * M1344 + M(4, 2) * M1334) + M(3, 1) * (M(1, 2) * M2344 - M(2, 2) * M1344 + M(4, 2) * M1324) - M(4, 1) * (M(1, 2) * M2334 - M(2, 2) * M1334 + M(3, 2) * M1324); }

float Matrix::det(float m11) { return m11; }

float Matrix::det(float m11, float m12, float m21, float m22) { return m11 * m22 - m12 * m21; }

float Matrix::det(float m11, float m12, float m13, float m21, float m22, float m23, float m31, float m32, float m33) { return m11 * (m22 * m33 - m32 * m23) - m21 * (m12 * m33 - m32 * m13) + m31 * (m12 * m23 - m22 * m13); }

float Matrix::det(float m11, float m12, float m13, float m14, float m21, float m22, float m23, float m24, float m31, float m32, float m33, float m34, float m41, float m42, float m43, float m44) { float M3344 = m33 * m44 - m43 * m34; float M2344 = m23 * m44 - m43 * m24; float M2334 = m23 * m34 - m33 * m24; float M1344 = m13 * m44 - m43 * m14; float M1334 = m13 * m34 - m33 * m14; float M1324 = m13 * m24 - m23 * m14;

return m11 * (m22 * M3344 - m32 * M2344 + m42 * M2334) - m21 * (m12 * M3344 - m32 * M1344 + m42 * M1334) + m31 * (m12 * M2344 - m22 * M1344 + m42 * M1324) - m41 * (m12 * M2334 - m22 * M1334 + m32 * M1324); }

float Matrix::det(const Vector &v1, const Vector &v2, const Vector &v3) { return v1 * (v2 % v3); }

float Matrix::tr(const Matrix &M) { return M(1, 1) + M(2, 2) + M(3, 3) + M(4, 4); }

Matrix &Matrix::orthogonalise() { // NOTE: Numerically instable, won't return exact the same result when already orhtogonal Matrix &M = *this;

Vector v1(M(1, 1), M(2, 1), M(3, 1)); Vector v2(M(1, 2), M(2, 2), M(3, 2)); Vector v3(M(1, 3), M(2, 3), M(3, 3));

v2 -= v1 * (v1 * v2) / (v1 * v1); v3 -= v1 * (v1 * v3) / (v1 * v1); v3 -= v2 * (v2 * v3) / (v2 * v2);

v1 /= N(v1); v2 /= N(v2); v3 /= N(v3);

M(1, 1) = v1.x; M(1, 2) = v2.x; M(1, 3) = v3.x; M(2, 1) = v1.y; M(2, 2) = v2.y; M(2, 3) = v3.y; M(3, 1) = v1.z; M(3, 2) = v2.z; M(3, 3) = v3.z;

return *this; }

Matrix Matrix::eulerRotate(float x, float y, float z) { float cz = cosf(z); float sz = sinf(z); float cx = cosf(x); float sx = sinf(x); float cy = cosf(y); float sy = sinf(y);

float sxsy = sx * sy; float sxcy = sx * cy;

return Matrix(cy * cz - sy * sx * sz, -cy * sz - sy * sx * cz, -sy * cx, cx * sz, cx * cz, -sx, sy * cz + cy * sx * sz, -sy * sz + cy * sx * cz, cy * cx); }

Matrix Matrix::eulerRotate(const Vector &v) { return eulerRotate(v.x, v.y, v.z); }

Matrix Matrix::translate(float x, float y, float z) { return Matrix(1, 0, 0, x, 0, 1, 0, y, 0, 0, 1, z, 0, 0, 0, 1); }

Matrix Matrix::translate(const Vector& V) { return translate(V.x, V.y, V.z); } }

Currently browsing [MiniBall.zip] (54,486 bytes) - [Plane.cpp] - (1,574 bytes)

#include "Plane.h"

#include "Matrix.h"

#include <math.h>

namespace MiniBall { Plane &Plane::operator+() { return *this; }

Plane Plane::operator-() const { return Plane(-n, -D); }

Plane &Plane::operator*=(const Matrix &M) { return *this = *this * M; }

Plane operator*(const Plane &p, const Matrix &T) { Matrix M = ~T^-1;

return Plane(p.A * M(1, 1) + p.B * M(2, 1) + p.C * M(3, 1) + p.D * M(4, 1), p.A * M(1, 2) + p.B * M(2, 2) + p.C * M(3, 2) + p.D * M(4, 2), p.A * M(1, 3) + p.B * M(2, 3) + p.C * M(3, 3) + p.D * M(4, 3), p.A * M(1, 4) + p.B * M(2, 4) + p.C * M(3, 4) + p.D * M(4, 4)); }

Plane operator*(const Matrix &T, const Plane &p) { Matrix M = ~T^-1;

return Plane(M(1, 1) * p.A + M(1, 2) * p.B + M(1, 3) * p.C + M(1, 4) * p.D, M(2, 1) * p.A + M(2, 2) * p.B + M(2, 3) * p.C + M(2, 4) * p.D, M(3, 1) * p.A + M(3, 2) * p.B + M(3, 3) * p.C + M(3, 4) * p.D, M(4, 1) * p.A + M(4, 2) * p.B + M(4, 3) * p.C + M(4, 4) * p.D); }

float operator^(const Plane &p1, const Plane &p2) { return acosf(p1.n * p2.n / (Vector::N(p1.n) * Vector::N(p2.n))); }

float Plane::d(const Point &P) const { return P * n + D; }

float Plane::d(const Point &P, const Plane &p) { return P * p.n + p.D; }

float Plane::d(const Plane &p, const Point &P) { return p.n * P + p.D; }

Plane &Plane::normalise() { float l = Vector::N(n);

n /= l; D /= l;

return *this; } }

Currently browsing [MiniBall.zip] (54,486 bytes) - [Point.cpp] - (942 bytes)

#include "Point.h"

namespace MiniBall { Point &Point::operator+=(const Vector &v) { x += v.x; y += v.y; z += v.z;

return *this; }

Point &Point::operator-=(const Vector &v) { x -= v.x; y -= v.y; z -= v.z;

return *this; }

Point operator+(const Point &P, const Vector &v) { return Point(P.x + v.x, P.y + v.y, P.z + v.z); }

Point operator-(const Point &P, const Vector &v) { return Point(P.x - v.x, P.y - v.y, P.z - v.z); }

Vector operator-(const Point &P, const Point &Q) { return Vector(P.x - Q.x, P.y - Q.y, P.z - Q.z); }

float Point::d(const Point &P) const { return Vector::N(*this - P); }

float Point::d2(const Point &P) const { return Vector::N2(*this - P); }

float Point::d(const Point &P, const Point &Q) { return Vector::N(P - Q); }

float Point::d2(const Point &P, const Point &Q) { return Vector::N2(P - Q); } }

Currently browsing [MiniBall.zip] (54,486 bytes) - [Quaternion.cpp] - (3,313 bytes)

#include "Quaternion.h"

namespace MiniBall { Quaternion::operator float() const { return w; }

Quaternion::operator Vector() const { return Vector(x, y, z); }

Quaternion::operator Matrix() const { return Matrix(1 - 2 * (y * y + z * z), 2 * (x * y - w * z), 2 * (x * z + w * y), 2 * (x * y + w * z), 1 - 2 * (x * x + z * z), 2 * (y * z - w * x), 2 * (x * z - w * y), 2 * (y * z + w * x), 1 - 2 * (x * x + y * y)); }

Quaternion Quaternion::operator+() const { return *this; }

Quaternion Quaternion::operator-() const { return Quaternion(-x, -y, -z, -w); }

Quaternion Quaternion::operator!() const { return ~(*this) / N2(*this); }

Quaternion Quaternion::operator~() const { return Quaternion(-x, -y, -z, w); }

Quaternion &Quaternion::operator+=(const Quaternion &Q) { x += Q.x; y += Q.y; z += Q.z; w += Q.w;

return *this; }

Quaternion &Quaternion::operator-=(const Quaternion &Q) { x -= Q.x; y -= Q.y; z -= Q.z; w -= Q.w;

return *this; }

Quaternion &Quaternion::operator*=(float s) { x *= s; y *= s; z *= s; w *= s;

return *this; }

Quaternion &Quaternion::operator*=(const Quaternion &Q) { float Qx = w * Q.x + x * Q.w + y * Q.z - z * Q.y; float Qy = w * Q.y - x * Q.z + y * Q.w + z * Q.x; float Qz = w * Q.z + x * Q.y - y * Q.x + z * Q.w; float Qw = w * Q.w - x * Q.x - y * Q.y - z * Q.z;

w = Qw; x = Qx; y = Qy; z = Qz;

return (*this); }

Quaternion &Quaternion::operator/=(float s) { float r = 1.0f / s;

return *this *= r; }

bool operator==(const Quaternion &Q, const Quaternion &R) { if(Q.x == R.x && Q.y == R.y && Q.z == R.z && Q.w == R.w) return true; else return false; }

bool operator!=(const Quaternion &Q, const Quaternion &R) { if(Q.x != R.x || Q.y != R.y || Q.z != R.z || Q.w != R.w) return true; else return false; }

Quaternion operator+(const Quaternion &Q, const Quaternion &R) { return Quaternion(Q.x + R.x, Q.y + R.y, Q.z + R.z, Q.w + R.w); }

Quaternion operator-(const Quaternion &Q, const Quaternion &R) { return Quaternion(Q.x - R.x, Q.y - R.y, Q.z - R.z, Q.w - R.w); }

Quaternion operator*(const Quaternion &Q, const Quaternion &R) { float x = +Q.x * R.w + Q.y * R.z - Q.z * R.y + Q.w * R.x; float y = -Q.x * R.z + Q.y * R.w + Q.z * R.x + Q.w * R.y; float z = +Q.x * R.y - Q.y * R.x + Q.z * R.w + Q.w * R.z; float w = -Q.x * R.x - Q.y * R.y - Q.z * R.z + Q.w * R.w;

return Quaternion(x, y, z, w); }

Quaternion operator*(float s, const Quaternion &Q) { return Quaternion(s * Q.x, s * Q.y, s * Q.z, s * Q.w); }

Quaternion operator*(const Quaternion &Q, float s) { return Quaternion(Q.x * s, Q.y * s, Q.z * s, Q.w * s); }

Quaternion operator/(const Quaternion &Q, float s) { float r = 1.0f / s;

return Q * r; }

float N(const Quaternion &Q) { return sqrtf(Q.x*Q.x + Q.y*Q.y + Q.z*Q.z + Q.w*Q.w); }

float N2(const Quaternion &Q) { return Q.x*Q.x + Q.y*Q.y + Q.z*Q.z + Q.w*Q.w; }

Quaternion slerp(const Quaternion &Q, const Quaternion &R, float t) { return Quaternion(Q, R, t); } }

Currently browsing [MiniBall.zip] (54,486 bytes) - [Sphere.cpp] - (3,612 bytes)

#include "Sphere.h"

#include "Vector.h" #include "Matrix.h"

#include "stdlib.h"

namespace MiniBall { const float radiusEpsilon = 1e-4f; // NOTE: To avoid numerical inaccuracies Sphere::Sphere() { radius = -1; }

Sphere::Sphere(const Sphere &S) { radius = S.radius; center = S.center; }

Sphere::Sphere(const Point &O) { radius = 0 + radiusEpsilon; center = O; }

Sphere::Sphere(const Point &O, float R) { radius = R; center = O; }

Sphere::Sphere(const Point &O, const Point &A) { Vector a = A - O;

Vector o = 0.5f * a;

radius = Vector::N(o) + radiusEpsilon; center = O + o; }

Sphere::Sphere(const Point &O, const Point &A, const Point &B) { Vector a = A - O; Vector b = B - O;

float Denominator = 2.0f * ((a % b) * (a % b));

Vector o = ((b^2) * ((a % b) % a) + (a^2) * (b % (a % b))) / Denominator;

radius = Vector::N(o) + radiusEpsilon; center = O + o; }

Sphere::Sphere(const Point &O, const Point &A, const Point &B, const Point &C) { Vector a = A - O; Vector b = B - O; Vector c = C - O;

float Denominator = 2.0f * Matrix::det(a.x, a.y, a.z, b.x, b.y, b.z, c.x, c.y, c.z);

Vector o = ((c^2) * (a % b) + (b^2) * (c % a) + (a^2) * (b % c)) / Denominator;

radius = Vector::N(o) + radiusEpsilon; center = O + o; }

Sphere &Sphere::operator=(const Sphere &S) { radius = S.radius; center = S.center;

return *this; }

float Sphere::d(const Point &P) const { return Point::d(P, center) - radius; }

float Sphere::d2(const Point &P) const { return Point::d2(P, center) - radius*radius; }

float Sphere::d(const Sphere &S, const Point &P) { return Point::d(P, S.center) - S.radius; }

float Sphere::d(const Point &P, const Sphere &S) { return Point::d(P, S.center) - S.radius; }

float Sphere::d2(const Sphere &S, const Point &P) { return Point::d2(P, S.center) - S.radius*S.radius; }

float Sphere::d2(const Point &P, const Sphere &S) { return Point::d2(P, S.center) - S.radius*S.radius; }

Sphere Sphere::recurseMini(Point *P[], unsigned int p, unsigned int b) { Sphere MB;

switch(b) { case 0: MB = Sphere(); break; case 1: MB = Sphere(*P[-1]); break; case 2: MB = Sphere(*P[-1], *P[-2]); break; case 3: MB = Sphere(*P[-1], *P[-2], *P[-3]); break; case 4: MB = Sphere(*P[-1], *P[-2], *P[-3], *P[-4]); return MB; }

for(unsigned int i = 0; i < p; i++) if(MB.d2(*P[i]) > 0) // Signed square distance to sphere { for(unsigned int j = i; j > 0; j--) { Point *T = P[j]; P[j] = P[j - 1]; P[j - 1] = T; }

MB = recurseMini(P + 1, i, b + 1); }

return MB; }

Sphere Sphere::miniBall(Point P[], unsigned int p) { Point **L = new Point*[p];

for(unsigned int i = 0; i < p; i++) L[i] = &P[i];

Sphere MB = recurseMini(L, p);

delete[] L;

return MB; }

Sphere Sphere::smallBall(Point P[], unsigned int p) { Vector center; float radius = -1;

if(p > 0) { center = 0;

for(unsigned int i = 0; i < p; i++) center += P[i];

center /= (float)p;

for(i = 0; i < p; i++) { float d2 = Point::d2(P[i], center); if(d2 > radius) radius = d2; }

radius = sqrtf(radius) + radiusEpsilon; }

return Sphere(center, radius); } }

Currently browsing [MiniBall.zip] (54,486 bytes) - [Timer.cpp] - (536 bytes)

#include "Timer.h"

#include <windows.h>

namespace MiniBall { Timer::Timer() { }

double Timer::seconds() { __int64 currentTime; __int64 frequency;

QueryPerformanceFrequency((LARGE_INTEGER*)&frequency); QueryPerformanceCounter((LARGE_INTEGER*)¤tTime);

return (double)currentTime / (double)frequency; }

double Timer::ticks() { return 4800.0 * seconds(); }

void Timer::pause(double p) { double stopTime = Timer::seconds() + p;

while(Timer::seconds() < stopTime); } }

Currently browsing [MiniBall.zip] (54,486 bytes) - [Tuple.cpp] - (46 bytes)

#include "Tuple.h"

namespace MiniBall { }

Currently browsing [MiniBall.zip] (54,486 bytes) - [Vector.cpp] - (1,982 bytes)

#include "Vector.h"

#include <math.h>

namespace MiniBall { Vector Vector::operator+() const { return *this; }

Vector Vector::operator-() const { return Vector(-x, -y, -z); }

Vector &Vector::operator+=(const Vector &v) { x += v.x; y += v.y; z += v.z;

return *this; }

Vector &Vector::operator-=(const Vector &v) { x -= v.x; y -= v.y; z -= v.z;

return *this; }

Vector &Vector::operator*=(float s) { x *= s; y *= s; z *= s;

return *this; }

Vector &Vector::operator/=(float s) { float r = 1.0f / s;

return *this *= r; }

bool operator==(const Vector &U, const Vector &v) { if(U.x == v.x && U.y == v.y && U.z == v.z) return true; else return false; }

bool operator!=(const Vector &U, const Vector &v) { if(U.x != v.x || U.y != v.y || U.z != v.z) return true; else return false; }

Vector operator+(const Vector &u, const Vector &v) { return Vector(u.x + v.x, u.y + v.y, u.z + v.z); }

Vector operator-(const Vector &u, const Vector &v) { return Vector(u.x - v.x, u.y - v.y, u.z - v.z); }

float operator*(const Vector &u, const Vector &v) { return u.x * v.x + u.y * v.y + u.z * v.z; }

Vector operator*(float s, const Vector &v) { return Vector(s * v.x, s * v.y, s * v.z); }

Vector operator*(const Vector &v, float s) { return Vector(v.x * s, v.y * s, v.z * s); }

Vector operator/(const Vector &v, float s) { float r = 1.0f / s;

return Vector(v.x * r, v.y * r, v.z * r); }

float operator^(const Vector &u, const Vector &v) { return acosf(u / Vector::N(u) * v / Vector::N(v)); }

Vector operator%(const Vector &u, const Vector &v) { return Vector(u.y * v.z - u.z * v.y, u.z * v.x - u.x * v.z, u.x * v.y - u.y * v.x); }

float Vector::N(const Vector &v) { return sqrtf(v^2); }

float Vector::N2(const Vector &v) { return v^2; } }

Currently browsing [MiniBall.zip] (54,486 bytes) - [ViewFrustum.cpp] - (1,300 bytes)

#include "ViewFrustum.h"

#include "Matrix.h"

namespace MiniBall { ViewFrustum::ViewFrustum() { }

ViewFrustum::ViewFrustum(const ViewFrustum &viewFrustum) { nearPlane = viewFrustum.nearPlane; farPlane = viewFrustum.farPlane; leftPlane = viewFrustum.leftPlane; rightPlane = viewFrustum.rightPlane; topPlane = viewFrustum.topPlane; bottomPlane = viewFrustum.bottomPlane; }

ViewFrustum::~ViewFrustum() { }

ViewFrustum &ViewFrustum::operator=(const ViewFrustum &viewFrustum) { nearPlane = viewFrustum.nearPlane; farPlane = viewFrustum.farPlane; leftPlane = viewFrustum.leftPlane; rightPlane = viewFrustum.rightPlane; topPlane = viewFrustum.topPlane; bottomPlane = viewFrustum.bottomPlane;

return *this; }

ViewFrustum operator*(const Matrix &M, const ViewFrustum &viewFrustum) { ViewFrustum transformedFrustum;

transformedFrustum.nearPlane = M * viewFrustum.nearPlane; transformedFrustum.farPlane = M * viewFrustum.farPlane; transformedFrustum.leftPlane = M * viewFrustum.leftPlane; transformedFrustum.rightPlane = M * viewFrustum.rightPlane; transformedFrustum.topPlane = M * viewFrustum.topPlane; transformedFrustum.bottomPlane = M * viewFrustum.bottomPlane;

return transformedFrustum; } }

Currently browsing [MiniBall.zip] (54,486 bytes) - [BoundingSphere.h] - (758 bytes)

#ifndef BOUNDINGSPHERE_H
#define BOUNDINGSPHERE_H

#include "BoundingVolume.h" #include "Sphere.h"

namespace MiniBall { class Matrix; class ViewFrustum; class Viewport; class Plane;

class BoundingSphere : public BoundingVolume, public Sphere { public: BoundingSphere(); BoundingSphere(const BoundingSphere &B); BoundingSphere(const Sphere &S);

~BoundingSphere();

BoundingSphere &operator=(const BoundingSphere &B); BoundingSphere &operator=(const Sphere &B);

bool visible(const ViewFrustum &Frustum) const; bool intersect(const Plane &ClippingPlane) const;

bool frontSide(const Plane &ClippingPlane) const; bool backSide(const Plane &ClippingPlane) const; }; }

#endif // BOUNDINGSPHERE_H

The zip file viewer built into the Developer Toolbox made use of the zlib library, as well as the zlibdll source additions.

 

Copyright 1999-2008 (C) FLIPCODE.COM and/or the original content author(s). All rights reserved.
Please read our Terms, Conditions, and Privacy information.