---

point3d.H

Go to the documentation of this file.

/*
 * Copyright 2000, Brown University, Providence, RI.
 * 
 *                         All Rights Reserved
 * 
 * Permission to use, copy, modify, and distribute this software and its
 * documentation for any purpose other than its incorporation into a
 * commercial product is hereby granted without fee, provided that the
 * above copyright notice appear in all copies and that both that
 * copyright notice and this permission notice appear in supporting
 * documentation, and that the name of Brown University not be used in
 * advertising or publicity pertaining to distribution of the software
 * without specific, written prior permission.
 * 
 * BROWN UNIVERSITY DISCLAIMS ALL WARRANTIES WITH REGARD TO THIS SOFTWARE,
 * INCLUDING ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR ANY
 * PARTICULAR PURPOSE.  IN NO EVENT SHALL BROWN UNIVERSITY BE LIABLE FOR
 * ANY SPECIAL, INDIRECT OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
 * WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN
 * ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF
 * OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.
 */
#ifndef JOT_POINT3D_H
#define JOT_POINT3D_H

// DESCRIPTION :
//    the '%' operator allows points to be added together

#include <limits.h>
#include "vector3d.H"
#include "std/config.H"
#include "nearest_pt.H"
#include "std/list.H"

template <class P, class V>
class _point3d
{
protected:
    Greal _x, _y, _z;

public:

    _point3d()                             : _x(0),    _y(0),    _z(0)   {}
    _point3d(Greal x, Greal y, Greal z) : _x(x),    _y(y),    _z(z)   {}

    void          set(Greal x, Greal y, Greal z)  { _x=x; _y=y; _z=z; }
    const Greal *data()            const{ return &_x; }

    P  operator  *(Greal s)        const{ return P(_x*s, _y*s, _z*s);}
    P  operator  /(Greal s)        const{ return P(_x/s, _y/s, _z/s);}
    // Point + Point
    P  operator  %(const P &p)      const{ return P(_x+p[0], _y+p[1], _z+p[2]);}
    P  operator  +(const V &v)      const{ return P(_x+v[0], _y+v[1], _z+v[2]);}
    V  operator  -(const P &p)      const{ return V(_x-p[0], _y-p[1], _z-p[2]);}
    P  operator  -(const V &v)      const{ return P(_x-v[0], _y-v[1], _z-v[2]);}
    P  operator  -()                const{ return P(-_x, -_y, -_z);}
    
    Greal  operator [](int index)  const{ return (&_x)[index]; }
    Greal& operator [](int index)       { return (&_x)[index]; }

    Greal  distSqrd   (const P &p) const{ return (_x-p._x)*(_x-p._x) + 
                                                     (_y-p._y)*(_y-p._y) + 
                                                     (_z-p._z)*(_z-p._z); }
    Greal  dist       (const P &p) const{ return sqrt(distSqrd(p));}

    bool    isEqual    (const P &p, Greal epsSqrd = epsAbsSqrdMath()) const 
                                         { return distSqrd(p) <= epsSqrd; }

    bool    operator ==(const P &p) const{ return _x==p._x&&_y==p._y&&_z==p._z;}
    bool    operator !=(const P &p) const{ return _x!=p._x||_y!=p._y||_z!=p._z;}
    // Point = Point + Point
    void    operator %=(const P &p)      { _x += p[0]; _y += p[1]; _z += p[2]; }
    // Point = Point + Vector
    void    operator +=(const V &v)      { _x += v[0]; _y += v[1]; _z += v[2]; }
    void    operator -=(const V &v)      { _x -= v[0]; _y -= v[1]; _z -= v[2]; }
    void    operator *=(Greal   s)       { _x *= s; _y *= s; _z *= s; }
    void    operator /=(Greal   s)       { _x /= s; _y /= s; _z /= s; }
}; // class point3d


/* ---------- inlined global functions using _point3d template ------ */

