d3/point.h

   1 // Copyright 2002, 2004 David Hilvert <dhilvert@auricle.dyndns.org>,
   2 //                                    <dhilvert@ugcs.caltech.edu>
   3
   4 /*  This file is part of the Anti-Lamenessing Engine.
   5
   6     The Anti-Lamenessing Engine is free software; you can redistribute it and/or modify
   7     it under the terms of the GNU General Public License as published by
   8     the Free Software Foundation; either version 3 of the License, or
   9     (at your option) any later version.
  10
  11     The Anti-Lamenessing Engine is distributed in the hope that it will be useful,
  12     but WITHOUT ANY WARRANTY; without even the implied warranty of
  13     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  14     GNU General Public License for more details.
  15
  16     You should have received a copy of the GNU General Public License
  17     along with the Anti-Lamenessing Engine; if not, write to the Free Software
  18     Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
  19 */
  20
  21 #ifndef __d3point_h__
  22 #define __d3point_h__
  23
  24 /*
  25  * Structure to describe a point in three dimensions.
  26  */
  27 class point {
  28 private:
  29         ale_pos x[3];
  30
  31 public:
  32         point() {
  33                 x[0] = 0;
  34                 x[1] = 0;
  35                 x[2] = 0;
  36         }
  37
  38         point(ale_pos x0, ale_pos x1, ale_pos x2) {
  39                 x[0] = x0;
  40                 x[1] = x1;
  41                 x[2] = x2;
  42         }
  43
  44         point(const point &p) {
  45                 x[0] = p[0];
  46                 x[1] = p[1];
  47                 x[2] = p[2];
  48         }
  49
  50         static point unit(int dimension) {
  51                 if (dimension == 0)
  52                         return point(1, 0, 0);
  53                 if (dimension == 1)
  54                         return point(0, 1, 0);
  55                 if (dimension == 2)
  56                         return point(0, 0, 1);
  57
  58                 assert(0);
  59         }
  60
  61         static point undefined() {
  62                 double a = 0;
  63
  64                 point p(0, 0, 0);
  65
  66                 return p / a;
  67         }
  68
  69         static point posinf() {
  70                 double a = +1;
  71                 double z = +0;
  72
  73                 a = a / z;
  74
  75                 assert (isinf(a) && a > 0);
  76
  77                 return point(a, a, a);
  78         }
  79
  80         static point neginf() {
  81                 point n = -posinf();
  82
  83                 assert (isinf(n[0]) && n[0] < 0);
  84
  85                 return n;
  86         }
  87
  88         void accumulate_max(point p) {
  89                 for (int d = 0; d < 3; d++)
  90                         if (p[d] > x[d])
  91                                 x[d] = p[d];
  92         }
  93
  94         void accumulate_min(point p) {
  95                 for (int d = 0; d < 3; d++)
  96                         if (p[d] < x[d])
  97                                 x[d] = p[d];
  98         }
  99
 100         int defined() const {
 101                 return (!isnan(x[0])
 102                      && !isnan(x[1])
 103                      && !isnan(x[2]));
 104         }
 105
 106         int finite() const {
 107                 return (::finite(x[0])
 108                      && ::finite(x[1])
 109                      && ::finite(x[2]));
 110         }
 111
 112         static int defined(const point &p) {
 113                 return p.defined();
 114         }
 115
 116         /*
 117          * Z-values of zero are almost never the right thing to do ...
 118          */
 119         point(d2::point p, ale_pos z = 0) {
 120                 x[0] = p[0];
 121                 x[1] = p[1];
 122                 x[2] = z;
 123         }
 124
 125         const ale_pos &operator[](int i) const {
 126                 assert (i >= 0);
 127                 assert (i < 3);
 128
 129                 return x[i];
 130         }
 131
 132         ale_pos &operator[](int i) {
 133                 assert (i >= 0);
 134                 assert (i < 3);
 135
 136                 return x[i];
 137         }
 138
 139         d2::point xy() const {
 140                 d2::point result;
 141
 142                 result[0] = x[0];
 143                 result[1] = x[1];
 144
 145                 return result;
 146         }
 147
 148         point operator+(point p) const {
 149                 return point(x[0] + p[0], x[1] + p[1], x[2] + p[2]);
 150         }
 151
 152         point operator-(point p) const {
 153                 return point(x[0] - p[0], x[1] - p[1], x[2] - p[2]);
 154         }
 155
 156         point operator-() const {
 157                 return point(-x[0], -x[1], -x[2]);
 158         }
 159
 160         point operator/(ale_pos r) const {
 161                 return point(x[0] / r, x[1] / r, x[2] / r);
 162         }
 163
 164         point operator /=(ale_pos r) {
 165                 x[0] /= r;
 166                 x[1] /= r;
 167                 x[2] /= r;
 168
 169                 return *this;
 170         }
 171
 172         point operator *=(ale_pos r) {
 173                 x[0] *= r;
 174                 x[1] *= r;
 175                 x[2] *= r;
 176
 177                 return *this;
 178         }
 179
 180         point operator +=(point p) {
 181                 x[0] += p[0];
 182                 x[1] += p[1];
 183                 x[2] += p[2];
 184
 185                 return *this;
 186         }
 187
 188         point operator -=(point p) {
 189                 x[0] -= p[0];
 190                 x[1] -= p[1];
 191                 x[2] -= p[2];
 192
 193                 return *this;
 194         }
 195
 196         int operator !=(point p) {
 197                 return (x[0] != p[0]
 198                      || x[1] != p[1]
 199                      || x[2] != p[2]);
 200         }
 201
 202         point mult(ale_pos d) const {
 203                 return point(x[0] * d, x[1] * d, x[2] * d);
 204         }
 205
 206         point operator*(point p) const {
 207                 return point(x[0] * p[0], x[1] * p[1], x[2] * p[2]);
 208         }
 209
 210         ale_pos normsq() const {
 211                 return x[0] * x[0] + x[1] * x[1] + x[2] * x[2];
 212         }
 213
 214         ale_pos norm() const {
 215                 return sqrt(normsq());
 216         }
 217
 218         point normalize() const {
 219                 return operator/(norm());
 220         }
 221
 222         ale_pos lengthtosq(point p) const {
 223                 point diff = operator-(p);
 224                 return diff.normsq();
 225         }
 226
 227         ale_pos lengthto(point p) const {
 228                 return sqrt(lengthtosq(p));
 229         }
 230
 231         ale_pos anglebetw(point p, point q) {
 232                 /*
 233                  * by the law of cosines, the cosine is equal to:
 234                  *
 235                  *      (lengthtosq(p) + lengthtosq(q) - p.lengthtosq(q))
 236                  *    / (2 * lengthto(p) * lengthto(q))
 237                  */
 238
 239                 ale_pos to_p = lengthtosq(p);
 240                 ale_pos to_q = lengthtosq(q);
 241
 242                 ale_pos cos_of = (double) (to_p + to_q - p.lengthtosq(q))
 243                                / (double) (2 * sqrt(to_p) * sqrt(to_q));
 244
 245                 /*
 246                  * XXX: is the fabs() required?
 247                  */
 248
 249                 return fabs(acos(cos_of));
 250         }
 251
 252
 253         /*
 254          * Determine the cross product
 255          */
 256         point xproduct(point p, point q) {
 257                 point pp = p;
 258                 point qq = q;
 259
 260                 pp -= *this;
 261                 qq -= *this;
 262
 263                 return point(pp[1] * qq[2] - pp[2] * qq[1],
 264                              pp[2] * qq[0] - pp[0] * qq[2],
 265                              pp[0] * qq[1] - pp[1] * qq[0]);
 266         }
 267
 268         /*
 269          * Determine the dot product
 270          */
 271         ale_pos dproduct(const point &p) {
 272                 return x[0] * p[0] + x[1] * p[1] + x[2] * p[2];
 273         }
 274
 275         /*
 276          * Determine whether the point is inside a given volume
 277          */
 278         int inside(const point &min, const point &max) {
 279                 for (int d = 0; d < 3; d++) {
 280                         if (min[d] > x[d])
 281                                 return 0;
 282                         if (max[d] < x[d])
 283                                 return 0;
 284                 }
 285
 286         }
 287 };
 288
 289 inline point operator*(const point &p, double d) {
 290         return p.mult(d);
 291 }
 292 inline point operator*(double d, const point &p) {
 293         return p.mult(d);
 294 }
 295 #endif