/*****
 * runpath.in
 *
 * Runtime functions for path operations.
 *
 *****/

pen      => primPen()
pair     => primPair()
path     => primPath()
transform => primTransform()
realarray* => realArray()
realarray2* => realArray2()
patharray* => pathArray()
penarray* => penArray()

#include "path.h"
#include "arrayop.h"
#include "predicates.h"

using namespace camp;
using namespace vm;

typedef array realarray;
typedef array realarray2;
typedef array patharray;

using types::realArray;
using types::realArray2;
using types::pathArray;

Int windingnumber(array *p, camp::pair z)
{
  size_t size=checkArray(p);
  Int count=0;
  for(size_t i=0; i < size; i++)
    count += read<path *>(p,i)->windingnumber(z);
  return count;
}

// Autogenerated routines:


path :nullPath()
{
  return nullpath;
}

bool ==(path a, path b)
{
  return a == b;
}

bool !=(path a, path b)
{
  return !(a == b);
}

pair point(path p, Int t)
{
  return p.point((Int) t);
}

pair point(path p, real t)
{
  return p.point(t);
}

pair precontrol(path p, Int t)
{
  return p.precontrol((Int) t);
}

pair precontrol(path p, real t)
{
  return p.precontrol(t);
}

pair postcontrol(path p, Int t)
{
  return p.postcontrol((Int) t);
}

pair postcontrol(path p, real t)
{
  return p.postcontrol(t);
}

pair dir(path p, Int t, Int sign=0, bool normalize=true)
{
  return p.dir(t,sign,normalize);
}

pair dir(path p, real t, bool normalize=true)
{
  return p.dir(t,normalize);
}

pair accel(path p, Int t, Int sign=0)
{
  return p.accel(t,sign);
}

pair accel(path p, real t)
{
  return p.accel(t);
}

real radius(path p, real t)
{
  pair v=p.dir(t,false);
  pair a=p.accel(t);
  real d=dot(a,v);
  real v2=v.abs2();
  real a2=a.abs2();
  real denom=v2*a2-d*d;
  real r=v2*sqrt(v2);
  return denom > 0 ? r/sqrt(denom) : 0.0;
}

path reverse(path p)
{
  return p.reverse();
}

path subpath(path p, Int a, Int b)
{
  return p.subpath((Int) a, (Int) b);
}

path subpath(path p, real a, real b)
{
  return p.subpath(a,b);
}

path nurb(pair z0, pair z1, pair z2, pair z3,
          real w0, real w1, real w2, real w3, Int m)
{
  return nurb(z0,z1,z2,z3,w0,w1,w2,w3,m);
}

Int length(path p)
{
  return p.length();
}

bool cyclic(path p)
{
  return p.cyclic();
}

bool straight(path p, Int t)
{
  return p.straight(t);
}

path unstraighten(path p)
{
  return p.unstraighten();
}

bool piecewisestraight(path p)
{
  return p.piecewisestraight();
}

real arclength(path p)
{
  return p.arclength();
}

real arclength(pair z0, pair c0, pair c1, pair z1)
{
  return arcLength(z0,c0,c1,z1);
}

real arctime(path p, real L)
{
  return p.arctime(L);
}

real dirtime(path p, pair z)
{
  return p.directiontime(z);
}

realarray* intersect(path p, path q, real fuzz=-1)
{
  bool exact=fuzz <= 0.0;
  if(fuzz < 0)
    fuzz=BigFuzz*::max(::max(length(p.max()),length(p.min())),
                       ::max(length(q.max()),length(q.min())));
  std::vector<real> S,T;
  real s,t;
  if(intersections(s,t,S,T,p,q,fuzz,true,exact)) {
    array *V=new array(2);
    (*V)[0]=s;
    (*V)[1]=t;
    return V;
  }
  return new array(0);
}

realarray2* intersections(path p, path q, real fuzz=-1)
{
  bool exact=fuzz <= 0.0;
  if(fuzz < 0.0)
    fuzz=BigFuzz*::max(::max(length(p.max()),length(p.min())),
                       ::max(length(q.max()),length(q.min())));
  real s,t;
  std::vector<real> S,T;
  intersections(s,t,S,T,p,q,fuzz,false,true);
  size_t n=S.size();
  if(n == 0 && !exact) {
    if(intersections(s,t,S,T,p,q,fuzz,true,false)) {
      array *V=new array(1);
      array *Vi=new array(2);
      (*V)[0]=Vi;
      (*Vi)[0]=s;
      (*Vi)[1]=t;
      return V;
    }
  }
  array *V=new array(n);
  for(size_t i=0; i < n; ++i) {
    array *Vi=new array(2);
    (*V)[i]=Vi;
    (*Vi)[0]=S[i];
    (*Vi)[1]=T[i];
  }
  stable_sort(V->begin(),V->end(),run::compare2<real>());
  return V;
}

realarray* intersections(path p, explicit pair a, explicit pair b, real fuzz=-1)
{
  if(fuzz < 0)
    fuzz=BigFuzz*::max(::max(length(p.max()),length(p.min())),
                       ::max(length(a),length(b)));
  std::vector<real> S;
  intersections(S,p,a,b,fuzz);
  sort(S.begin(),S.end());
  size_t n=S.size();
  array *V=new array(n);
  for(size_t i=0; i < n; ++i)
    (*V)[i]=S[i];
  return V;
}

