Shamusworld >> Repos - architektonas/blob - src/vector.cpp

   1 //\r
   2 // vector.cpp: Various structures used for 3 dimensional imaging\r
   3 //\r
   4 // by James Hammons\r
   5 // (C) 2006 Underground Software\r
   6 //\r
   7 // JLH = James L. Hammons <jlhamm@acm.org>\r
   8 //\r
   9 // WHO  WHEN        WHAT\r
  10 // ---  ----------  ------------------------------------------------------------\r
  11 // JLH  09/19/2006  Created this file\r
  12 // JLH  03/22/2011  Moved implementation of constructor from header to here\r
  13 // JLH  04/02/2011  Fixed divide-by-zero bug in Unit(), added Angle() function\r
  14 // JLH  08/04/2013  Added Parameter() function\r
  15 //\r
  16 \r
  17 #include "vector.h"\r
  18 \r
  19 #include <math.h>                                                               // For sqrt()\r
  20 #include "mathconstants.h"\r
  21 \r
  22 // Vector implementation\r
  23 \r
  24 Vector::Vector(double xx/*= 0*/, double yy/*= 0*/, double zz/*= 0*/): x(xx), y(yy), z(zz)\r
  25 {\r
  26 }\r
  27 \r
  28 \r
  29 Vector::Vector(Vector head, Vector tail): x(head.x - tail.x), y(head.y - tail.y), z(head.z - tail.z)\r
  30 {\r
  31 }\r
  32 \r
  33 \r
  34 Vector Vector::operator=(Vector const v)\r
  35 {\r
  36         x = v.x, y = v.y, z = v.z;\r
  37 \r
  38         return *this;\r
  39 }\r
  40 \r
  41 \r
  42 Vector Vector::operator+(Vector const v)\r
  43 {\r
  44         return Vector(x + v.x, y + v.y, z + v.z);\r
  45 }\r
  46 \r
  47 \r
  48 Vector Vector::operator-(Vector const v)\r
  49 {\r
  50         return Vector(x - v.x, y - v.y, z - v.z);\r
  51 }\r
  52 \r
  53 \r
  54 // Unary negation\r
  55 \r
  56 Vector Vector::operator-(void)\r
  57 {\r
  58         return Vector(-x, -y, -z);\r
  59 }\r
  60 \r
  61 \r
  62 // Vector x constant\r
  63 \r
  64 Vector Vector::operator*(double const v)\r
  65 {\r
  66         return Vector(x * v, y * v, z * v);\r
  67 }\r
  68 \r
  69 \r
  70 // Vector x constant\r
  71 \r
  72 Vector Vector::operator*(float const v)\r
  73 {\r
  74         return Vector(x * v, y * v, z * v);\r
  75 }\r
  76 \r
  77 \r
  78 // Vector / constant\r
  79 \r
  80 Vector Vector::operator/(double const v)\r
  81 {\r
  82         return Vector(x / v, y / v, z / v);\r
  83 }\r
  84 \r
  85 \r
  86 // Vector / constant\r
  87 \r
  88 Vector Vector::operator/(float const v)\r
  89 {\r
  90         return Vector(x / v, y / v, z / v);\r
  91 }\r
  92 \r
  93 \r
  94 // Vector (cross) product\r
  95 \r
  96 Vector Vector::operator*(Vector const v)\r
  97 {\r
  98         // a x b = [a2b3 - a3b2, a3b1 - a1b3, a1b2 - a2b1]\r
  99         return Vector((y * v.z) - (z * v.y), (z * v.x) - (x * v.z), (x * v.y) - (y * v.x));\r
 100 }\r
 101 \r
 102 \r
 103 // Dot product\r
 104 \r
 105 double Vector::Dot(Vector const v)\r
 106 {\r
 107         return (x * v.x) + (y * v.y) + (z * v.z);\r
 108 }\r
 109 \r
 110 \r
 111 // Vector x constant, self assigned\r
 112 \r
 113 Vector& Vector::operator*=(double const v)\r
 114 {\r
 115         x *= v, y *= v, z *= v;\r
 116 \r
 117         return *this;\r
 118 }\r
 119 \r
 120 \r
 121 // Vector / constant, self assigned\r
 122 \r
 123 Vector& Vector::operator/=(double const v)\r
 124 {\r
 125         x /= v, y /= v, z /= v;\r
 126 \r
 127         return *this;\r
 128 }\r
 129 \r
 130 // Vector + vector, self assigned\r
 131 \r
 132 Vector& Vector::operator+=(Vector const v)\r
 133 {\r
 134         x += v.x, y += v.y, z += v.z;\r
 135 \r
 136         return *this;\r
 137 }\r
 138 \r
 139 \r
 140 // Vector + constant, self assigned\r
 141 \r
 142 Vector& Vector::operator+=(double const v)\r
 143 {\r
 144         x += v, y += v, z += v;\r
 145 \r
 146         return *this;\r
 147 }\r
 148 \r
 149 \r
 150 // Vector - vector, self assigned\r
 151 \r
 152 Vector& Vector::operator-=(Vector const v)\r
 153 {\r
 154         x -= v.x, y -= v.y, z -= v.z;\r
 155 \r
 156         return *this;\r
 157 }\r
 158 \r
 159 \r
 160 // Vector - constant, self assigned\r
 161 \r
 162 Vector& Vector::operator-=(double const v)\r
 163 {\r
 164         x -= v, y -= v, z -= v;\r
 165 \r
 166         return *this;\r
 167 }\r
 168 \r
 169 \r
 170 // Check for equality\r
 171 bool Vector::operator==(Vector const v)\r
 172 {\r
 173         return (x == v.x && y == v.y && z == v.z ? true : false);\r
 174 }\r
 175 \r
 176 \r
 177 // Check for inequality\r
 178 bool Vector::operator!=(Vector const v)\r
 179 {\r
 180         return (x != v.x || y != v.y || z != v.z ? true : false);\r
 181 }\r
 182 \r
 183 \r
 184 Vector Vector::Unit(void)\r
 185 {\r
 186         double mag = Magnitude();\r
 187 \r
 188         // If the magnitude of the vector is zero, then the Unit vector is undefined...\r
 189         if (mag == 0)\r
 190                 return Vector(0, 0, 0);\r
 191 \r
 192         return Vector(x / mag, y / mag, z / mag);\r
 193 }\r
 194 \r
 195 \r
 196 double Vector::Magnitude(void)\r
 197 {\r
 198         return sqrt(x * x + y * y + z * z);\r
 199 }\r
 200 \r
 201 \r
 202 double Vector::Angle(void)\r
 203 {\r
 204         // acos returns a value between zero and PI, which means we don't know which\r
 205         // quadrant the angle is in... Though, if the y-coordinate of the vector is\r
 206         // negative, that means that the angle is in quadrants III - IV.\r
 207         double rawAngle = acos(Unit().x);\r
 208         double correctedAngle = (y < 0 ? (2.0 * PI) - rawAngle : rawAngle);\r
 209 \r
 210         return correctedAngle;\r
 211 }\r
 212 \r
 213 \r
 214 bool Vector::isZero(double epsilon/*= 1e-6*/)\r
 215 {\r
 216         return (fabs(x) < epsilon && fabs(y) < epsilon && fabs(z) < epsilon ? true : false);\r
 217 }\r
 218 \r
 219 \r
 220 // Class methods\r
 221 \r
 222 double Vector::Dot(Vector v1, Vector v2)\r
 223 {\r
 224         return (v1.x * v2.x) + (v1.y * v2.y) + (v1.z * v2.z);\r
 225 }\r
 226 \r
 227 \r
 228 double Vector::Magnitude(Vector v1, Vector v2)\r
 229 {\r
 230         double xx = v1.x - v2.x;\r
 231         double yy = v1.y - v2.y;\r
 232         double zz = v1.z - v2.z;\r
 233         return sqrt(xx * xx + yy * yy + zz * zz);\r
 234 }\r
 235 \r
 236 \r
 237 // Returns the parameter of a point in space to this vector. If the parameter\r
 238 // is between 0 and 1, the normal of the vector to the point is on the vector.\r
 239 double Vector::Parameter(Vector v1, Vector v2, Vector p)\r
 240 {\r
 241         // Geometric interpretation:\r
 242         // The parameterized point on the vector lineSegment is where the normal of\r
 243         // the lineSegment to the point intersects lineSegment. If the pp < 0, then\r
 244         // the perpendicular lies beyond the 1st endpoint. If pp > 1, then the\r
 245         // perpendicular lies beyond the 2nd endpoint.\r
 246 \r
 247         Vector lineSegment = v1 - v2;\r
 248         double magnitude = lineSegment.Magnitude();\r
 249         Vector pointSegment = p - v2;\r
 250         double t = lineSegment.Dot(pointSegment) / (magnitude * magnitude);\r
 251         return t;\r
 252 }\r
 253 \r
 254 \r
 255 // Return the normal to the linesegment formed by the passed in points.\r
 256 // (Not sure which is head or tail, or which hand the normal lies)\r
 257 /*static*/ Vector Vector::Normal(Vector v1, Vector v2)\r
 258 {\r
 259         Vector v = (v1 - v2).Unit();\r
 260         return Vector(-v.y, v.x);\r
 261 }\r
 262 \r