template <class P, class V>
inline ostream &
operator<<(ostream &os, const _point3d<P,V> &p) 
{ 
   return os << "< " << p[0] << " " << p[1] << " " << p[2] << " > ";
}


template <class P, class V>
inline istream &
operator>>(istream &is, _point3d<P,V> &p) 
{ 
   char dummy;
   return is >> dummy >> p[0] >> p[1] >> p[2] >> dummy;
}


template <class P, class V>
inline Greal 
det(const _point3d<P,V> &a, const _point3d<P,V> &b, const _point3d<P,V> &c)
{
   return (a[0] * (b[1]*c[2] - b[2]*c[1]) +
           a[1] * (b[2]*c[0] - b[0]*c[2]) +
           a[2] * (b[0]*c[1] - b[1]*c[0]));
} 


template <class P>
extern bool     areCoplanar(const P &, const P &, const P &, const P &);


// -----------------------------------------------------------------------
//
//   point3d_list -
//
// -----------------------------------------------------------------------
template <class L, class T, class P, class V>
class _point3d_list : public ARRAY<P> {
 protected:
   ARRAY<Greal>    _partial_length;

 public :
   _point3d_list(int max=16)       :ARRAY<P>(max),_partial_length(0) { }
   _point3d_list(const ARRAY<P> &p):ARRAY<P>(p),  _partial_length(0) { }
   _point3d_list(const ARRAY<P> &p, const T &t):ARRAY<P>(p),_partial_length(0) 
                                                           { xform(t); }

   void   xform   (const T &t);
   // computes everything about the point on the line closest to p.
   // returns the point, its distance to the line,
   // and the index of the point that starts its segment.
   // (that is, the point is on the line segment between
   // points [seg_index] and [seg_index+1])
   void   closest (const P &p, P &nearpt, Greal &neardist, int &segind) const{
      neardist = DBL_MAX;
     
      for (int i = 0; i < num()-1; i++) {
         const P a((*this)[i]);
         const P b((*this)[i+1]);
         const P npt = nearest_pt_to_line_seg(p,a,b);
         const Greal d = (npt-p).lengthSqrd();
         if (d < neardist) {
            neardist = d;
            nearpt   = npt;
            segind   = i;
         }
      }
      neardist = sqrt(neardist);
   }
   Greal closest (const P &p, P  &nearpt, int    &nearind) const {
      Greal d;
      int    i;
      closest(p, nearpt, d, i);

      P a((*this)[i]);
      P b((*this)[i+1]);
      if ((p-a).length() < (p-b).length())
           nearind = i;
      else nearind = i+1;
      return d;
   }
   P closest (const P &p) const {
      P nearpt(DBL_MAX,DBL_MAX,DBL_MAX);
      Greal dist;
      int    index;
      closest (p, nearpt, dist, index);
      return nearpt;
   }
   V      tan     (int i) const { P *ar = _array;
                                  const int n=num()-1;
                                  if (i<0 || i>n || n<1) return V();
                                  if (i==0) return (ar[1]-ar[  0]).normalize();
                                  if (i==n) return (ar[n]-ar[n-1]).normalize();
                                  return ((ar[i+1]-ar[i]).normalize() +
                                          (ar[i]-ar[i-1]).normalize()).
                                                                   normalize();
                                  }

