path.h

/*****
 * path.h
 * Andy Hammerlindl 2002/05/16
 *
 * Stores a piecewise cubic spline with known control points.
 *
 * When changing the path algorithms, also update the corresponding
 * three-dimensional algorithms in path3.cc and three.asy.
 *****/

#ifndef PATH_H
#define PATH_H

#include <cfloat>

#include "mod.h"
#include "pair.h"
#include "transform.h"
#include "bbox.h"

inline double Intcap(double t) {
  if(t <= (double) Int_MIN) return (double) Int_MIN;
  if(t >= (double) Int_MAX) return (double) Int_MAX;
  return t;
}

// The are like floor and ceil, except they return an integer;
// if the argument cannot be converted to a valid integer, they return
// Int_MAX (for positive arguments) or Int_MIN (for negative arguments).

inline Int Floor(double t) {return (Int) floor(Intcap(t));}
inline Int Ceil(double t) {return (Int) ceil(Intcap(t));}

bool simpson(double& integral, double (*)(double), double a, double b,
             double acc, double dxmax);

bool unsimpson(double integral, double (*)(double), double a, double& b,
               double acc, double& area, double dxmax, double dxmin=0);

namespace camp {

void checkEmpty(Int n);

inline Int adjustedIndex(Int i, Int n, bool cycles)
{
  checkEmpty(n);
  if(cycles)
    return imod(i,n);
  else if(i < 0)
    return 0;
  else if(i >= n)
    return n-1;
  else
    return i;
}

// Used in the storage of solved path knots.
struct solvedKnot : public gc {
  pair pre;
  pair point;
  pair post;
  bool straight;
  solvedKnot() : straight(false) {}

  friend bool operator== (const solvedKnot& p, const solvedKnot& q)
  {
    return p.pre == q.pre && p.point == q.point && p.post == q.post;
  }
};

extern const double Fuzz;
extern const double Fuzz2;
extern const double Fuzz4;
extern const double sqrtFuzz;
extern const double BigFuzz;
extern const double fuzzFactor;

class path : public gc {
  bool cycles;  // If the path is closed in a loop

  Int n; // The number of knots

  mem::vector<solvedKnot> nodes;
  mutable double cached_length; // Cache length since path is immutable.

  mutable bbox box;
  mutable bbox times; // Times where minimum and maximum extents are attained.

public:
  path()
    : cycles(false), n(0), nodes(), cached_length(-1) {}

  // Create a path of a single point
  path(pair z, bool = false)
    : cycles(false), n(1), nodes(1), cached_length(-1)
  {
    nodes[0].pre = nodes[0].point = nodes[0].post = z;
    nodes[0].straight = false;
  }

  // Creates path from a list of knots.  This will be used by camp
  // methods such as the guide solver, but should probably not be used by a
  // user of the system unless he knows what he is doing.
  path(mem::vector<solvedKnot>& nodes, Int n, bool cycles = false)
    : cycles(cycles), n(n), nodes(nodes), cached_length(-1)
  {
  }

  friend bool operator== (const path& p, const path& q)
  {
    return p.cycles == q.cycles && p.nodes == q.nodes;
  }

public:
  path(solvedKnot n1, solvedKnot n2)
    : cycles(false), n(2), nodes(2), cached_length(-1)
  {
    nodes[0] = n1;
    nodes[1] = n2;
    nodes[0].pre = nodes[0].point;
    nodes[1].post = nodes[1].point;
  }

  // Copy constructor
  path(const path& p)
    : cycles(p.cycles), n(p.n), nodes(p.nodes), cached_length(p.cached_length),
      box(p.box), times(p.times)
  {}

  path unstraighten() const
  {
    path P=path(*this);
    for(int i=0; i < n; ++i)
      P.nodes[i].straight=false;
    return P;
  }

  virtual ~path()
  {
  }

  // Getting control points
  Int size() const
  {
    return n;
  }

  bool empty() const
  {
    return n == 0;
  }

  Int length() const
  {
    return cycles ? n : n-1;
  }

  bool cyclic() const
  {
    return cycles;
  }

  mem::vector<solvedKnot>& Nodes() {
    return nodes;
  }

  bool straight(Int t) const
  {
    if (cycles) return nodes[imod(t,n)].straight;
    return (t >= 0 && t < n) ? nodes[t].straight : false;
  }

