src/2geom/point.cpp

   1 #include <2geom/point.h>
   2 #include <assert.h>
   3 #include <2geom/coord.h>
   4 #include <2geom/isnan.h> //temporary fix for isnan()
   5 #include <2geom/transforms.h>
   6
   7 namespace Geom {
   8
   9 /** \brief Scales this vector to make it a unit vector (within rounding error).
  10  *
  11  *  The current version tries to handle infinite coordinates gracefully,
  12  *  but it's not clear that any callers need that.
  13  *
  14  *  \pre \f$this \neq (0, 0)\f$
  15  *  \pre Neither component is NaN.
  16  *  \post \f$-\epsilon<\left|this\right|-1<\epsilon\f$
  17  */
  18 void Point::normalize() {
  19     double len = hypot(_pt[0], _pt[1]);
  20     if(len == 0) return;
  21     if(IS_NAN(len)) return;
  22     static double const inf = 1e400;
  23     if(len != inf) {
  24         *this /= len;
  25     } else {
  26         unsigned n_inf_coords = 0;
  27         /* Delay updating pt in case neither coord is infinite. */
  28         Point tmp;
  29         for ( unsigned i = 0 ; i < 2 ; ++i ) {
  30             if ( _pt[i] == inf ) {
  31                 ++n_inf_coords;
  32                 tmp[i] = 1.0;
  33             } else if ( _pt[i] == -inf ) {
  34                 ++n_inf_coords;
  35                 tmp[i] = -1.0;
  36             } else {
  37                 tmp[i] = 0.0;
  38             }
  39         }
  40         switch (n_inf_coords) {
  41             case 0: {
  42                 /* Can happen if both coords are near +/-DBL_MAX. */
  43                 *this /= 4.0;
  44                 len = hypot(_pt[0], _pt[1]);
  45                 assert(len != inf);
  46                 *this /= len;
  47                 break;
  48             }
  49             case 1: {
  50                 *this = tmp;
  51                 break;
  52             }
  53             case 2: {
  54                 *this = tmp * sqrt(0.5);
  55                 break;
  56             }
  57         }
  58     }
  59 }
  60
  61 /** Compute the L1 norm, or manhattan distance, of \a p. */
  62 Coord L1(Point const &p) {
  63     Coord d = 0;
  64     for ( int i = 0 ; i < 2 ; i++ ) {
  65         d += fabs(p[i]);
  66     }
  67     return d;
  68 }
  69
  70 /** Compute the L infinity, or maximum, norm of \a p. */
  71 Coord LInfty(Point const &p) {
  72     Coord const a(fabs(p[0]));
  73     Coord const b(fabs(p[1]));
  74     return ( a < b || IS_NAN(b)
  75              ? b
  76              : a );
  77 }
  78
  79 /** Returns true iff p is a zero vector, i.e.\ Point(0, 0).
  80  *
  81  *  (NaN is considered non-zero.)
  82  */
  83 bool
  84 is_zero(Point const &p)
  85 {
  86     return ( p[0] == 0 &&
  87              p[1] == 0   );
  88 }
  89
  90 bool
  91 is_unit_vector(Point const &p)
  92 {
  93     return fabs(1.0 - L2(p)) <= 1e-4;
  94     /* The tolerance of 1e-4 is somewhat arbitrary.  Point::normalize is believed to return
  95        points well within this tolerance.  I'm not aware of any callers that want a small
  96        tolerance; most callers would be ok with a tolerance of 0.25. */
  97 }
  98
  99 Coord atan2(Point const p) {
 100     return std::atan2(p[Y], p[X]);
 101 }
 102
 103 /** compute the angle turning from a to b.  This should give \f$\pi/2\f$ for angle_between(a, rot90(a));
 104  * This works by projecting b onto the basis defined by a, rot90(a)
 105  */
 106 Coord angle_between(Point const a, Point const b) {
 107     return std::atan2(cross(b,a), dot(b,a));
 108 }
 109
 110
 111
 112 /** Returns a version of \a a scaled to be a unit vector (within rounding error).
 113  *
 114  *  The current version tries to handle infinite coordinates gracefully,
 115  *  but it's not clear that any callers need that.
 116  *
 117  *  \pre a != Point(0, 0).
 118  *  \pre Neither coordinate is NaN.
 119  *  \post L2(ret) very near 1.0.
 120  */
 121 Point unit_vector(Point const &a)
 122 {
 123     Point ret(a);
 124     ret.normalize();
 125     return ret;
 126 }
 127
 128 Point abs(Point const &b)
 129 {
 130     Point ret;
 131     for ( int i = 0 ; i < 2 ; i++ ) {
 132         ret[i] = fabs(b[i]);
 133     }
 134     return ret;
 135 }
 136
 137 Point operator*(Point const &v, Matrix const &m) {
 138     Point ret;
 139     for(int i = 0; i < 2; i++) {
 140         ret[i] = v[X] * m[i] + v[Y] * m[i + 2] + m[i + 4];
 141     }
 142     return ret;
 143 }
 144
 145 Point operator/(Point const &p, Matrix const &m) { return p * m.inverse(); }
 146
 147 Point &Point::operator*=(Matrix const &m)
 148 {
 149     *this = *this * m;
 150     return *this;
 151 }
 152
 153 Point constrain_angle(Point const &A, Point const &B, unsigned int n, Point const &dir)
 154 {
 155     // for special cases we could perhaps use explicit testing (which might be faster)
 156     if (n == 0.0) {
 157         return B;
 158     }
 159     Point diff(B - A);
 160     double angle = -angle_between(diff, dir);
 161     double k = round(angle * (double)n / (2.0*M_PI));
 162     return A + dir * Rotate(k * 2.0 * M_PI / (double)n) * L2(diff);
 163 }
 164
 165 }  //namespace Geom
 166
 167 /*
 168   Local Variables:
 169   mode:c++
 170   c-file-style:"stroustrup"
 171   c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))
 172   indent-tabs-mode:nil
 173   fill-column:99
 174   End:
 175 */
 176 // vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:fileencoding=utf-8:textwidth=99 :