   // Interpolation
   P      interpolate       (Greal s, V *tan=0, int*segp=0, Greal*tp=0) const{
                 int seg;
                 Greal t;
                 interpolate_length(s, seg, t);
                 const V v = _array[seg+1]-_array[seg];
                 if (tan)  *tan  = v.normalize();
                 if (segp) *segp = seg;
                 if (tp)   *tp   = t;

                 return _array[seg]+v*t;
                 }
         // assume update_length() has been called.
         // note if there are 0 or 1 points, _length will be 0, we return early.
         // note also, generally returns 0<=t<1,
         // except at the extreme right, returns t=1;
   void   interpolate_length(Greal s, int &seg, Greal &t)  const {

     const Greal val = s*length();
     if (val <= 0) { 
        seg = 0; 
        t   = 0; 
        return; 
     }
     if (val >= length()) { 
        seg = _num-2; 
        t = 1; 
        return; 
     }

     int l=0, r=_num-1, m;
     while ((m = (l+r) >> 1) != l) {
        if (val < _partial_length[m])       
         r = m;
        else l = m;
     }
     seg = m;
     t = (val-_partial_length[m])/(_partial_length[m+1]-_partial_length[m]);
      }
   Greal invert             (const P &pt) const {
      P      nearpt;
      Greal neardist;
      int    seg_index;

      closest(pt, nearpt, neardist, seg_index);
      return invert(pt, seg_index);
   }
   Greal invert             (const P &pt, int seg ) const {
      const V v     = _array[seg + 1] - _array[seg];
      const V vnorm = v.normalize();
      const P pp    = _array[seg] + ( (pt - _array[seg]) * vnorm ) * vnorm;

      Greal len = _partial_length[seg] + (pp-_array[seg]).length();
      return len/length();
   }

   // Parameterization / Lengths
   
   // Update partial lengths
   void   update_length      () { Greal len = 0;
                                  _partial_length.clear();
                                  if (_num > 0)
                                     _partial_length.realloc(_num);
                                  _partial_length += 0;
                                  for ( int i=1; i< _num; i++ ) {
                                     len += segment_length(i-1);
                                     _partial_length += len;
                                  }
                               }
   Greal partial_length     (int i) const { return i>=_partial_length.num() ? 
                                                   0 : _partial_length[i]; }
   // Partial lengths must be up to date
   Greal length             ()      const { return _partial_length.last(); }
   V      segment            (int i) const { return _array[i+1]-_array[i];}
   Greal segment_length     (int i) const { return segment(i).length(); }

   // Editing
   
   // Prepend poly to our list
   void   prepend( L * poly ) { L pts(*poly); 
                                for (int i = 0; i < _num; i++)
                                   pts += (*this)[i];

                                clear();   // Is this necessary?
                                *this = pts;
                              }
   // Append poly to our list
   void   append ( L *poly ) { *this += *poly; }
   L*     clone_piece( int k1, int k2 ) const;
   V      normal() const {
               V  ret(V::Y);
               if (num() > 2) {
                  const int delta = num() / (num() == 3 ? 3 : 4);
                  const P a = _array[0];
                  const P b = _array[delta];
                  const P c = _array[2 * delta];
                  const P d = _array[(3 * delta) % num()];
                  ret = V((c[1]-a[1])*(d[2]-b[2]) + (c[2]-a[2])*(b[1]-d[1]),
                          (c[2]-a[2])*(d[0]-b[0]) + (c[0]-a[0])*(b[2]-d[2]),
                          (c[0]-a[0])*(d[1]-b[1]) + (c[1]-a[1])*(b[0]-d[0]));
               }
               return ret;
          }

   P      average () const { P sum; for(int i=0;i<num();i++) sum=sum%(*this)[i];
                             return sum/(num() ? num():1); }
   L      inset(Greal ins) const { 
              L cur;
              T rot(T::rotation(V::Y,1.57));
              for (int i = 0; i < num(); i++)
                 cur += inset_pt(i, rot, ins);
              return cur;
          }
   P      inset_pt(int ind, const T &rot, Greal inset) const {
               P p1((*this)[(ind-1+num())%num()]);
               P p2((*this)[ ind               ]);
               P p3((*this)[(ind+1)      %num()]);

               if (num() == 2) {
                  if (ind == 1)
                       return P(p2 - (rot * (p3-p2)).normalize() * inset);
                  else return P(p2 - (rot * (p2-p1)).normalize() * inset);
               }

               return P(p2 + ((rot * (p3-p2)).normalize() + 
                              (rot * (p2-p1)).normalize() ).normalize()*inset);
            }
};

#ifdef GLUE_NEEDS_TEMPLATES_IN_H_FILE
#include "point3d.C"
#endif
#endif

---