// Copyright 2003, softSurfer (www.softsurfer.com) // This code may be freely used and modified for any purpose // providing that this copyright notice is included with it. // SoftSurfer makes no warranty for this code, and cannot be held // liable for any real or imagined damage resulting from its use. // Users of this code must verify correctness for their application. //-------------------------------------------------------------------- // common.h - common definition header file //-------------------------------------------------------------------- #ifndef SS_Common_H #define SS_Common_H #include <iostream.h> enum {FALSE=0, TRUE=1, ERROR=(-1)}; // Error codes enum Error { Enot, // no error Edim, // error: dim of point invalid for operation Esum // error: sum not affine (cooefs add to 1) }; // utility macros #define abs(x) ((x) >= 0 ? x : -(x)); #define min(x,y) ((x) < (y) ? (x) : (y)); #define max(x,y) ((x) > (y) ? (x) : (y)); #endif SS_Common_H //-------------------------------------------------------------------- // point.h - Point class definition //-------------------------------------------------------------------- #ifndef SS_Point_H #define SS_Point_H #include "common.h" class Point { friend class Vector; private: int dimn; // # coords (1, 2, or 3 max here) Error err; // error indicator public: double x, y, z; // z=0 for 2D, y=z=0 for 1D //---------------------------------------------------------- // Lots of Constructors (add more as needed) Point() { dimn=3; x=y=z=0; err=Enot; } // 1D Point Point( int a) { dimn=1; x=a; y=z=0; err=Enot; } Point( double a) { dimn=1; x=a; y=z=0; err=Enot; } // 2D Point Point( int a, int b) { dimn=2; x=a; y=b; z=0; err=Enot; } Point( double a, double b) { dimn=2; x=a; y=b; z=0; err=Enot; } // 3D Point Point( int a, int b, int c) { dimn=3; x=a; y=b; z=c; err=Enot; } Point( double a, double b, double c) { dimn=3; x=a; y=b; z=c; err=Enot; } // n-dim Point Point( int n, int a[]); Point( int n, double a[]); // Destructor ~Point() {}; //---------------------------------------------------------- // Input/Output streams friend istream& operator>>( istream&, Point&); friend ostream& operator<<( ostream&, Point); //---------------------------------------------------------- // Assignment "=": use the default to copy all members int dim() { return dimn; } // get dimension int setdim( int); // set new dimension //---------------------------------------------------------- // Comparison (dimension must match, or not) int operator==( Point); int operator!=( Point); //---------------------------------------------------------- // Point and Vector Operations (always valid) Vector operator-( Point); // Vector difference Point operator+( Vector); // +translate Point operator-( Vector); // -translate Point& operator+=( Vector); // inc translate Point& operator-=( Vector); // dec translate //---------------------------------------------------------- // Point Scalar Operations (convenient but often illegal) // using any type of scalar (int, float, or double) // are not valid for points in general, // unless they are 'affine' as coeffs of // a sum in which all the coeffs add to 1, // such as: the sum (a*P + b*Q) with (a+b == 1). // The programmer must enforce this (if they want to). // Scalar Multiplication friend Point operator*( int, Point); friend Point operator*( double, Point); friend Point operator*( Point, int); friend Point operator*( Point, double); // Scalar Division friend Point operator/( Point, int); friend Point operator/( Point, double); //---------------------------------------------------------- // Point Addition (also convenient but often illegal) // is not valid unless part of an affine sum. // The programmer must enforce this (if they want to). friend Point operator+( Point, Point); // add points // Affine Sum // Returns weighted sum, even when not affine, but... // Tests if coeffs add to 1. If not, sets: err = Esum. friend Point asum( int, int[], Point[]); friend Point asum( int, double[], Point[]); //---------------------------------------------------------- // Point Relations friend double d( Point, Point); // Distance friend double d2( Point, Point); // Distance^2 double isLeft( Point, Point); // 2D only //---------------------------------------------------------- // Error Handling void clerr() { err = Enot;} // clear error int geterr() { return err;} // get error char* errstr(); // return error string }; #endif SS_Point_H //-------------------------------------------------------------------- // vector.h - Vector class definition //-------------------------------------------------------------------- #ifndef SS_Vector_H #define SS_Vector_H #include "common.h" class Vector : public Point { public: // Constructors same as Point class Vector() : Point() {}; Vector( int a) : Point(a) {}; Vector( double a) : Point(a) {}; Vector( int a, int b) : Point(a,b) {}; Vector( double a, double b) : Point(a,b) {}; Vector( int a, int b, int c) : Point(a,b,c) {}; Vector( double a, double b, double c) : Point(a,b,c) {}; Vector( int n, int a[]) : Point(n,a) {}; Vector( int n, double a[]) : Point(n,a) {}; ~Vector() {}; //---------------------------------------------------------- // IO streams and Comparisons: inherit from Point class //---------------------------------------------------------- // Vector Unary Operations Vector operator-(); // unary minus Vector operator~(); // unary 2D perp operator //---------------------------------------------------------- // Scalar Multiplication friend Vector operator*( int, Vector); friend Vector operator*( double, Vector); friend Vector operator*( Vector, int); friend Vector operator*( Vector, double); // Scalar Division friend Vector operator/( Vector, int); friend Vector operator/( Vector, double); //---------------------------------------------------------- // Vector Arithmetic Operations Vector operator+( Vector); // vector add Vector operator-( Vector); // vector subtract double operator*( Vector); // inner dot product double operator|( Vector); // 2D exterior perp product Vector operator^( Vector); // 3D exterior cross product Vector& operator*=( double); // vector scalar mult Vector& operator/=( double); // vector scalar div Vector& operator+=( Vector); // vector increment Vector& operator-=( Vector); // vector decrement Vector& operator^=( Vector); // 3D exterior cross product //---------------------------------------------------------- // Vector Properties double len() { // vector length return sqrt(x*x + y*y + z*z); } double len2() { // vector length squared (faster) return (x*x + y*y + z*z); } //---------------------------------------------------------- // Special Operations void normalize(); // convert vector to unit length friend Vector sum( int, int[], Vector[]); // vector sum friend Vector sum( int, double[], Vector[]); // vector sum }; #endif SS_Vector_H point.c: //================================================================== // Copyright 2003, softSurfer (www.softsurfer.com) // This code may be freely used and modified for any purpose // providing that this copyright notice is included with it. // SoftSurfer makes no warranty for this code, and cannot be held // liable for any real or imagined damage resulting from it's use. // Users of this code must verify correctness for their application. //================================================================== #include "point.h" #include "vector.h" //================================================================== // Point Class Methods //================================================================== //------------------------------------------------------------------ // Constructors (add more as needed) //------------------------------------------------------------------ // n-dim Point Point::Point( int n, int a[]) { x = y = z = 0; err = Enot; switch (dimn = n) { case 3: z = a[2]; case 2: y = a[1]; case 1: x = a[0]; break; default: err=Edim; } } Point::Point( int n, double a[]) { x = y = z = 0.0; err = Enot; switch (dimn = n) { case 3: z = a[2]; case 2: y = a[1]; case 1: x = a[0]; break; default: err=Edim; } } //------------------------------------------------------------------ // IO streams //------------------------------------------------------------------ // Read input Point format: "(%f)", "(%f, %f)", or "(%f, %f, %f)" istream& operator>>( istream& input, Point& P) { char c; input >> c; // skip '(' input >> P.x; input >> c; if (c == ')') { P.setdim(1); // 1D coord return input; } // else // skip ',' input >> P.y; input >> c; if (c == ')') { P.setdim(2); // 2D coord return input; } // else // skip ',' input >> P.z; P.setdim(3); // 3D coord input >> c; // skip ')' return input; } // Write output Point in format: "(%f)", "(%f, %f)", or "(%f, %f, %f)" ostream& operator<<( ostream& output, Point P) { switch (P.dim()) { case 1: output << "(" << P.x << ")"; break; case 2: output << "(" << P.x << ", " << P.y << ")"; break; case 3: output << "(" << P.x << ", " << P.y << ", " << P.z << ")"; break; default: output << "Error: P.dim = " << P.dim(); } return output; } //------------------------------------------------------------------ // Assign (set) dimension //------------------------------------------------------------------ int Point::setdim( int n) { switch (n) { case 1: y = 0; case 2: z = 0; case 3: return dimn = n; default: // out of range value err = Edim; // just flag the error return ERROR; } } //------------------------------------------------------------------ // Comparison (note: dimension must compare) //------------------------------------------------------------------ int Point::operator==( Point Q) { if (dimn != Q.dim()) return FALSE; switch (dimn) { case 1: return (x==Q.x); case 2: return (x==Q.x && y==Q.y); case 3: default: return (x==Q.x && y==Q.y && z==Q.z); } } int Point::operator!=( Point Q) { if (dimn != Q.dim()) return TRUE; switch (dimn) { case 1: return (x!=Q.x); case 2: return (x!=Q.x || y!=Q.y); case 3: default: return (x!=Q.x || y!=Q.y || z!=Q.z); } } //------------------------------------------------------------------ // Point Vector Operations //------------------------------------------------------------------ Vector Point::operator-( Point Q) // Vector diff of Points { Vector v; v.x = x - Q.x; v.y = y - Q.y; v.z = z - Q.z; v.dimn = max( dimn, Q.dim()); return v; } Point Point::operator+( Vector v) // +ve translation { Point P; P.x = x + v.x; P.y = y + v.y; P.z = z + v.z; P.dimn = max( dimn, v.dim()); return P; } Point Point::operator-( Vector v) // -ve translation { Point P; P.x = x - v.x; P.y = y - v.y; P.z = z - v.z; P.dimn = max( dimn, v.dim()); return P; } Point& Point::operator+=( Vector v) // +ve translation { x += v.x; y += v.y; z += v.z; dimn = max( dimn, v.dim()); return *this; } Point& Point::operator-=( Vector v) // -ve translation { x -= v.x; y -= v.y; z -= v.z; dimn = max( dimn, v.dim()); return *this; } //------------------------------------------------------------------ // Point Scalar Operations (convenient but often illegal) // are not valid for points in general, // unless they are 'affine' as coeffs of // a sum in which all the coeffs add to 1, // such as: the sum (a*P + b*Q) with (a+b == 1). // The programmer must enforce this (if they want to). //------------------------------------------------------------------ Point operator*( int c, Point Q) { Point P; P.x = c * Q.x; P.y = c * Q.y; P.z = c * Q.z; P.dimn = Q.dim(); return P; } Point operator*( double c, Point Q) { Point P; P.x = c * Q.x; P.y = c * Q.y; P.z = c * Q.z; P.dimn = Q.dim(); return P; } Point operator*( Point Q, int c) { Point P; P.x = c * Q.x; P.y = c * Q.y; P.z = c * Q.z; P.dimn = Q.dim(); return P; } Point operator*( Point Q, double c) { Point P; P.x = c * Q.x; P.y = c * Q.y; P.z = c * Q.z; P.dimn = Q.dim(); return P; } Point operator/( Point Q, int c) { Point P; P.x = Q.x / c; P.y = Q.y / c; P.z = Q.z / c; P.dimn = Q.dim(); return P; } Point operator/( Point Q, double c) { Point P; P.x = Q.x / c; P.y = Q.y / c; P.z = Q.z / c; P.dimn = Q.dim(); return P; } //------------------------------------------------------------------ // Point Addition (also convenient but often illegal) // is not valid unless part of an affine sum. // The programmer must enforce this (if they want to). //------------------------------------------------------------------ Point operator+( Point Q, Point R) { Point P; P.x = Q.x + R.x; P.y = Q.y + R.y; P.z = Q.z + R.z; P.dimn = max( Q.dim(), R.dim()); return P; } //------------------------------------------------------------------ // Affine Sums // Returns weighted sum, even when not affine, but... // Tests if coeffs add to 1. If not, sets: err = Esum. //------------------------------------------------------------------ Point asum( int n, int c[], Point Q[]) { int maxd = 0; int cs = 0; Point P; for (int i=0; i<n; i++) { cs += c[i]; if (Q[i].dim() > maxd) maxd = Q[i].dim(); } if (cs != 1) // not an affine sum P.err = Esum; // flag error, but compute sum anyway for (int i=0; i<n; i++) { P.x += c[i] * Q[i].x; P.y += c[i] * Q[i].y; P.z += c[i] * Q[i].z; } P.dimn = maxd; return P; } Point asum( int n, double c[], Point Q[]) { int maxd = 0; double cs = 0.0; Point P; for (int i=0; i<n; i++) { cs += c[i]; if (Q[i].dim() > maxd) maxd = Q[i].dim(); } if (cs != 1) // not an affine sum P.err = Esum; // flag error, but compute sum anyway for (int i=0; i<n; i++) { P.x += c[i] * Q[i].x; P.y += c[i] * Q[i].y; P.z += c[i] * Q[i].z; } P.dimn = maxd; return P; } //------------------------------------------------------------------ // Distance between Points //------------------------------------------------------------------ double d( Point P, Point Q) { // Euclidean distance double dx = P.x - Q.x; double dy = P.y - Q.y; double dz = P.z - Q.z; return sqrt(dx*dx + dy*dy + dz*dz); } double d2( Point P, Point Q) { // squared distance (more efficient) double dx = P.x - Q.x; double dy = P.y - Q.y; double dz = P.z - Q.z; return (dx*dx + dy*dy + dz*dz); } //------------------------------------------------------------------ // Sidedness of a Point wrt a directed line P1->P2 // - makes sense in 2D only //------------------------------------------------------------------ double Point::isLeft( Point P1, Point P2) { if (dimn != 2 || P1.dim() != 2 || P2.dim() != 2) { err = Edim; // flag error, but compute anyway } return ((P1.x - x) * (P2.y - y) - (P2.x - x) * (P1.y - y)); } //------------------------------------------------------------------ // Error Routines //------------------------------------------------------------------ char* Point::errstr() { // return error string switch (err) { case Enot: return "no error"; case Edim: return "error: invalid dimension for operation"; case Esum: return "error: Point sum is not affine"; default: return "error: unknown err value"; } } vector.c //================================================================== // Copyright 2003, softSurfer (www.softsurfer.com) // This code may be freely used and modified for any purpose // providing that this copyright notice is included with it. // SoftSurfer makes no warranty for this code, and cannot be held // liable for any real or imagined damage resulting from it's use. // Users of this code must verify correctness for their application. //================================================================== #include "point.h" #include "vector.h" //================================================================== // Vector Class Methods //================================================================== //------------------------------------------------------------------ // Unary Ops //------------------------------------------------------------------ // unary minus Vector Vector::operator-() { Vector v; v.x = -x; v.y = -y; v.z = -z; v.dimn = dimn; return v; } // unary 2D perp operator Vector Vector::operator~() { if (dimn != 2) err = Edim; // and continue anyway Vector v; v.x = -y; v.y = x; v.z = z; v.dimn = dimn; return v; } //------------------------------------------------------------------ // Scalar Ops //------------------------------------------------------------------ // Scalar Multiplies Vector operator*( int c, Vector w ) { Vector v; v.x = c * w.x; v.y = c * w.y; v.z = c * w.z; v.setdim( w.dim()); return v; } Vector operator*( double c, Vector w ) { Vector v; v.x = c * w.x; v.y = c * w.y; v.z = c * w.z; v.setdim( w.dim()); return v; } Vector operator*( Vector w, int c ) { Vector v; v.x = c * w.x; v.y = c * w.y; v.z = c * w.z; v.setdim( w.dim()); return v; } Vector operator*( Vector w, double c ) { Vector v; v.x = c * w.x; v.y = c * w.y; v.z = c * w.z; v.setdim( w.dim()); return v; } // Scalar Divides Vector operator/( Vector w, int c ) { Vector v; v.x = w.x / c; v.y = w.y / c; v.z = w.z / c; v.setdim( w.dim()); return v; } Vector operator/( Vector w, double c ) { Vector v; v.x = w.x / c; v.y = w.y / c; v.z = w.z / c; v.setdim( w.dim()); return v; } //------------------------------------------------------------------ // Arithmetic Ops //------------------------------------------------------------------ Vector Vector::operator+( Vector w ) { Vector v; v.x = x + w.x; v.y = y + w.y; v.z = z + w.z; v.dimn = max( dimn, w.dim()); return v; } Vector Vector::operator-( Vector w ) { Vector v; v.x = x - w.x; v.y = y - w.y; v.z = z - w.z; v.dimn = max( dimn, w.dim()); return v; } //------------------------------------------------------------------ // Products //------------------------------------------------------------------ // Inner Dot Product double Vector::operator*( Vector w ) { return (x * w.x + y * w.y + z * w.z); } // 2D Exterior Perp Product double Vector::operator|( Vector w ) { if (dimn != 2) err = Edim; // and continue anyway return (x * w.y - y * w.x); } // 3D Exterior Cross Product Vector Vector::operator^( Vector w ) { Vector v; v.x = y * w.z - z * w.y; v.y = z * w.x - x * w.z; v.z = x * w.y - y * w.x; v.dimn = 3; return v; } //------------------------------------------------------------------ // Shorthand Ops //------------------------------------------------------------------ Vector& Vector::operator*=( double c ) { // vector scalar mult x *= c; y *= c; z *= c; return *this; } Vector& Vector::operator/=( double c ) { // vector scalar div x /= c; y /= c; z /= c; return *this; } Vector& Vector::operator+=( Vector w ) { // vector increment x += w.x; y += w.y; z += w.z; dimn = max(dimn, w.dim()); return *this; } Vector& Vector::operator-=( Vector w ) { // vector decrement x -= w.x; y -= w.y; z -= w.z; dimn = max(dimn, w.dim()); return *this; } Vector& Vector::operator^=( Vector w ) { // 3D exterior cross product double ox=x, oy=y, oz=z; x = oy * w.z - oz * w.y; y = oz * w.x - ox * w.z; z = ox * w.y - oy * w.x; dimn = 3; return *this; } //------------------------------------------------------------------ // Special Operations //------------------------------------------------------------------ void Vector::normalize() { // convert to unit length double ln = sqrt( x*x + y*y + z*z ); if (ln == 0) return; // do nothing for nothing x /= ln; y /= ln; z /= ln; } Vector sum( int n, int c[], Vector w[] ) { // vector sum int maxd = 0; Vector v; for (int i=0; i<n; i++) { if (w[i].dim() > maxd) maxd = w[i].dim(); } v.setdim( maxd); for (int i=0; i<n; i++) { v.x += c[i] * w[i].x; v.y += c[i] * w[i].y; v.z += c[i] * w[i].z; } return v; } Vector sum( int n, double c[], Vector w[] ) { // vector sum int maxd = 0; Vector v; for (int i=0; i<n; i++) { if (w[i].dim() > maxd) maxd = w[i].dim(); } v.setdim( maxd); for (int i=0; i<n; i++) { v.x += c[i] * w[i].x; v.y += c[i] * w[i].y; v.z += c[i] * w[i].z; } return v; }