  bool piecewisestraight() const
  {
    Int L=length();
    for(Int i=0; i < L; ++i)
      if(!straight(i)) return false;
    return true;
  }

  pair point(Int t) const
  {
    return nodes[adjustedIndex(t,n,cycles)].point;
  }

  pair point(double t) const;

  pair precontrol(Int t) const
  {
    return nodes[adjustedIndex(t,n,cycles)].pre;
  }

  pair precontrol(double t) const;

  pair postcontrol(Int t) const
  {
    return nodes[adjustedIndex(t,n,cycles)].post;
  }

  pair postcontrol(double t) const;

  inline double norm(const pair& z0, const pair& c0, const pair& c1,
                     const pair& z1) const {
    return Fuzz2*camp::max((c0-z0).abs2(),
                           camp::max((c1-z0).abs2(),(z1-z0).abs2()));
  }

  pair predir(Int t, bool normalize=true) const {
    if(!cycles && t <= 0) return pair(0,0);
    pair z1=point(t);
    pair c1=precontrol(t);
    pair dir=3.0*(z1-c1);
    if(!normalize) return dir;
    pair z0=point(t-1);
    pair c0=postcontrol(t-1);
    double epsilon=norm(z0,c0,c1,z1);
    if(dir.abs2() > epsilon) return unit(dir);
    dir=2.0*c1-c0-z1;
    if(dir.abs2() > epsilon) return unit(dir);
    return unit(z1-z0+3.0*(c0-c1));
  }

  pair postdir(Int t, bool normalize=true) const {
    if(!cycles && t >= n-1) return pair(0,0);
    pair c0=postcontrol(t);
    pair z0=point(t);
    pair dir=3.0*(c0-z0);
    if(!normalize) return dir;
    pair z1=point(t+1);
    pair c1=precontrol(t+1);
    double epsilon=norm(z0,c0,c1,z1);
    if(dir.abs2() > epsilon) return unit(dir);
    dir=z0-2.0*c0+c1;
    if(dir.abs2() > epsilon) return unit(dir);
    return unit(z1-z0+3.0*(c0-c1));
  }

  pair dir(Int t, Int sign, bool normalize=true) const {
    if(sign == 0) {
      pair v=predir(t,normalize)+postdir(t,normalize);
      return normalize ? unit(v) : 0.5*v;
    }
    if(sign > 0) return postdir(t,normalize);
    return predir(t,normalize);
  }

  pair dir(double t, bool normalize=true) const {
    if(!cycles) {
      if(t <= 0) return postdir((Int) 0,normalize);
      if(t >= n-1) return predir(n-1,normalize);
    }
    Int i=Floor(t);
    t -= i;
    if(t == 0) return dir(i,0,normalize);
    pair z0=point(i);
    pair c0=postcontrol(i);
    pair c1=precontrol(i+1);
    pair z1=point(i+1);
    pair a=3.0*(z1-z0)+9.0*(c0-c1);
    pair b=6.0*(z0+c1)-12.0*c0;
    pair c=3.0*(c0-z0);
    pair dir=a*t*t+b*t+c;
    if(!normalize) return dir;
    double epsilon=norm(z0,c0,c1,z1);
    if(dir.abs2() > epsilon) return unit(dir);
    dir=2.0*a*t+b;
    if(dir.abs2() > epsilon) return unit(dir);
    return unit(a);
  }

  pair postaccel(Int t) const {
    if(!cycles && t >= n-1) return pair(0,0);
    pair z0=point(t);
    pair c0=postcontrol(t);
    pair c1=precontrol(t+1);
    return 6.0*(z0+c1)-12.0*c0;
  }

  pair preaccel(Int t) const {
    if(!cycles && t <= 0) return pair(0,0);
    pair c0=postcontrol(t-1);
    pair c1=precontrol(t);
    pair z1=point(t);
    return 6.0*(z1+c0)-12.0*c1;
  }

  pair accel(Int t, Int sign) const {
    if(sign == 0) return 0.5*(preaccel(t)+postaccel(t));
    if(sign > 0) return postaccel(t);
    return preaccel(t);
  }

  pair accel(double t) const {
    if(!cycles) {
      if(t <= 0) return postaccel((Int) 0);
      if(t >= n-1) return preaccel(n-1);
    }
    Int i=Floor(t);
    t -= i;
    if(t == 0) return 0.5*(postaccel(i)+preaccel(i));
    pair z0=point(i);
    pair c0=postcontrol(i);
    pair c1=precontrol(i+1);
    pair z1=point(i+1);
    return 6.0*t*(z1-z0+3.0*(c0-c1))+6.0*(z0+c1)-12.0*c0;
  }