// Return the intersection point of the extensions of the line segments
// PQ and pq.
pair extension(pair P, pair Q, pair p, pair q)
{
  pair ac=P-Q;
  pair bd=q-p;
  real det=ac.getx()*bd.gety()-ac.gety()*bd.getx();
  if(det == 0) return pair(infinity,infinity);
  return P+((p.getx()-P.getx())*bd.gety()-(p.gety()-P.gety())*bd.getx())*ac/det;
}

Int size(path p)
{
  return p.size();
}

path &(path p, path q)
{
  return camp::concat(p,q);
}

pair min(explicit path p)
{
  return p.min();
}

pair max(explicit path p)
{
  return p.max();
}

Int size(patharray *p)
{
  size_t size=checkArray(p);
  Int count=0;
  for (size_t i = 0; i < size; i++)
    count += read<path *>(p,i)->size();
  return count;
}

pair min(patharray *p)
{
  size_t size=checkArray(p);

  if(size == 0)
    error(nopoints);

  path *g = p->read<path *>(0);
  pair z = g->min();
  double minx = z.getx(), miny = z.gety();

  for (size_t i = 1; i < size; ++i) {
    path *g = p->read<path *>(i);
    pair z = g->min();
    double x = z.getx(), y = z.gety();
    if (x < minx)
      minx = x;
    if (y < miny)
      miny = y;
  }

  return pair(minx, miny);
}

pair max(patharray *p)
{
  size_t size=checkArray(p);

  if(size == 0)
    error(nopoints);

  path *g = p->read<path *>(0);
  pair z = g->max();
  double maxx = z.getx(), maxy = z.gety();

  for (size_t i = 1; i < size; ++i) {
    path *g = p->read<path *>(i);
    pair z = g->max();
    double x = z.getx(), y = z.gety();
    if (x > maxx)
      maxx = x;
    if (y > maxy)
      maxy = y;
  }

  return pair(maxx, maxy);
}

pair minAfterTransform(transform t, patharray *p)
{
  size_t size=checkArray(p);

  if(size == 0)
    error(nopoints);

  path g = p->read<path *>(0)->transformed(t);
  pair z = g.min();
  double minx = z.getx(), miny = z.gety();

  for (size_t i = 1; i < size; ++i) {
    path g = p->read<path *>(i)->transformed(t);
    pair z = g.min();
    double x = z.getx(), y = z.gety();
    if (x < minx)
      minx = x;
    if (y < miny)
      miny = y;
  }

  return pair(minx, miny);
}

pair maxAfterTransform(transform t, patharray *p)
{
  size_t size=checkArray(p);

  if(size == 0)
    error(nopoints);

  path g = p->read<path *>(0)->transformed(t);
  pair z = g.max();
  double maxx = z.getx(), maxy = z.gety();

  for (size_t i = 1; i < size; ++i) {
    path g = p->read<path *>(i)->transformed(t);
    pair z = g.max();
    double x = z.getx(), y = z.gety();
    if (x > maxx)
      maxx = x;
    if (y > maxy)
      maxy = y;
  }

  return pair(maxx, maxy);
}

realarray *mintimes(path p)
{
  array *V=new array(2);
  pair z=p.mintimes();
  (*V)[0]=z.getx();
  (*V)[1]=z.gety();
  return V;
}

realarray *maxtimes(path p)
{
  array *V=new array(2);
  pair z=p.maxtimes();
  (*V)[0]=z.getx();
  (*V)[1]=z.gety();
  return V;
}

real relativedistance(real theta, real phi, real t, bool atleast)
{
  return camp::velocity(theta,phi,tension(t,atleast));
}

Int windingnumber(patharray *p, pair z)
{
  return windingnumber(p,z);
}

bool inside(explicit patharray *g, pair z, pen fillrule=CURRENTPEN)
{
  return fillrule.inside(windingnumber(g,z));
}

bool inside(path g, pair z, pen fillrule=CURRENTPEN)
{
  return fillrule.inside(g.windingnumber(z));
}

// Return a positive (negative) value if a--b--c--cycle is oriented
// counterclockwise (clockwise) or zero if all three points are colinear.
// Equivalently, return a positive (negative) value if c lies to the
// left (right) of the line through a and b or zero if c lies on this line.
// The value returned is the determinant
// |a.x a.y 1|
// |b.x b.y 1|
// |c.x c.y 1|
//
real orient(pair a, pair b, pair c)
{
  return orient2d(a,b,c);
}

// Return a positive (negative) value if d lies inside (outside)
// the circle passing through the counterclockwise-oriented points a,b,c
// or zero if d lies on this circle.
// The value returned is the determinant
// |a.x a.y a.x^2+a.y^2 1|
// |b.x b.y b.x^2+b.y^2 1|
// |c.x c.y c.x^2+c.y^2 1|
// |d.x d.y d.x^2+d.y^2 1|
real incircle(pair a, pair b, pair c, pair d)
{
  return incircle(a.getx(),a.gety(),b.getx(),b.gety(),c.getx(),c.gety(),
                  d.getx(),d.gety());
}