  // Returns the path traced out in reverse.
  path reverse() const;

  // Generates a path that is a section of the old path, using the time
  // interval given.
  path subpath(Int start, Int end) const;
  path subpath(double start, double end) const;

  // Special case of subpath used by intersect.
  void halve(path &first, path &second) const;

  // Used by picture to determine bounding box.
  bbox bounds() const;

  pair mintimes() const {
    checkEmpty(n);
    bounds();
    return camp::pair(times.left,times.bottom);
  }

  pair maxtimes() const {
    checkEmpty(n);
    bounds();
    return camp::pair(times.right,times.top);
  }

  template<class T>
  void addpoint(bbox& box, T i) const {
    box.addnonempty(point(i),times,(double) i);
  }

  template<class T>
  void addpoint(bbox& box, T i, double min, double max) const {
    static const pair I(0,1);
    pair v=I*dir(i);
    pair z=point(i);
    box.add(z+min*v);
    box.addnonempty(z+max*v);
  }

  // Return bounding box accounting for padding perpendicular to path.
  bbox bounds(double min, double max) const;

  // Return bounding box accounting for internal pen padding (but not pencap).
  bbox internalbounds(const bbox &padding) const;

  double cubiclength(Int i, double goal=-1) const;
  double arclength () const;
  double arctime (double l) const;
  double directiontime(const pair& z) const;

  pair max() const {
    checkEmpty(n);
    return bounds().Max();
  }

  pair min() const {
    checkEmpty(n);
    return bounds().Min();
  }

  // Debugging output
  friend std::ostream& operator<< (std::ostream& out, const path& p);

// Increment count if the path has a vertical component at t.
  bool Count(Int& count, double t) const;

// Count if t is in (begin,end] and z lies to the left of point(i+t).
  void countleft(Int& count, double x, Int i, double t,
                 double begin, double end, double& mint, double& maxt) const;

// Return the winding number of the region bounded by the (cyclic) path
// relative to the point z.
  Int windingnumber(const pair& z) const;

  // Transformation
  path transformed(const transform& t) const;

};

double arcLength(const pair& z0, const pair& c0, const pair& c1,
                 const pair& z1);

extern path nullpath;
extern const unsigned maxdepth;
extern const unsigned mindepth;
extern const char *nopoints;

bool intersect(double& S, double& T, path& p, path& q, double fuzz,
               unsigned depth=maxdepth);
bool intersections(double& s, double& t, std::vector<double>& S,
                   std::vector<double>& T, path& p, path& q,
                   double fuzz, bool single, bool exact,
                   unsigned depth=maxdepth);
void intersections(std::vector<double>& S, path& g,
                   const pair& p, const pair& q, double fuzz);


// Concatenates two paths into a new one.
path concat(const path& p1, const path& p2);

// Applies a transformation to the path
path transformed(const transform& t, const path& p);

inline double quadratic(double a, double b, double c, double x)
{
  return a*x*x+b*x+c;
}

class quadraticroots {
public:
  enum {NONE=0, ONE=1, TWO=2, MANY} distinct; // Number of distinct real roots.
  unsigned roots; // Total number of real roots.
  double t1,t2;   // Real roots

  quadraticroots(double a, double b, double c);
};

class Quadraticroots {
public:
  unsigned roots; // Total number of roots.
  pair z1,z2;     // Complex roots
  Quadraticroots(pair a, pair b, pair c);
};

class cubicroots {
public:
  unsigned roots; // Total number of real roots.
  double t1,t2,t3;
  cubicroots(double a, double b, double c, double d);
};

path nurb(pair z0, pair z1, pair z2, pair z3,
          double w0, double w1, double w2, double w3, Int m);

double orient2d(const pair& a, const pair& b, const pair& c);

void roots(std::vector<double> &roots, double a, double b, double c, double d);
void roots(std::vector<double> &r, double x0, double c0, double c1, double x1,
           double x);

inline bool goodroot(double t)
{
  return 0.0 <= t && t <= 1.0;
}

extern const double third;

}

#ifndef BROKEN_COMPILER
// Delete the following line to work around problems with old broken compilers.
GC_DECLARE_PTRFREE(camp::solvedKnot);
#endif

#endif