path2d.c

Go to the documentation of this file.

#include "path2d.h"
BOOL RectContainsPoint(RECTDBL*rc,double x,double y){
 return (x>=rc->left&&y>=rc->top&&x<rc->right&&y<rc->bottom);
}
BOOL RectIntersectsRect(RECTDBL*rc,double x,double y,double x1,double y1){
 if(IsRectEmpty(rc)||x>=x1||y>=y1)
  return FALSE;
 return (x1>rc->left&&y1>rc->top&&x<rc->right&&y<rc->bottom);
}
BOOL RectContainsRect(RECTDBL*rc,double x,double y,double x1,double y1){
 if(IsRectEmpty(rc)||x>=x1||y>=y1)
  return FALSE;
 return (x>=rc->left&&y>=rc->top&&x1<=rc->right&&y1<=rc->bottom);
}


BOOL ViewSpaceCalcTransformEx(
 MATRIX*xfmDst,
 RECTDBL*rcViewport,
 RECTDBL*rcWindow,
 int align,
 int meetOrSlice
){
 double xScale,yScale;
 double result[6];
 double vcx,vcy,wcx,wcy;
 if(!xfmDst||!rcViewport||!rcWindow)
  return FALSE;
 vcx=RectWidth(rcViewport);
 vcy=RectHeight(rcViewport);
 wcx=RectWidth(rcWindow);
 wcy=RectHeight(rcWindow);
 if(vcx<0 || vcy<0 || wcx<0 || wcy<0){
  MatrixSetIdentity(xfmDst);
  return FALSE;
 }
 if(ISZERO(wcx)||ISZERO(wcy)){
  xScale=0.0;
  yScale=0.0;
 } else {
  xScale=(double)vcx/wcx;
  yScale=(double)vcy/wcy;
 }
 result[MatrixM12]=result[MatrixM21]=0.0;
 if(align==0){
  result[MatrixM11]=xScale;
  result[MatrixM22]=yScale;
  result[MatrixDx]=(double)rcViewport->left-xScale*(double)rcWindow->left;
  result[MatrixDy]=(double)rcViewport->top-yScale*(double)rcWindow->top;
 }
 else {
  double xt,yt;
  double xw,xh,x,y;
  align-=1;
  if(meetOrSlice){ // if slice
   result[MatrixM11]=max(xScale,yScale);
  } else {
   result[MatrixM11]=min(xScale,yScale);
  }
  result[MatrixM22]=result[MatrixM11];
  xw=wcx*result[MatrixM11];
  xh=wcy*result[MatrixM11];
  switch (align % 3){ //x
  case 0:
   x=0.0;
   break;
  case 1:
   x=0.5*(vcx-xw);
   break;
  case 2:
   x=(vcx-xw);
   break;
  }
  switch (align/3){ //y
  case 0:
   y=0.0;
   break;
  case 1:
   y=0.5*(vcy-xh);
   break;
  case 2:
   y=(vcy-xh);
   break;
  }
  result[MatrixDx] =
   (double)rcViewport->left-result[MatrixM11]*(double)rcWindow->left+x;
  result[MatrixDy]=
   (double)rcViewport->top-result[MatrixM22]*(double)rcWindow->top+y;
 }
 MatrixSetValueArray(xfmDst,result);
 return TRUE;
}

inline double CubeRoot(double x){
 if( x>=0.0 )return pow( x, 1.0/3.0 );
 else return -pow( -x, 1.0/3.0 );
}

double PointDistance(double x1,double y1,double x2,double y2){
 double dx=x1-x2;
 double dy=y1-y2;
 return sqrt(dx*dx+dy*dy);
}
double SegmentDistanceSq(
  double x0, double y0, double x1, double y1,
  double ptx, double pty
){
 double a = x1-x0;
 double b = y1-y0;
 double px = ptx-x0;
 double py = pty-y0;
 double len = a*px+b*py;
 double ret;
 double c;
 if(len<=0.0){
  return px*px+py*py;
 } else {
  px=a-px;
  py=b-py;
  len=a*px+b*py;
  ret=px*px+py*py;
  if(len>0.0){
   ret-=len*len/(a*a+b*b);
  }
  return ret;
 }
}


double QuadFlatnessSq(double *c){
 return SegmentDistanceSq(c[0],c[1],c[4],c[5],c[2],c[3]);
}

double CubicFlatnessSq(double *c){
 double a=SegmentDistanceSq(c[0],c[1],c[6],c[7],c[2],c[3]);
 double b=SegmentDistanceSq(c[0],c[1],c[6],c[7],c[4],c[5]);
 return max(a,b);
}

void CubicCurveFromQuadratic(double *pSrc,double *pDst){
 double a,b,c,d,e,f,g,h;
 a=pSrc[0];
 b=pSrc[1];
 c=pSrc[0]+2*(pSrc[2]-pSrc[0])/3;
 d=pSrc[1]+2*(pSrc[3]-pSrc[1])/3;
 e=pSrc[2]+2*(pSrc[4]-pSrc[2])/3;
 f=pSrc[3]+2*(pSrc[5]-pSrc[3])/3;
 g=pSrc[4];
 h=pSrc[5];
 pDst[0]=a;
 pDst[1]=b;
 pDst[2]=c;
 pDst[3]=d;
 pDst[4]=e;
 pDst[5]=f;
 pDst[6]=g;
 pDst[7]=h;
}

double PerpendicularDistance(
  double x0, double y0, double x1, double y1,
  double ptx, double pty
){
 double a = -(y0 - y1);
 double b = (x0 - x1);
 double len= sqrt(a*a+b*b);
 double c;
 if(!ISZERO(len)){
  a/=len;
  b/=len;
 }
 c=-(a*x0+b*y0);
 return a*ptx+b*pty+c;
}

void LineSubdivide(double *src,double *left,double *right,double t){
 double x0=src[0];
 double y0=src[1];
 double x1=src[2];
 double y1=src[3];
 if(left){
  left[0]=x0;
  left[1]=y0;
 }
 if(right){
  right[2]=x1;
  right[3]=y1;
 }
 x1=x0+(x1-x0)*t;
 y1=y0+(y1-y0)*t;
 if(left){
  left[2]=x1;
  left[3]=y1;
 }
 if(right){
  right[0]=x1;
  right[1]=y1;
 }
}

void QuadSubdivide(double *src,double *left,double *right,double t){
 double x0=src[0];
 double y0=src[1];
 double cx=src[2];
 double cy=src[3];
 double x1=src[4];
 double y1=src[5];
 if(left){
  left[0]=x0;
  left[1]=y0;
 }
 if(right){
  right[4]=x1;
  right[5]=y1;
 }
 x1=cx+(x1-cx)*t;
 y1=cy+(y1-cy)*t;
 x0=x0+(cx-x0)*t;
 y0=y0+(cy-y0)*t;
 cx=x0+(x1-x0)*t;
 cy=y0+(y1-y0)*t;
 if(left){
  left[2]=x0;
  left[3]=y0;
  left[4]=cx;
  left[5]=cy;
 }
 if(right){
  right[0]=cx;
  right[1]=cy;
  right[2]=x1;
  right[3]=y1;
 }
}
void CubicSubdivide(double *src,double *left,double *right,double t){
 double x0=src[0];
 double y0=src[1];
 double cx0=src[2];
 double cy0=src[3];
 double cx1=src[4];
 double cy1=src[5];
 double x1=src[6];
 double y1=src[7];
 double cx,cy;
 if(left){
  left[0]=x0;
  left[1]=y0;
 }
 if(right){
  right[6]=x1;
  right[7]=y1;
 }
 x1=cx1+(x1-cx1)*t;
 y1=cy1+(y1-cy1)*t;
 x0=x0+(cx0-x0)*t;
 y0=y0+(cy0-y0)*t;
 cx0=cx0+(cx1-cx0)*t;
 cy0=cy0+(cy1-cy0)*t;
 cx1=cx0+(x1-cx0)*t;
 cy1=cy0+(y1-cy0)*t;
 cx0=x0+(cx0-x0)*t;
 cy0=y0+(cy0-y0)*t;
 cx=cx0+(cx1-cx0)*t;
 cy=cy0+(cy1-cy0)*t;
 if(left){
  left[2]=x0;
  left[3]=y0;
  left[4]=cx0;
  left[5]=cy0;
  left[6]=cx;
  left[7]=cy;
 }
 if(right){
  right[0]=cx;
  right[1]=cy;
  right[2]=cx1;
  right[3]=cy1;
  right[4]=x1;
  right[5]=y1;
 }
}

double *PathSegmentEndPoint(PATHSEGMENT*ps){
 if(!ps)
  return NULL;
 switch (ps->type){
 default:
  return NULL;
 case PATH_MOVETO:
 case PATH_LINETO:
  return &ps->coords[0];
 case PATH_QUADTO:
  return &ps->coords[2];
 case PATH_CUBICTO:
  return &ps->coords[4];
 }
}
BOOL PathMoveTo(ITEMARRAY*ia,double x,double y){
 PATHSEGMENT *ps=NULL;
 if(!ISARRAY(ia,PATHSEGMENT))return FALSE;
 ps=ItemArrayAllocOnePtr(ia);
 if(!ps)return FALSE;
 ps->type=PATH_MOVETO;
 ps->coords[0]=x;
 ps->coords[1]=y;
 return TRUE;
}
BOOL PathAddPoint(ITEMARRAY *ia, double x, double y){
 double endpoint[2];
 BOOL haveendpoint;
 if(!ISARRAY(ia,PATHSEGMENT))return FALSE;
 haveendpoint=PathGetEndPoint(ia,endpoint);
 if(!haveendpoint||endpoint[0]!=x||endpoint[1]!=y){
  return PathMoveTo(ia,x,y);
 }
 return TRUE;
}
BOOL PathLineTo(ITEMARRAY*ia,double x,double y){
 PATHSEGMENT *ps;
 if(!ISARRAY(ia,PATHSEGMENT)||ia->length==0)return 0;
 ps=ItemArrayAllocOnePtr(ia);
 if(!ps)return 0;
 ps->type=PATH_LINETO;
 ps->coords[0]=x;
 ps->coords[1]=y;
 return TRUE;
}
BOOL PathQuadTo(ITEMARRAY*ia,double x,double y,double x1,double y1){
 PATHSEGMENT *ps;
 if(!ISARRAY(ia,PATHSEGMENT)||ia->length==0)return FALSE;
 ps=ItemArrayAllocOnePtr(ia);
 if(!ps)return FALSE;
 ps->type=PATH_QUADTO;
 ps->coords[0]=x;
 ps->coords[1]=y;
 ps->coords[2]=x1;
 ps->coords[3]=y1;
 return TRUE;
}
BOOL PathCubicTo(ITEMARRAY*ia,
 double x,double y,double x1,double y1,double x2,double y2){
 PATHSEGMENT *ps;
 if(!ISARRAY(ia,PATHSEGMENT)||ia->length==0)return 0;
 ps=ItemArrayAllocOnePtr(ia);
 if(!ps)return FALSE;
 ps->type=PATH_CUBICTO;
 ps->coords[0]=x;
 ps->coords[1]=y;
 ps->coords[2]=x1;
 ps->coords[3]=y1;
 ps->coords[4]=x2;
 ps->coords[5]=y2;
 return TRUE;
}

BOOL PathSmoothQuadTo(ITEMARRAY *ia, double x, double y){
 PATHSEGMENT *ps;
 if(!ISARRAY(ia,PATHSEGMENT)||ia->length==0)return 0; 
 ps=ItemArrayPtr(ia,ia->length-1);
 if(ps->type==PATH_QUADTO){
  double newx=ps->coords[2]*2 - ps->coords[0];
  double newy=ps->coords[3]*2 - ps->coords[1];
  return PathQuadTo(ia,newx,newy,x,y);
 } else {
  double endpt[2];
  if(!PathGetEndPoint(ia,endpt))
   return FALSE;
  return PathQuadTo(ia,endpt[0],endpt[1],x,y);
 }
}


BOOL PathSmoothCubicTo(ITEMARRAY *ia, double x1, double y1, double x2, double y2){
 PATHSEGMENT *ps;
 if(!ISARRAY(ia,PATHSEGMENT)||ia->length==0)return 0; 
 ps=ItemArrayPtr(ia,ia->length-1);
 if(ps->type==PATH_CUBICTO){
  double newx=ps->coords[4]*2 - ps->coords[2];
  double newy=ps->coords[5]*2 - ps->coords[3];
  return PathCubicTo(ia,newx,newy,x1,y1,x2,y2);
 } else {
  double endpt[2];
  if(!PathGetEndPoint(ia,endpt))
   return FALSE;
  return PathCubicTo(ia,endpt[0],endpt[1],x1,y1,x2,y2);
 }
}


BOOL PathClose(ITEMARRAY*ia){
 PATHSEGMENT *ps=NULL;
 if(!ISARRAY(ia,PATHSEGMENT)||ia->length==0)return 0;
 ps=ItemArrayPtr(ia,ia->length-1);
 if(ps->type==PATH_CLOSE)
  return TRUE;
 ps=ItemArrayAllocOnePtr(ia);
 if(!ps)return FALSE;
 ps->type=PATH_CLOSE;
 return TRUE;
}
BOOL PathAddSegment(ITEMARRAY*ia,PATHSEGMENT*ps){
 if(!ISARRAY(ia,PATHSEGMENT)||!ps)
  return FALSE;
 switch (ps->type){
 case PATH_MOVETO:
  return PathMoveTo(ia,ps->coords[0],ps->coords[1]);
 case PATH_CLOSE:
  return PathClose(ia);
 case PATH_LINETO:
  return PathLineTo(ia,ps->coords[0],ps->coords[1]);
 case PATH_QUADTO:
  return PathQuadTo(ia,ps->coords[0],ps->coords[1],
     ps->coords[2],ps->coords[3]);
 case PATH_CUBICTO:
  return PathCubicTo(ia,ps->coords[0],ps->coords[1],
     ps->coords[2],ps->coords[3],ps->coords[4],ps->coords[5]);
 default:
  return FALSE;
 }
}
BOOL PathAppendPath(
 ITEMARRAY*ia,
 PATHITERATOR *pcb,
 LPVOID data,
 MATRIX*xfm,
 BOOL connect
){
 PATHSEGMENT ps;
 LPVOID lParam;
 int i;
 PATHSEGMENT *pps;
 double *endpt;
 if(!pcb||!ISARRAY(ia,PATHSEGMENT))
  return FALSE;
 lParam=pcb->Init(data);
 if(!lParam)
  return FALSE;
 while(pcb->Next(lParam,&ps)>0){
  PathSegmentTransform(&ps,xfm);
  switch (ps.type){
  case PATH_MOVETO:
   if(!connect||ia->length==0){
    PathMoveTo(ia,ps.coords[0],ps.coords[1]);
    break;
   }
   pps=ItemArrayPtr(ia,ia->length-1);
   endpt=PathSegmentEndPoint(pps);
   if(endpt&&endpt[0]==ps.coords[0]&&endpt[1]==ps.coords[1]){
   } else {
    PathLineTo(ia,ps.coords[0],ps.coords[1]);
   }
   break;
  case PATH_LINETO:
   PathLineTo(ia,ps.coords[0],ps.coords[1]);
   break;
  case PATH_QUADTO:
   PathQuadTo(ia,ps.coords[0],ps.coords[1],ps.coords[2],ps.coords[3]);
   break;
  case PATH_CUBICTO:
   PathCubicTo(ia,ps.coords[0],ps.coords[1],
    ps.coords[2],ps.coords[3],ps.coords[4],ps.coords[5]);
   break;
  case PATH_CLOSE:
   PathClose(ia);
   break;
  }
  connect=FALSE;
 }
 pcb->Free(lParam);
 return TRUE;
}


void DebugOutFullPathSegment(FULLPATHSEGMENT*ps,LPCTSTR add){
#ifdef DEBUG
 if(add){
  DebugOutF(add);
  DebugOutF(" ");
 }
 switch (ps->type){
 case PATH_MOVETO:
 case PATH_LINETO:
 case PATH_CLOSE:
  DebugOut("[%d] %f %f %f %f",ps->type,ps->coords[0],ps->coords[1],ps->coords[2],ps->coords[3]);
  break;
 case PATH_QUADTO:
  DebugOut("[%d] %f %f %f %f %f %f",ps->type,ps->coords[0],ps->coords[1],
   ps->coords[2],ps->coords[3],ps->coords[4],ps->coords[5]);
  break;
 case PATH_CUBICTO:
  DebugOut("[%d] %f %f %f %f %f %f %f %f",ps->type,ps->coords[0],ps->coords[1],
   ps->coords[2],ps->coords[3],ps->coords[4],ps->coords[5],ps->coords[6],ps->coords[7]);
  break;
 }
#endif
}

void DebugOutPathSegment(PATHSEGMENT*ps,LPCTSTR add){
#ifdef DEBUG
 if(add){
  DebugOutF(add);
  DebugOutF(" ");
 }
 switch (ps->type){
 case PATH_MOVETO:
 case PATH_LINETO:
  DebugOut("[%d] %f %f",ps->type,ps->coords[0],ps->coords[1]);
  break;
 case PATH_QUADTO:
  DebugOut("[%d] %f %f %f %f",ps->type,ps->coords[0],ps->coords[1],
   ps->coords[2],ps->coords[3]);
  break;
 case PATH_CUBICTO:
  DebugOut("[%d] %f %f %f %f %f %f",ps->type,ps->coords[0],ps->coords[1],
   ps->coords[2],ps->coords[3],ps->coords[4],ps->coords[5]);
  break;
 case PATH_CLOSE:
  DebugOut("[%d]",ps->type);
  break;
 }
#endif
}

void DebugOutPathSegmentAsSvg(PATHSEGMENT *ps){
#ifdef DEBUG
 switch(ps->type){
  case PATH_CLOSE:
   DebugOutF("Z ");
   break;
  case PATH_MOVETO:
   DebugOutF("M%f,%f ",ps->coords[0],ps->coords[1]);
   break;
  case PATH_LINETO:
   DebugOutF("L%f,%f ",ps->coords[0],ps->coords[1]);
   break;
  case PATH_QUADTO:
   DebugOutF("Q%f,%f,%f,%f ",ps->coords[0],ps->coords[1],
      ps->coords[2],ps->coords[3]);
   break;
  case PATH_CUBICTO:
   DebugOutF("C%f,%f,%f,%f,%f,%f ",ps->coords[0],ps->coords[1],
      ps->coords[2],ps->coords[3],
      ps->coords[4],ps->coords[5]);
   break;
 }
#endif
}


BOOL PathGetEndPoint(ITEMARRAY *ia,double *endpt){
 LONG realtype;
 PATHSEGMENT *ps;
 LONG i;
 if(!endpt||!ISARRAY(ia,PATHSEGMENT)||ia->length==0){
  return FALSE;
 }
 i=ia->length-1;
 ps=(PATHSEGMENT*)ItemArrayPtr(ia,i);
 if(ps->type==PATH_CLOSE){
  int j;
  for(j=i;j>=0;j--){
   ps=(PATHSEGMENT*)ItemArrayPtr(ia,j);
   if(ps->type==PATH_MOVETO)
    break;
  }
  if(j<0){
   return FALSE;
  }
 }
 switch(ps->type){
  case PATH_MOVETO:
  case PATH_LINETO:
   endpt[0]=ps->coords[0];
   endpt[1]=ps->coords[1];
   break;
  case PATH_QUADTO:
   endpt[0]=ps->coords[2];
   endpt[1]=ps->coords[3];
   break;
  case PATH_CUBICTO:
   endpt[0]=ps->coords[4];
   endpt[1]=ps->coords[5];
   break;
  default:
   return FALSE;
 }
 return TRUE;
}


void PathCardinalSplineTo(
  ITEMARRAY *path,
  double *coords, 
  LONG numpoints, 
  double tension
){
 PTDBL *pts=(PTDBL*)coords;
 PTDBL point0;
 LONG i;
 double x0,x1,y0,y1;
 if(numpoints<1)return;
 if(!PathGetEndPoint(path,(double*)&point0))
  return;
 tension/=3;
 for(i=0;i<numpoints;i++){
  LONG iprev=i-1;
  LONG ithispt=i;
  LONG inext1=i+1;
  LONG inext2=i+2;
  PTDBL *prev,*thispt,*next1,*next2;
  prev=(iprev==0) ? &point0 : &pts[iprev-1];
  thispt=(ithispt==0) ? &point0 : &pts[ithispt-1];
  next1=(inext1==0) ? &point0 : &pts[inext1-1];
  next2=(inext2==0) ? &point0 : &pts[inext2-1];
  if(iprev<0){
   // At start of curve
   ASSERT(ithispt==0);
   x0=thispt->x+(next1->x-thispt->x)*tension;
   y0=thispt->y+(next1->y-thispt->y)*tension;
  } else {
   x0=thispt->x+(next1->x-prev->x)*tension;
   y0=thispt->y+(next1->y-prev->y)*tension;
  }
  if(inext2>=numpoints+1){
   // At end of curve
   x1=next1->x+(thispt->x-next1->x)*tension;
   y1=next1->y+(thispt->y-next1->y)*tension;  
  } else {
   x1=next1->x-(next2->x-thispt->x)*tension;
   y1=next1->y-(next2->y-thispt->y)*tension;
  }
  PathCubicTo(path,x0,y0,x1,y1,next1->x,next1->y);
 }
}

void PathAddCardinalSpline(
  ITEMARRAY *path, 
  double *coords, 
  LONG numpoints, 
  double tension,
  BOOL closefigure
){
 PTDBL *pts=(PTDBL*)coords;
 LONG i;
 double x0,x1,y0,y1;
 if(closefigure && numpoints<3)return;
 if(numpoints<2)return;
 tension/=3;
 PathMoveTo(path,pts[0].x,pts[0].y);
 for(i=0;i<numpoints;i++){
  LONG prev=i-1;
  LONG thispt=i;
  LONG next1=i+1;
  LONG next2=i+2;
  if(thispt==0 && closefigure)
   prev=numpoints-1;
  if(thispt==numpoints-2 && closefigure)
   next2=0;
  if(thispt==numpoints-1){
   if(!closefigure)break;
   next1=0;
   next2=1;
  }
  if(prev<0){
   // At start of curve
   ASSERT(thispt==0);
   x0=pts[thispt].x+(pts[next1].x-pts[thispt].x)*tension;
   y0=pts[thispt].y+(pts[next1].y-pts[thispt].y)*tension;
  } else {
   x0=pts[thispt].x+(pts[next1].x-pts[prev].x)*tension;
   y0=pts[thispt].y+(pts[next1].y-pts[prev].y)*tension;
  }
  if(next2>=numpoints){
   // At end of curve
   x1=pts[next1].x+(pts[thispt].x-pts[next1].x)*tension;
   y1=pts[next1].y+(pts[thispt].y-pts[next1].y)*tension;  
  } else {
   x1=pts[next1].x-(pts[next2].x-pts[thispt].x)*tension;
   y1=pts[next1].y-(pts[next2].y-pts[thispt].y)*tension;
  }
  PathCubicTo(path,x0,y0,x1,y1,pts[next1].x,pts[next1].y);
 }
 if(closefigure){
  PathClose(path);
 }
}

BOOL PathTransform(ITEMARRAY *path,MATRIX*xfm){
 DWORD i;
 if(!ISARRAY(path,PATHSEGMENT)||!xfm)return FALSE;
 for(i=0;i<path->length;i++){
  PATHSEGMENT *seg=ItemArrayPtr(path,i);
  PathSegmentTransform(seg,xfm);
 }
 return TRUE;
}

BOOL PathSegmentTransform(PATHSEGMENT*ps,MATRIX*xfm){
 if(!ps)
  return FALSE;
 if(!xfm)
  return TRUE;
 switch (ps->type){
 case PATH_MOVETO:
 case PATH_LINETO:
  MatrixTransformPoints(xfm,ps->coords,1);
  return TRUE;
 case PATH_QUADTO:
  MatrixTransformPoints(xfm,ps->coords,2);
  return TRUE;
 case PATH_CUBICTO:
  MatrixTransformPoints(xfm,ps->coords,3);
  return TRUE;
 case PATH_CLOSE:
  return TRUE;
 default:
  return FALSE;
 }
}

//     Rectangles
typedef struct{
 int stage;
 double coords[8];
} SEGITERATOR;

static int SegIterator_GetWinding(LPVOID h){
 return WIND_NONZERO; 
}
static LPVOID SegIterator_Init(double *data, DWORD numpoints){
 double *rect=(double*)data;
 SEGITERATOR *iterator;
 if(!data)return NULL;
 iterator=malloc(sizeof(SEGITERATOR));
 if(iterator){
  DWORD i;
  iterator->stage=0;
  // Copy left, top, right, bottom of source rect
  for(i=0;i<numpoints;i++){
   iterator->coords[i]=data[i];
  }
 }
 return iterator;
}
static void SegIterator_Free(LPVOID h){
 free(h);
}
static LPVOID LineIterator_Init(double *data){
 return SegIterator_Init(data,4);
}
static int LineIterator_Next(LPVOID h, PATHSEGMENT *ps){
 SEGITERATOR *rc=(SEGITERATOR *)h;
 if(!h||!ps)return -1;
 switch (rc->stage){
  case 0:
   ps->type=PATH_MOVETO;
   ps->coords[0]=rc->coords[0];
   ps->coords[1]=rc->coords[1];
   break;
  case 1:
   ps->type=PATH_LINETO;
   ps->coords[0]=rc->coords[2];
   ps->coords[1]=rc->coords[3];
   break;
  default:
   return 0;
 }
 rc->stage++;
 return 1;
}
static LPVOID QuadIterator_Init(double *data){
 return SegIterator_Init(data,6);
}
static int QuadIterator_Next(LPVOID h, PATHSEGMENT *ps){
 SEGITERATOR *rc=(SEGITERATOR *)h;
 if(!h||!ps)return -1;
 switch (rc->stage){
  case 0:
   ps->type=PATH_MOVETO;
   ps->coords[0]=rc->coords[0];
   ps->coords[1]=rc->coords[1];
   break;
  case 1:
   ps->type=PATH_QUADTO;
   ps->coords[0]=rc->coords[2];
   ps->coords[1]=rc->coords[3];
   ps->coords[2]=rc->coords[4];
   ps->coords[3]=rc->coords[5];
   break;
  default:
   return 0;
 }
 rc->stage++;
 return 1;
}
static LPVOID CubicIterator_Init(double *data){
 return SegIterator_Init(data,8);
}
static int CubicIterator_Next(LPVOID h, PATHSEGMENT *ps){
 SEGITERATOR *rc=(SEGITERATOR *)h;
 if(!h||!ps)return -1;
 switch (rc->stage){
  case 0:
   ps->type=PATH_MOVETO;
   ps->coords[0]=rc->coords[0];
   ps->coords[1]=rc->coords[1];
   break;
  case 1:
   ps->type=PATH_CUBICTO;
   ps->coords[0]=rc->coords[2];
   ps->coords[1]=rc->coords[3];
   ps->coords[2]=rc->coords[4];
   ps->coords[3]=rc->coords[5];
   ps->coords[2]=rc->coords[6];
   ps->coords[3]=rc->coords[7];
   break;
  default:
   return 0;
 }
 rc->stage++;
 return 1;
}
typedef struct{
 int stage;
 double rect[4];
} RECTITERATOR;
static int RectPath_GetWinding(LPVOID h){
 return WIND_NONZERO; 
}
static LPVOID RectPath_Init(LPVOID data){
 double *rect=(double*)data;
 RECTITERATOR *rectiter;
 if(!data)return NULL;
 rectiter=malloc(sizeof(RECTITERATOR));
 if(rectiter){
  rectiter->stage=0;
  // Copy left, top, right, bottom of source rect
  rectiter->rect[0]=rect[0];
  rectiter->rect[1]=rect[1];
  rectiter->rect[2]=rect[2];
  rectiter->rect[3]=rect[3];
 }
 return rectiter;
}
static void RectPath_Free(LPVOID h){
 free(h);
}
static int RectPath_Next(LPVOID h, PATHSEGMENT *ps){
 RECTITERATOR *rc=(RECTITERATOR *)h;
 if(!h||!ps)return -1;
 switch (rc->stage){
  case 0:
   ps->type=PATH_MOVETO;
   ps->coords[0]=rc->rect[0];
   ps->coords[1]=rc->rect[1];
   break;
  case 1:
   ps->type=PATH_LINETO;
   ps->coords[0]=rc->rect[2];
   ps->coords[1]=rc->rect[1];
   break;
  case 2:
   ps->type=PATH_LINETO;
   ps->coords[0]=rc->rect[2];
   ps->coords[1]=rc->rect[3];
   break;
  case 3:
   ps->type=PATH_LINETO;
   ps->coords[0]=rc->rect[0];
   ps->coords[1]=rc->rect[3];
   break;
  case 4:
   ps->type=PATH_CLOSE;
   break;
  default:
   return 0;
 }
 rc->stage++;
 return 1;
}

typedef struct{
 int stage;
 RECTDBL rect;
} ELLIPSEITERATOR;

#define PCV 0.77614237491539666
#define NCV 0.22385762508460333
static double ctrlpts[4][6]={
 {1.0,PCV,PCV,1.0,0.5,1.0},
 {NCV,1.0,0.0,PCV,0.0,0.5},
 {0.0,NCV,NCV,0.0,0.5,0.0},
 {PCV,0.0,1.0,NCV,1.0,0.5}
};
static int EllipsePath_GetWinding(LPVOID h){
 return WIND_NONZERO; 
}
static LPVOID EllipsePath_Init(LPVOID data){
 double *rect=(double*)data;
 ELLIPSEITERATOR *rectiter;
 if(!data)return NULL;
 rectiter=malloc(sizeof(ELLIPSEITERATOR));
 if(rectiter){
  rectiter->stage=0;
  // Copy left, top, right, bottom of source rect
  rectiter->rect.left=rect[0];
  rectiter->rect.top=rect[1];
  rectiter->rect.right=rect[2];
  rectiter->rect.bottom=rect[3];
 }
 return rectiter;
}
static void EllipsePath_Free(LPVOID h){
 free(h);
}
static int EllipsePath_Next(LPVOID handle, PATHSEGMENT *ps){
 ELLIPSEITERATOR *rc=(ELLIPSEITERATOR *)handle;
 if(!handle||!ps)return -1;
   double *ctrls;
   double w,h;
   if(rc->stage>=6)
    return 0;
   if(rc->stage==5){
    ps->type=PATH_CLOSE;
    rc->stage++;
    return 1;
   }
   w=RectWidth(&rc->rect);
   h=RectHeight(&rc->rect);
   if(rc->stage==0){
    ps->type=PATH_MOVETO;
    ctrls=ctrlpts[3];
    ps->coords[0]=rc->rect.left+ctrls[4]*w;
    ps->coords[1]=rc->rect.top+ctrls[5]*h;
    rc->stage++;
    return 1;
   }
   ctrls=ctrlpts[rc->stage-1];
   ps->type=PATH_CUBICTO;
   ps->coords[0]=rc->rect.left+ctrls[0]*w;
   ps->coords[1]=rc->rect.top+ctrls[1]*h;
   ps->coords[2]=rc->rect.left+ctrls[2]*w;
   ps->coords[3]=rc->rect.top+ctrls[3]*h;
   ps->coords[4]=rc->rect.left+ctrls[4]*w;
   ps->coords[5]=rc->rect.top+ctrls[5]*h;
   rc->stage++;
   return 1;
}

#define NCV 0.22385762508460333
static double rrctrlpts[]={
 0.0,0.0,0.0,0.5,
 0.0,0.0,1.0,-0.5,
 0.0,0.0,1.0,-NCV,
 0.0,NCV,1.0,0.0,
 0.0,0.5,1.0,0.0,
 1.0,-0.5,1.0,0.0,
 1.0,-NCV,1.0,0.0,
 1.0,0.0,1.0,-NCV,
 1.0,0.0,1.0,-0.5,
 1.0,0.0,0.0,0.5,
 1.0,0.0,0.0,NCV,
 1.0,-NCV,0.0,0.0,
 1.0,-0.5,0.0,0.0,
 0.0,0.5,0.0,0.0,
 0.0,NCV,0.0,0.0,
 0.0,0.0,0.0,NCV,
 0.0,0.0,0.0,0.5,
};
static int rrindices[]={
 0,4,8,20,24,36,40,52,56,68,68
};
static int rrtypes[]={
 PATH_MOVETO,
 PATH_LINETO,PATH_CUBICTO,
 PATH_LINETO,PATH_CUBICTO,
 PATH_LINETO,PATH_CUBICTO,
 PATH_LINETO,PATH_CUBICTO,
 PATH_CLOSE
};
typedef struct{
 int stage;
 ROUNDRECT rect;
} RRITERATOR;

static int RoundRectPath_GetWinding(LPVOID h){
 return WIND_NONZERO; 
}
static LPVOID RoundRectPath_Init(LPVOID data){
 ROUNDRECT *rect=(ROUNDRECT *)data;
 RRITERATOR *rectiter;
 if(!data)return NULL;
 rectiter=malloc(sizeof(RRITERATOR));
 if(rectiter){
  rectiter->stage=0;
  // Copy left, top, right, bottom of source rect
  rectiter->rect=*rect;
 }
 return rectiter;
}
static void RoundRectPath_Free(LPVOID h){
 free(h);
}
static int RoundRectPath_Next(LPVOID handle, PATHSEGMENT *ps){
 RRITERATOR *rc=(RRITERATOR *)handle;
 ROUNDRECT *rrc;
 double *ctrls;
 double w,h,aw,ah;
 int nc=0,i;
 if(!handle||!ps)return -1;
 rrc=&rc->rect;
 if(rc->stage>=DIMOF(rrtypes))
  return 0;
 w=RectWidth(rrc);
 h=RectHeight(rrc);
 aw=min(w,rrc->cxArc);
 ah=min(h,rrc->cyArc);
 ps->type=rrtypes[rc->stage];
 for(i=rrindices[rc->stage];i<rrindices[rc->stage+1];i+=4){
  ps->coords[nc++]=(rrc->left+rrctrlpts[i+1]*aw+rrctrlpts[i]*w);
  ps->coords[nc++]=(rrc->top+rrctrlpts[i+3]*ah+rrctrlpts[i+2]*h);
 }
 rc->stage++;
 return 1;
}

PATHITERATOR RoundRectPath={
 RoundRectPath_Init,
 RoundRectPath_Next,
 RoundRectPath_Free,
 RoundRectPath_GetWinding
};
PATHITERATOR EllipsePath={
 EllipsePath_Init,
 EllipsePath_Next,
 EllipsePath_Free,
 EllipsePath_GetWinding
};
PATHITERATOR RectPath={
 RectPath_Init,
 RectPath_Next,
 RectPath_Free,
 RectPath_GetWinding
};
PATHITERATOR LinePath={
 LineIterator_Init,
 LineIterator_Next,
 SegIterator_Free,
 SegIterator_GetWinding
};
PATHITERATOR QuadPath={
 QuadIterator_Init,
 QuadIterator_Next,
 SegIterator_Free,
 SegIterator_GetWinding
};
PATHITERATOR CubicPath={
 CubicIterator_Init,
 CubicIterator_Next,
 SegIterator_Free,
 SegIterator_GetWinding
};

//     Arcs

static BOOL CubicCurveFromArc(
  double x0, double y0, double rx, double ry,
  double startangle, double sweepangle, double *P
){
  startangle=fmod(startangle,M_PITIMESTWO);
  double end=startangle+sweepangle;
  end=fmod(end,M_PITIMESTWO);
  double bcp = (4.0/3 * (1 - cos(sweepangle/2)) / sin(sweepangle/2));
  double sin_start = sin(startangle);
  double sin_end =  sin(end);
  double cos_start = cos(startangle);
  double cos_end =  cos(end);
  P[0] = x0 + rx * cos_start;
  P[1] = y0 + ry * sin_start;
  P[2] = x0 + rx* (cos_start - bcp * sin_start);
  P[3] = y0 + ry* (sin_start + bcp * cos_start);
  P[4] = x0 + rx* (cos_end + bcp * sin_end);
  P[5] = y0 + ry* (sin_end - bcp * cos_end);
  P[6] = x0 + rx* cos_end;
  P[7] = y0 + ry* sin_end;
  return TRUE;
}

// Based on code in newsgroup post "Elliptical Arc" 
// in microsoft.public.dotnet.framework.drawing
// http://www.dotnet247.com/247reference/msgs/55/276471.aspx
static double FixAngle(double radiusX, double radiusY, double angle){
 double myAngle = angle;
 double diff = 0;
 while(myAngle < 0){
  myAngle += M_PITIMESTWO;
  diff -= M_PITIMESTWO;
 }
 while(myAngle >= M_PITIMESTWO){
  myAngle -= M_PITIMESTWO;
  diff += M_PITIMESTWO;
 }
 double temp = myAngle;
 myAngle = atan((radiusX / radiusY) * tan(myAngle));
 if(myAngle < 0){
   myAngle += M_HALFPI;
 }
 if(temp > M_PI*1.5){
  myAngle += M_PI*1.5;
 } else if(temp > M_PI){
  myAngle += M_PI;
 } else if(temp > M_HALFPI){
  myAngle += M_HALFPI;
 }
 myAngle += diff;
 return myAngle;
}

static int CubicCurvesFromArc(double x,
 double y,double rx,double ry,double start,double sweep,double* vertices){
 int ret=2;
 double total_sweep=0.0;
 double local_sweep=0.0;
 BOOL done=FALSE;
 if(FALSE){
  // for compatibility with GDI Plus
  double end;
  end=start+sweep;
  start=FixAngle(rx,ry,start);
  end=FixAngle(rx,ry,end);
  sweep=end-start;
 }
 if(sweep>=M_PITIMESTWO)
  sweep=M_PITIMESTWO;
 if(sweep<=-M_PITIMESTWO)
  sweep=-M_PITIMESTWO;
 do {
  if(sweep<0.0){
   local_sweep=-M_HALFPI;
   total_sweep-=M_HALFPI;
   if(total_sweep<=sweep){
    local_sweep=sweep-(total_sweep+M_HALFPI);
    done=TRUE;
   }
  }
  else {
   local_sweep=M_HALFPI;
   total_sweep+=M_HALFPI;
   if(total_sweep>=sweep){
    local_sweep=sweep-(total_sweep-M_HALFPI);
    done=TRUE;
   }
  }
  CubicCurveFromArc(x,y,rx,ry,start,local_sweep,vertices+ret-2);
  ret+=6;
  start+=local_sweep;
 } while(!done&&ret<26);
 return ret;
}

typedef enum {
 ArcCheckEndPoint=1,
 ArcDrawLine=2,
 ArcMoveTo=4,
 ArcHaveEndCoords=8
} ArcInternalFlags;

static BOOL PathAddArcInternal(
 ITEMARRAY *ia,
 ArcInternalFlags mode,
 double centerX,
 double centerY,
 double radiusX,
 double radiusY,
 double start, //in radians
 double sweep, //in radians
 double phi,   //in radians
 double x0,
 double y0,
 double x,
 double y
){
 double vertices[26];
 double endpoint[2];
 int i,numvertices;
 MATRIX xfm;
 if(!ISARRAY(ia,PATHSEGMENT))
  return FALSE;
 radiusX=abs(radiusX);
 radiusY=abs(radiusY);
 if(radiusX==0||radiusY==0){
  return TRUE;
 }
 numvertices=CubicCurvesFromArc(centerX,centerY,radiusX,radiusY,start,sweep,vertices);
 if(phi!=0){
  MatrixSetIdentity(&xfm);
  MatrixTranslate(&xfm,-centerX,-centerY,MatrixAppend);
  MatrixRotate(&xfm,phi*RADTODEG,MatrixAppend);
  MatrixTranslate(&xfm,centerX,centerY,MatrixAppend);
  MatrixTransformPoints(&xfm,vertices,numvertices>>1);
 }
 if(mode&ArcHaveEndCoords){
  vertices[0]=x0;
  vertices[1]=y0;
  vertices[numvertices-2]=x;
  vertices[numvertices-1]=y;  
 }
 if(mode&ArcCheckEndPoint){
  if(!PathGetEndPoint(ia,endpoint)){
   return FALSE;
  } 
  if(endpoint[0]!=vertices[0]||
     endpoint[1]!=vertices[1]){
   if(mode&ArcDrawLine){
    if(!PathLineTo(ia,vertices[0],vertices[1]))
     return FALSE;
   }
  }
 }
 if(mode&ArcMoveTo){
  // Move to
  if(!PathMoveTo(ia,vertices[0],vertices[1]))
   return FALSE;
 }
 if(sweep==0){
  return TRUE;
 }
 for(i=2;i<numvertices;i+=6){
  if(!PathCubicTo(ia,
    vertices[i],vertices[i+1],
    vertices[i+2],vertices[i+3],
    vertices[i+4],vertices[i+5]
  )){
   return FALSE;
  }
 }
 return TRUE;
}

static BOOL PathAddSvgArcInternal(
 ITEMARRAY *ia,
 BOOL arcto,
 double x0,   // Start point
 double y0,
 double x,    // End point
 double y,
 double rx,   // Radii
 double ry,   
 double phi,  // Ellipse rotation, in radians
 BOOL large_arc,
 BOOL sweep
){
// In SVG Essentials:
// http://commons.oreilly.com/wiki/index.php/SVG_Essentials/Paths
    double cx, cy, theta, delta;
    double dx2, dy2, x1, y1,
        rx_sq, ry_sq,
        x1_sq, y1_sq,
        sign, sq, coef,
        cx1, cy1, sx2, sy2,
        p, n,
        ux, uy, vx, vy, tmp;
    double radius_check;
    double cosPhiR, sinPhiR;
    if(arcto){
     double endpoint[2];
     if(!PathGetEndPoint(ia,endpoint)){
      return FALSE;
     }
     x0=endpoint[0];
     y0=endpoint[1];
    } else {
     PathMoveTo(ia,x0,y0);
    }
    if(x0==x&&y0==y)
     return TRUE;
    rx = abs(rx);
    ry = abs(ry);
    if(rx==0||ry==0){
//     startAngle=sweepAngle=atan2(y - y0, x - x0);
     PathLineTo(ia,x,y);
     return TRUE;
    }
    // Compute 1/2 distance between current and final point
    dx2 = (x0 - x) / 2.0;
    dy2 = (y0 - y) / 2.0;
    phi = fmod(phi,M_PITIMESTWO);
    // Compute (x1, y1)
    cosPhiR = cos(phi);
    sinPhiR = sin(phi);
    large_arc=TOBOOL(large_arc);
    sweep=TOBOOL(sweep);
    x1 = cosPhiR * dx2 + sinPhiR * dy2;
    y1 = -sinPhiR * dx2 + cosPhiR * dy2;
    // Make sure radii are large enough
    rx_sq = rx * rx;
    ry_sq = ry * ry;
    x1_sq = x1 * x1;
    y1_sq = y1 * y1;
    radius_check = (x1_sq / rx_sq) + (y1_sq / ry_sq);
    if (radius_check > 1){
        rx *= sqrt(radius_check);
        ry *= sqrt(radius_check);
        rx_sq = rx * rx;
        ry_sq = ry * ry;
    }
    // Step 2: Compute (cx1, cy1)
    sign = (large_arc == sweep) ? -1 : 1;
    sq = ((rx_sq * ry_sq) - (rx_sq * y1_sq) - (ry_sq * x1_sq)) /
        ((rx_sq * y1_sq) + (ry_sq * x1_sq));
    sq = (sq < 0) ? 0 : sq;
    coef = (sign * sqrt(sq));
    cx1 = coef * ((rx * y1) / ry);
    cy1 = coef * -((ry * x1) / rx);
    //   Step 3: Compute (cx, cy) from (cx1, cy1)
    sx2 = (x0 + x) / 2.0;
    sy2 = (y0 + y) / 2.0;
    cx = sx2 + (cosPhiR * cx1 - sinPhiR * cy1);
    cy = sy2 + (sinPhiR * cx1 + cosPhiR * cy1);
    //   Step 4: Compute angle start and angle extent
    ux = (x1 - cx1) / rx;
    uy = (y1 - cy1) / ry;
    vx = (-x1 - cx1) / rx;
    vy = (-y1 - cy1) / ry;
    n = atan2(uy, ux);
    theta=(n<0) ?  M_PITIMESTWO+n : n;
    p = atan2(vy, vx);
    delta=(p<n) ?  M_PITIMESTWO+p-n : p-n;
    if (!sweep && delta > 0){
        delta -= M_PITIMESTWO;
    } else if (sweep && delta < 0){
        delta += M_PITIMESTWO;
    }
    return PathAddArcInternal(ia, ArcHaveEndCoords, cx, cy, rx, ry, theta, delta, phi,x0,y0,x,y);
}

BOOL PathAddArcMultiSweep(
 ITEMARRAY*ia,
 double centerX, double centerY,
 double radiusX, double radiusY,
 double start, double sweep,
 double phi    
){
 while(sweep>360){
  if(!PathAddArc(ia,centerX,centerY,radiusX,radiusY,start,360,phi)){
   return FALSE;
  }
  sweep-=360;
 }
 while(sweep<-360){
  if(!PathAddArc(ia,centerX,centerY,radiusX,radiusY,start,-360,phi)){
   return FALSE;
  }
  sweep+=360;
 }
 return PathAddArc(ia,centerX,centerY,radiusX,radiusY,start,sweep,phi);
}

BOOL PathArcTo(
 ITEMARRAY*ia,
 double centerX, double centerY,
 double radiusX, double radiusY,
 double start, double sweep,//in degrees
 double phi//in degrees
){
 // Will check for previous move to and draw a line to the starting point
 return PathAddArcInternal(ia,
   ArcDrawLine|ArcCheckEndPoint,centerX,centerY,
   radiusX,radiusY,
   start*DEGTORAD,sweep*DEGTORAD,
   phi*DEGTORAD,0,0,0,0
 );
}


BOOL PathAddArc(
 ITEMARRAY*ia,
 double centerX, double centerY,
 double radiusX, double radiusY,
 double start, double sweep,//in degrees
 double phi//in degrees
){
 // Will move to starting point before drawing
 return PathAddArcInternal(ia,ArcMoveTo,centerX,centerY,
   radiusX,radiusY,
   start*DEGTORAD,sweep*DEGTORAD,
   phi*DEGTORAD,0,0,0,0
 );
}

BOOL PathSvgArcTo(
 ITEMARRAY*ia,
 double x, double y,
 double rx, double ry, double phi,//in degrees
 BOOL largeArc, BOOL sweep
){
 return PathAddSvgArcInternal(ia,TRUE,0,0,x,y,rx,ry,phi*DEGTORAD,largeArc,sweep);
}
BOOL PathAddSvgArc(
 ITEMARRAY*ia,
 double x0, double x1,
 double x, double y,
 double rx, double ry, double phi,//in degrees
 BOOL largeArc, BOOL sweep
){
 return PathAddSvgArcInternal(ia,FALSE,x0,x1,x,y,rx,ry,phi*DEGTORAD,largeArc,sweep);
}

BOOL EllipseContainsPoint(RECTDBL*rc,double x,double y){
 double w,dx,h,dy;
 w=RectWidth(rc);
 h=RectHeight(rc);
 if(w<=0||h<=0)return FALSE;
 dx=(x-rc->left)/w-0.5;
 dy=(y-rc->top)/h-0.5;
 return (dx*dx+dy*dy)<0.25;
}
BOOL EllipseContainsRect(RECTDBL*rc,double x0,double y0,
 double x1,double y1){
 return EllipseContainsPoint(rc,x0,y0)
  &&EllipseContainsPoint(rc,x0,y1)
  &&EllipseContainsPoint(rc,x1,y0)
  &&EllipseContainsPoint(rc,y1,y1);
}
//     Path segments

typedef struct{
 DWORD stage;
 ITEMARRAY *path;
} PATHSEGITERATOR;
static int PathSegmentsNonZero_GetWinding(LPVOID h){
 return WIND_NONZERO; 
}
static int PathSegmentsEvenOdd_GetWinding(LPVOID h){
 return WIND_EVENODD; 
}
static LPVOID PathSegments_Init(LPVOID data){
 PATHSEGITERATOR *iter;
 ITEMARRAY *ia=data;
 if(!ISARRAY(ia,PATHSEGMENT))return NULL;
 iter=malloc(sizeof(PATHSEGITERATOR));
 if(iter){
  iter->stage=0;
  iter->path=data;
 }
 return iter;
}
static void PathSegments_Free(LPVOID h){
 free(h);
}
static int PathSegments_Next(LPVOID h, PATHSEGMENT *ps){
 PATHSEGITERATOR *rc=(PATHSEGITERATOR *)h;
 PATHSEGMENT *pps;
 if(!h||!ps)return -1;
 if(rc->stage>=rc->path->length)
  return 0;
 pps=ItemArrayPtr(rc->path,rc->stage);
 *ps=*pps;
 rc->stage++;
 return 1;
}

PATHITERATOR PathSegmentsNonZero={
 PathSegments_Init,
 PathSegments_Next,
 PathSegments_Free,
 PathSegmentsNonZero_GetWinding
};

PATHITERATOR PathSegmentsEvenOdd={
 PathSegments_Init,
 PathSegments_Next,
 PathSegments_Free,
 PathSegmentsEvenOdd_GetWinding
};

//     Full path segments

void FullPathContextInit(FULLPATHCONTEXT *fpc){
 fpc->curx=0.0;
 fpc->cury=0.0;
 fpc->movx=0.0;
 fpc->movy=0.0;
 fpc->state=0;
}


BOOL FullPathSegmentTransform(FULLPATHSEGMENT*ps,MATRIX*xfm){
 int numpts,i;
 double x,y;
 if(!ps)
  return FALSE;
 if(!xfm)
  return TRUE;
 switch (ps->type){
 case PATH_MOVETO:
 case PATH_LINETO:
 case PATH_CLOSE:
  MatrixTransformPoints(xfm,ps->coords,2);
  return TRUE;
 case PATH_QUADTO:
  MatrixTransformPoints(xfm,ps->coords,3);
  return TRUE;
 case PATH_CUBICTO:
  MatrixTransformPoints(xfm,ps->coords,4);
  return TRUE;
 default:
  return FALSE;
 }
}

void FullToPathSegment(
  FULLPATHSEGMENT *fullSegment,
  PATHSEGMENT *pathSegment
){
  pathSegment->type=fullSegment->type;
  CopyMemory(pathSegment->coords,&fullSegment->coords[2],
      sizeof(double)*6);
}

typedef enum {
 PathNewMoveTo=0x01,
 PathEndOfPath=0x02,
 PathMask=0x07,
 PathHaveLineTo=0x08
} CLOSEDPATHFLAGS;

int ClosedPathIteratorNext(
  PATHITERATOR *iter,
  LPVOID h,
  FULLPATHCONTEXT *fpc,
  FULLPATHSEGMENT *fps,
  MATRIX *matrix
){
 PATHSEGMENT ps;
 int ret;
 if((fpc->state&PathMask)==PathNewMoveTo){
  // Return new move to
  fps->type=PATH_MOVETO;
  fps->coords[0]=fpc->curx;
  fps->coords[1]=fpc->cury;
  fps->coords[2]=fpc->movx;
  fps->coords[3]=fpc->movy;
  fpc->state&=~PathMask;
  return 1;
 } else if((fpc->state&PathMask)==PathEndOfPath){
  // Signal end of path
  fpc->state&=~PathMask;
  return 0;
 }
 // Clear closed path flag
 ret=iter->Next(h,&ps);
 if(ret<0)return ret;
 if(ret==0){
  if(fpc->state&PathHaveLineTo){
   fps->type=PATH_CLOSE;
   fps->coords[0]=fpc->curx;
   fps->coords[1]=fpc->cury;
   fps->coords[2]=fpc->movx;
   fps->coords[3]=fpc->movy;
   fpc->state&=~PathHaveLineTo;
   fpc->state|=PathEndOfPath;
   return 1;
  }
  // End of path
  return 0;
 }
 if(matrix)
  PathSegmentTransform(&ps,matrix);
 switch(ps.type){
   case PATH_MOVETO:
    if(fpc->state&PathHaveLineTo){
     fps->type=PATH_CLOSE;
     fps->coords[0]=fpc->curx;
     fps->coords[1]=fpc->cury;
     fps->coords[2]=fpc->curx=fpc->movx;
     fps->coords[3]=fpc->cury=fpc->movy;
     fpc->movx=ps.coords[0];
     fpc->movy=ps.coords[1];
     // Set closed path and move-to flags
     fpc->state&=~PathHaveLineTo;
     fpc->state|=PathNewMoveTo;
     return 1;
    } else {
     fps->type=PATH_MOVETO;
     fps->coords[0]=fpc->curx;
     fps->coords[1]=fpc->cury;
     fps->coords[2]=ps.coords[0];
     fps->coords[3]=ps.coords[1];
     fpc->movx=fpc->curx=ps.coords[0];
     fpc->movy=fpc->cury=ps.coords[1];
     fpc->state&=~PathHaveLineTo;
     break;
    }
   case PATH_LINETO:
    fps->type=PATH_LINETO;
    fps->coords[0]=fpc->curx;
    fps->coords[1]=fpc->cury;
    fpc->curx=fps->coords[2]=ps.coords[0];
    fpc->cury=fps->coords[3]=ps.coords[1];
    fpc->state|=PathHaveLineTo;
    break;
   case PATH_CLOSE:
    fps->type=PATH_CLOSE;
    fps->coords[0]=fpc->curx;
    fps->coords[1]=fpc->cury;
    fpc->curx=fps->coords[2]=fpc->movx;
    fpc->cury=fps->coords[3]=fpc->movy;
    fpc->state&=~PathHaveLineTo;
    break;
   case PATH_QUADTO:
    fps->type=PATH_QUADTO;
    fps->coords[0]=fpc->curx;
    fps->coords[1]=fpc->cury;
    fps->coords[2]=ps.coords[0];
    fps->coords[3]=ps.coords[1];
    fpc->curx=fps->coords[4]=ps.coords[2];
    fpc->cury=fps->coords[5]=ps.coords[3];
    fpc->state|=PathHaveLineTo;
    break;
   case PATH_CUBICTO:
    fps->type=PATH_CUBICTO;
    fps->coords[0]=fpc->curx;
    fps->coords[1]=fpc->cury;
    fps->coords[2]=ps.coords[0];
    fps->coords[3]=ps.coords[1];
    fps->coords[4]=ps.coords[2];
    fps->coords[5]=ps.coords[3];
    fpc->curx=fps->coords[6]=ps.coords[4];
    fpc->cury=fps->coords[7]=ps.coords[5];
    fpc->state|=PathHaveLineTo;
    break;
   default:
    return -1;
 }
 fpc->state|=16;
 return 1;
}

int FullPathIteratorNext(
  PATHITERATOR *iter,
  LPVOID h,
  FULLPATHCONTEXT *fpc,
  FULLPATHSEGMENT *fps,
  MATRIX *matrix
){
 PATHSEGMENT ps;
 int ret=iter->Next(h,&ps);
 if(ret<=0)return ret;
 if(matrix){
  PathSegmentTransform(&ps,matrix);
 }
 switch(ps.type){
   case PATH_MOVETO:
    fps->type=PATH_MOVETO;
    fps->coords[0]=fpc->curx;
    fps->coords[1]=fpc->cury;
    fps->coords[2]=ps.coords[0];
    fps->coords[3]=ps.coords[1];
    fpc->movx=fpc->curx=ps.coords[0];
    fpc->movy=fpc->cury=ps.coords[1];
    break;
   case PATH_LINETO:
    fps->type=PATH_LINETO;
    fps->coords[0]=fpc->curx;
    fps->coords[1]=fpc->cury;
    fpc->curx=fps->coords[2]=ps.coords[0];
    fpc->cury=fps->coords[3]=ps.coords[1];
    break;
   case PATH_CLOSE:
    fps->type=PATH_CLOSE;
    fps->coords[0]=fpc->curx;
    fps->coords[1]=fpc->cury;
    fpc->curx=fps->coords[2]=fpc->movx;
    fpc->cury=fps->coords[3]=fpc->movy;
    break;
   case PATH_QUADTO:
    fps->type=PATH_QUADTO;
    fps->coords[0]=fpc->curx;
    fps->coords[1]=fpc->cury;
    fps->coords[2]=ps.coords[0];
    fps->coords[3]=ps.coords[1];
    fpc->curx=fps->coords[4]=ps.coords[2];
    fpc->cury=fps->coords[5]=ps.coords[3];
    break;
   case PATH_CUBICTO:
    fps->type=PATH_CUBICTO;
    fps->coords[0]=fpc->curx;
    fps->coords[1]=fpc->cury;
    fps->coords[2]=ps.coords[0];
    fps->coords[3]=ps.coords[1];
    fps->coords[4]=ps.coords[2];
    fps->coords[5]=ps.coords[3];
    fpc->curx=fps->coords[6]=ps.coords[4];
    fpc->cury=fps->coords[7]=ps.coords[5];
    break;
   default:
    return -1;
 }
 return 1;
}


typedef struct{
 int phase;
 BOOL last;
 double left[8];
 double right[8];
} CURVEFLATSTATUS;

typedef struct{
 double startcurve[8];
 CURVEFLATSTATUS *stack;
 int mode;
 double flatness;
 double *thiscurve;
 LONG stackpos;
 LONG stacksize;
 double coords[8];
 DWORD numpoints;
 DWORD iPoint;
} CURVEFLATTENER;

static void CurveFlattenerPush(CURVEFLATTENER *f){
 f->stackpos++;
 if(!f->stack || f->stackpos>=f->stacksize){ 
  //If, e.g., stackpos is now 2, a stack size of 2 or more is not
  //enough
  f->stacksize+=10;
  f->stack=realloc(f->stack,f->stacksize*sizeof(CURVEFLATSTATUS));
 }
 f->stack[f->stackpos].phase=0;
 f->stack[f->stackpos].last=FALSE;
}

static LPVOID CurveFlattenerCreate(double flatness, int flatnessmode){
 CURVEFLATTENER *f=xmalloc(sizeof(CURVEFLATTENER));
 f->thiscurve=NULL;
 if(flatnessmode==0){
  // Used in curve flattener
  f->flatness=flatness;
 } else {
  // Used in cubic-to-quad iterator
  f->flatness=abs(flatness)/0.09622504486493;
 }
 f->stack=NULL;
 f->stacksize=0;
 f->iPoint=0;
 f->numpoints=0;
 f->stackpos=-1;
 return f;
}

static void CurveFlattenerFree(LPVOID pf){
 CURVEFLATTENER *f=(CURVEFLATTENER*)pf;
 if(f){
  free(f->stack); 
  ZeroMemory(f,sizeof(CURVEFLATTENER));
  free(f);
 }
}

static void CurveFlattenerReset(LPVOID pf, FULLPATHSEGMENT *fps){
 CURVEFLATTENER *f=(CURVEFLATTENER*)pf;
 if(fps->type==PATH_QUADTO){
  f->stackpos=-1;
  CopyMemory(f->startcurve,fps->coords,sizeof(double)*6);
  f->thiscurve=f->startcurve;
  f->mode=0;
  CurveFlattenerPush(f); // to allocate first stack position
 } else if(fps->type==PATH_CUBICTO){
  f->stackpos=-1;
  CopyMemory(f->startcurve,fps->coords,sizeof(double)*8);
  f->thiscurve=f->startcurve;
  f->mode=1;
  CurveFlattenerPush(f); // to allocate first stack position
 }
 f->stack[0].last=TRUE;
 f->iPoint=0;
 f->numpoints=0;
}


// See http://www.antigrain.com/research/adaptive_bezier/
BOOL QuadFlatEnough(double *curve, double tolerance, double *endpt){
    double collinearTolerance=1e-25;
    double angleTolerance=12*DEGTORAD;
    double x12   = (curve[0] + curve[2]) / 2;                
    double y12   = (curve[1] + curve[3]) / 2;
    double x23   = (curve[2] + curve[4]) / 2;
    double y23   = (curve[3] + curve[5]) / 2;
    double midpointX  = (x12 + x23) / 2;
    double midpointY  = (y12 + y23) / 2;
    double dx=curve[4]-curve[0];
    double dy=curve[5]-curve[1];
    double d2=(curve[2]-curve[4])*dy-(curve[3]-curve[5])*dx;
    d2=abs(d2);
    if(d2>collinearTolerance){
     // Regular, not collinear
     if(d2*d2<tolerance*(dx*dx+dy*dy)){
      double a23 = atan2(curve[5] - curve[3], curve[4] - curve[2]);
      double da1 = fabs(a23 - atan2(curve[3] - curve[1], curve[2] - curve[0]));
      if(da1 >= M_PI) da1 = M_PITIMESTWO - da1;
      if(da1<angleTolerance){
       endpt[0]=midpointX;
       endpt[1]=midpointY;
       return 1;
      }
     }
    } else {
     // Collinear
     double tol=2/(0.5/sqrt(tolerance));
     if(    fabs(curve[0] + curve[4] - curve[2] - curve[2]) +
            fabs(curve[1] + curve[5] - curve[3] - curve[3]) <= tol){
      endpt[0]=midpointX;
      endpt[1]=midpointY;
      return 1;
     }
    }
    return 0;
}
// See http://www.antigrain.com/research/adaptive_bezier/
BOOL CubicFlatEnough(double *curve, double tolerance, double *endpt){
    double collinearTolerance=1e-25;
    double distanceTolerance=tolerance;
    double angleTolerance=12*DEGTORAD;
    double cuspLimit=M_PI-(1*DEGTORAD);
    double midpointX = (curve[2] + curve[4]) / 2;
    double midpointY = (curve[3] + curve[5]) / 2;
    double dx=curve[6]-curve[0];
    double dy=curve[7]-curve[1];
    double d2=(curve[2]-curve[6])*dy-(curve[3]-curve[7])*dx;
    double d3=(curve[4]-curve[6])*dy-(curve[5]-curve[7])*dx;
    d2=abs(d2);
    d3=abs(d3);
    if(d2>collinearTolerance && d3>collinearTolerance){
     // Regular, not collinear
     if((d2+d3)*(d2+d3)<distanceTolerance*(dx*dx+dy*dy)){
      if(angleTolerance<1e-25){
       endpt[0]=midpointX;
       endpt[1]=midpointY;
       return 1;
      }
      double a23 = atan2(curve[5] - curve[3], curve[4] - curve[2]);
      double da1 = fabs(a23 - atan2(curve[3] - curve[1], curve[2] - curve[0]));
      double da2 = fabs(atan2(curve[7] - curve[5], curve[6] - curve[4]) - a23);
      if(da1 >= M_PI) da1 = M_PITIMESTWO - da1;
      if(da2 >= M_PI) da2 = M_PITIMESTWO - da2;
      if(da1+da2<angleTolerance){
       endpt[0]=midpointX;
       endpt[1]=midpointY;
       return 1;
      } else if(cuspLimit!=0 && da1>cuspLimit){
       endpt[0]=curve[2];
       endpt[1]=curve[3];
       return 1;
      } else if(cuspLimit!=0 && da2>cuspLimit){
       endpt[0]=curve[4];
       endpt[1]=curve[5];
       return 1;
      }
     }
    } else if(d2>collinearTolerance){
     if(d2*d2<distanceTolerance*(dx*dx+dy*dy)){
      if(angleTolerance<1e-25){
       endpt[0]=midpointX;
       endpt[1]=midpointY;
       return 1;
      }
      double a23 = atan2(curve[5] - curve[3], curve[4] - curve[2]);
      double da1 = fabs(a23 - atan2(curve[3] - curve[1], curve[2] - curve[0]));
      if(da1 >= M_PI) da1 = M_PITIMESTWO - da1;
      if(da1<angleTolerance){
       endpt[0]=curve[2];
       endpt[1]=curve[3];
       endpt[2]=curve[4];
       endpt[3]=curve[5];
       return 2;
      } else if(cuspLimit!=0 && da1>cuspLimit){
       endpt[0]=curve[2];
       endpt[1]=curve[3];
       return 1;
      }
     }
    } else if(d3>collinearTolerance){
     if(d3*d3<distanceTolerance*(dx*dx+dy*dy)){
      if(angleTolerance<1e-25){
       endpt[0]=midpointX;
       endpt[1]=midpointY;
       return 1;
      }
      double a23 = atan2(curve[5] - curve[3], curve[4] - curve[2]);
      double da2 = fabs(atan2(curve[7] - curve[5], curve[6] - curve[4]) - a23);
      if(da2 >= M_PI) da2 = M_PITIMESTWO - da2;
      if(da2<angleTolerance){
       endpt[0]=curve[2];
       endpt[1]=curve[3];
       endpt[2]=curve[4];
       endpt[3]=curve[5];
       return 2;
      } else if(cuspLimit!=0 && da2>cuspLimit){
       endpt[0]=curve[4];
       endpt[1]=curve[5];
       return 1;
      }
     }
    } else {
     // Collinear
     double tol=4/(0.5/sqrt(tolerance));
     if(    fabs(curve[0] + curve[4] - curve[2] - curve[2]) +
            fabs(curve[1] + curve[5] - curve[3] - curve[3]) +
            fabs(curve[2] + curve[6] - curve[4] - curve[4]) +
            fabs(curve[3] + curve[7] - curve[5] - curve[5]) <= tol){
      endpt[0]=midpointX;
      endpt[1]=midpointY;
      return 1;
     }
    }
    return 0;
}
static BOOL CurveFlattenerNextLine(CURVEFLATTENER *f,double *line){
 if(!f||!f->thiscurve){
  return FALSE;
 }
 while(1){
  if(f->numpoints>0){
   line[0]=f->coords[f->iPoint];
   line[1]=f->coords[f->iPoint+1];
   line[2]=f->coords[f->iPoint+2];
   line[3]=f->coords[f->iPoint+3];
   f->iPoint+=2;
   if(f->iPoint>=f->numpoints){
    f->numpoints=0;
    if(f->stackpos==0){
     // Break, we're done
     f->thiscurve=NULL;
    } else {
     // Pop the flattener stack
     f->stackpos--;
    }
   }
   return TRUE;
  }
  if(f->stack[f->stackpos].phase==0){
   // Phase 0: Check whether the curve is close to a line
   BOOL closeenough=FALSE;
   ASSERT(f->thiscurve);
   if(f->mode==0){ // Quadratic curve
    int points=QuadFlatEnough(f->thiscurve,f->flatness,&f->coords[2]);
    if(points>0){
     f->coords[0]=f->thiscurve[0];
     f->coords[1]=f->thiscurve[1];
     if(f->stack[f->stackpos].last){
      f->numpoints=(points+1)*2;
      f->coords[f->numpoints]=f->thiscurve[4];
      f->coords[f->numpoints+1]=f->thiscurve[5];
     } else {
      f->numpoints=points*2;
     }
     f->iPoint=0;
     continue;
    }
   } else { // Cubic curve
    int points=CubicFlatEnough(f->thiscurve,f->flatness,&f->coords[2]);
    if(points>0){
     f->coords[0]=f->thiscurve[0];
     f->coords[1]=f->thiscurve[1];
     if(f->stack[f->stackpos].last){
      f->numpoints=(points+1)*2;
      f->coords[f->numpoints]=f->thiscurve[6];
      f->coords[f->numpoints+1]=f->thiscurve[7];
     } else {
      f->numpoints=points*2;
     }
     f->iPoint=0;
     continue;
    }
   }
   f->stack[f->stackpos].phase=1;
  }
  if(f->stack[f->stackpos].phase==1){
   // Phase 1: Subdivide the curve
   CURVEFLATSTATUS *status=&f->stack[f->stackpos];
   ASSERT(f->thiscurve);
   if(f->mode==0){ // Quadratic curve
    QuadSubdivide(f->thiscurve,status->left,status->right,0.5);
   } else { // Cubic curve
    CubicSubdivide(f->thiscurve,status->left,status->right,0.5);
   }
   status->phase=2;
  }
  if(f->stack[f->stackpos].phase==2){
   // Phase 2: Push the left curve
   f->thiscurve=f->stack[f->stackpos].left;
   f->stack[f->stackpos].phase=3;
   CurveFlattenerPush(f);
   continue;
  }
  if(f->stack[f->stackpos].phase==3){
   // Phase 3: Push the right curve
   f->thiscurve=f->stack[f->stackpos].right;
   f->stack[f->stackpos].phase=4;
   CurveFlattenerPush(f);
   f->stack[f->stackpos].last=f->stack[f->stackpos-1].last;
   continue;
  }
  if(f->stack[f->stackpos].phase==4){
   // Phase 4: Pop the stack
   if(f->stackpos==0){
    // Break, we're done
    f->thiscurve=NULL;
    return FALSE;
   } else {
    // Pop the flattener stack
    f->stackpos--;
   }
  }
 }
 return FALSE;
}


static void MidpointCubicToQuad(double *cubic, double *quad){
    double midX=(3*cubic[4] - cubic[6] + 3*cubic[2] - cubic[0])/4;
    double midY=(3*cubic[5] - cubic[7] + 3*cubic[3] - cubic[1])/4;
    quad[0]=cubic[0];
    quad[1]=cubic[1];
    quad[2]=midX;
    quad[3]=midY;
    quad[4]=cubic[6];
    quad[5]=cubic[7];
}

//  Using algorithm for adaptive subdivision found at:
//  http://caffeineowl.com/graphics/2d/vectorial/cubic2quad01.html
static BOOL CubicIteratorNextQuadAlt(CURVEFLATTENER *f,double *quad){
 if(!f||!f->thiscurve){
  return FALSE;
 }
 while(1){
  double tmax=0.0;
  if(f->stack[f->stackpos].phase==0){
   // Phase 0: Check whether the curve is close enough
   // to a quadratic curve
   double *thiscurve=f->thiscurve;
   double dx=thiscurve[6] - 3*thiscurve[4] + 3*thiscurve[2] - thiscurve[0];
   double dy=thiscurve[7] - 3*thiscurve[5] + 3*thiscurve[3] - thiscurve[1];
   double dist=sqrt(dx*dx+dy*dy)/2;
   if(dist==0.0){
    tmax=1.0;
   } else {
    tmax=pow(f->flatness/dist,1.0/3.0);
   }
   if(tmax>=1.0){
    MidpointCubicToQuad(thiscurve,quad);
    if(f->stackpos==0){
     // Break, we're done
     f->thiscurve=NULL;
     return TRUE;
    } else {
     // Pop the flattener stack
     f->stackpos--;
     return TRUE;
    }
   }
   f->stack[f->stackpos].phase=1;
  }
  if(f->stack[f->stackpos].phase==1){
   // Phase 1: Subdivide the curve
   CURVEFLATSTATUS *status=&f->stack[f->stackpos];
   ASSERT(f->thiscurve);
   // Subdividing at tmax, found at phase 0
   CubicSubdivide(f->thiscurve,status->left,status->right,tmax);
   // Return this curve
   MidpointCubicToQuad(status->left,quad);
   status->phase=2;
   return TRUE;
  }
  if(f->stack[f->stackpos].phase==2){
   // Phase 2: Push the right curve
   f->thiscurve=f->stack[f->stackpos].right;
   f->stack[f->stackpos].phase=3;
   CurveFlattenerPush(f);
   continue;
  }
  if(f->stack[f->stackpos].phase==3){
   // Phase 3: Pop the stack
   if(f->stackpos==0){
    // Break, we're done
    f->thiscurve=NULL;
    return FALSE;
   } else {
    // Pop the flattener stack
    f->stackpos--;
   }
  }
 }
 return FALSE;
}
// Using adaptive symmetric division, also found at:
//  http://caffeineowl.com/graphics/2d/vectorial/cubic2quad01.html
static BOOL CubicIteratorNextQuad(CURVEFLATTENER *f,double *quad){
 if(!f||!f->thiscurve){
  return FALSE;
 }
 while(1){
  double tmax=0.0;
  if(f->stack[f->stackpos].phase==0){
   // Phase 0: Check whether the curve is close enough
   // to a quadratic curve
   double *thiscurve=f->thiscurve;
   double dx=thiscurve[6] - 3*thiscurve[4] + 3*thiscurve[2] - thiscurve[0];
   double dy=thiscurve[7] - 3*thiscurve[5] + 3*thiscurve[3] - thiscurve[1];
   double dist=sqrt(dx*dx+dy*dy)/2;
   if(dist==0.0){
    tmax=1.0;
   } else {
    tmax=pow(f->flatness/dist,1.0/3.0);
   }
   if(tmax>=1.0){
    MidpointCubicToQuad(thiscurve,quad);
    if(f->stackpos==0){
     // Break, we're done
     f->thiscurve=NULL;
     return TRUE;
    } else {
     // Pop the flattener stack
     f->stackpos--;
     return TRUE;
    }
   }
   f->stack[f->stackpos].phase=1;
  }
  if(f->stack[f->stackpos].phase==1){
   // Phase 1: Subdivide the curve (Can reach here only
   // directly after phase 0; this phase can modify "thiscurve")
   CURVEFLATSTATUS *status=&f->stack[f->stackpos];
   ASSERT(f->thiscurve);
   if(tmax<0.5){
    double tmax2=1-tmax;
    // Subdividing into left and mid/right
    CubicSubdivide(f->thiscurve,f->thiscurve,status->left,tmax);
    // Subdividing mid/right into mid and right
    tmax2=(tmax2-tmax)/(1.0-tmax);
    CubicSubdivide(status->left,status->left,status->right,tmax2);
    // Return the left curve
    MidpointCubicToQuad(f->thiscurve,quad);
    status->phase=2;
   } else {
    // Halving the curve
    CubicSubdivide(f->thiscurve,status->left,status->right,0.5);
    // Return the left curve
    MidpointCubicToQuad(status->left,quad);
    // Skip to phase 3, no middle curve to push
    status->phase=3;
   }
   return TRUE;
  }
  if(f->stack[f->stackpos].phase==2){
   // Phase 2: Push the middle curve
   f->thiscurve=f->stack[f->stackpos].left;
   f->stack[f->stackpos].phase=3;
   CurveFlattenerPush(f);
   continue;
  }
  if(f->stack[f->stackpos].phase==3){
   // Phase 3: Return the right curve
   f->thiscurve=f->stack[f->stackpos].right;
   f->stack[f->stackpos].phase=4;
   MidpointCubicToQuad(f->thiscurve,quad);
   return TRUE;
  }
  if(f->stack[f->stackpos].phase==4){
   // Phase 4: Pop the stack
   if(f->stackpos==0){
    // Break, we're done
    f->thiscurve=NULL;
    return FALSE;
   } else {
    // Pop the flattener stack
    f->stackpos--;
   }
  }
 }
 return FALSE;
}
typedef struct{
 PATHITERATOR *iter;
 LPVOID data;
 LPVOID handle;
 MATRIX matrix;
 BOOL haveMatrix;
 FULLPATHCONTEXT fpc;
 CURVEFLATTENER *flat;
} FLATPATHINTERNAL;

static LPVOID FlatPathIterator_InitInternal(FLATPATHPARAMS *fpp, int mode){
 FLATPATHINTERNAL *fpi;
 if(!fpp||!fpp->iterator||fpp->tolerance==0.0)return NULL;
 fpi=malloc(sizeof(FLATPATHINTERNAL));
 if(!fpi)return NULL;
 FullPathContextInit(&fpi->fpc);
 fpi->flat=CurveFlattenerCreate(fpp->tolerance,mode);
 if(!fpi->flat){
  free(fpi);
  return NULL;
 }
 fpi->iter=fpp->iterator;
 fpi->data=fpp->data;
 fpi->handle=fpi->iter->Init(fpi->data);
 if(!fpi->handle){
  CurveFlattenerFree(fpi->flat);
  free(fpi);
  return NULL;
 }
 fpi->haveMatrix=FALSE;
 if(fpp->matrix){
  fpi->haveMatrix=TRUE;
  fpi->matrix=*(fpp->matrix);
 }
 return (LPVOID)fpi;
}
static LPVOID FlatPathIterator_Init(FLATPATHPARAMS *fpp){
 return FlatPathIterator_InitInternal(fpp,0);
}
static LPVOID CubicToQuadIterator_Init(FLATPATHPARAMS *fpp){
 return FlatPathIterator_InitInternal(fpp,1);
}


static void FlatPathIterator_Free(LPVOID h){
 FLATPATHINTERNAL *fpi=(FLATPATHINTERNAL *)h;
 if(fpi){
  fpi->iter->Free(fpi->handle);
  CurveFlattenerFree(fpi->flat);
  free(fpi);
 }
}
static int FlatPathIterator_GetWinding(LPVOID h){
 FLATPATHINTERNAL *fpi=(FLATPATHINTERNAL *)h;
 if(!fpi)return WIND_NONZERO;
 return fpi->iter->GetWinding(fpi->handle);
}
static int FlatPathIterator_Next(LPVOID h, PATHSEGMENT *seg){
 FLATPATHINTERNAL *fpi=(FLATPATHINTERNAL *)h;
 FULLPATHSEGMENT fps;
 double line[4];
 int ret;
 if(!fpi||!seg)return -1;
 while(1){
  if(CurveFlattenerNextLine(fpi->flat,line)){
   seg->type=PATH_LINETO;
   seg->coords[0]=line[2];
   seg->coords[1]=line[3];
   return 1;
  }
  ret=FullPathIteratorNext(fpi->iter,fpi->handle,&fpi->fpc,&fps,
   fpi->haveMatrix ? &fpi->matrix : NULL);
  DebugOut("Returning seg",__LINE__);
  if(ret<=0)return ret;
  if(fps.type==PATH_LINETO|| 
     fps.type==PATH_CLOSE||
     fps.type==PATH_MOVETO){
   FullToPathSegment(&fps,seg);
   return 1;
  } else if(fps.type==PATH_CUBICTO||fps.type==PATH_QUADTO){
   CurveFlattenerReset(fpi->flat,&fps);
  }
 }
 return -1;
}

static int CubicToQuadIterator_Next(LPVOID h, PATHSEGMENT *seg){
 FLATPATHINTERNAL *fpi=(FLATPATHINTERNAL *)h;
 FULLPATHSEGMENT fps;
 double line[6];
 int ret;
 if(!fpi||!seg)return -1;
 while(1){
  if(CubicIteratorNextQuad(fpi->flat,line)){
   seg->type=PATH_QUADTO;
   seg->coords[0]=line[2];
   seg->coords[1]=line[3];
   seg->coords[2]=line[4];
   seg->coords[3]=line[5];
   return 1;
  }
  ret=FullPathIteratorNext(fpi->iter,fpi->handle,&fpi->fpc,&fps,
   fpi->haveMatrix ? &fpi->matrix : NULL);
  if(ret<=0)return ret;
  if(fps.type==PATH_LINETO|| 
     fps.type==PATH_CLOSE||
     fps.type==PATH_MOVETO||
     fps.type==PATH_QUADTO){
   FullToPathSegment(&fps,seg);
   return 1;
  } else if(fps.type==PATH_CUBICTO){
   CurveFlattenerReset(fpi->flat,&fps);
  }
 }
 return -1;
}
PATHITERATOR FlatPathIterator={
 FlatPathIterator_Init,
 FlatPathIterator_Next,
 FlatPathIterator_Free,
 FlatPathIterator_GetWinding
};

PATHITERATOR CubicToQuadIterator={
 CubicToQuadIterator_Init,
 CubicToQuadIterator_Next,
 FlatPathIterator_Free,
 FlatPathIterator_GetWinding
};

#define HORNERQUAD(x,a,b,c) \
  ( (a) + (x) * ( (b) + (x) * (c) ))
#define HORNERCUBIC(x,a,b,c,d) \
  ( (a) + (x) * ( (b) + (x) * ( (c) + (x) * (d) ) ))


int SolveLinear(double c[2], double s[1])
{
if (abs(c[1])<=1e-9)
    return 0;
s[0] = - c[0] / c[1];
if(ISONE(s[0]))  s[0]=1.0;
if(ISZERO(s[0])) s[1]=0.0;
return 1;
}

int SolveQuadratic(double eqn[3], double s[2]){
 double discrim,a,b,c,temp,sqrtDiscrim;
 int num,i;
 c=eqn[0];
 b=eqn[1];
 a=eqn[2];
 // make sure we have a d2 equation
 if (abs(eqn[2])<=1e-9)
  return SolveLinear(eqn, s);
 discrim = b*b - 4.0f*a*c;
 if (discrim < 0){
  return 0;
 } else if(discrim == 0){
  s[0]=s[1]=-0.5 * b / a ;
  num=1;
 } else {
  sqrtDiscrim = sqrt(discrim);
  if (b < 0){
   temp = -0.5f * (b - sqrtDiscrim);
  } else {
   temp = -0.5f * (b + sqrtDiscrim);
  }
  s[0] = temp / a;
  s[1] = c / temp;
  num=2;
 }
 for(i=0;i<num;i++){
  if(ISZERO(s[i]))s[i]=0.0;
  if(ISONE(s[i]))s[i]=1.0;
 }
 return num;
}

static double NewtonFixup(double s, double c0, double c1, double c2, double c3){
 double c2_2 = c2 * 2.0;
 double c3_3 = c3 * 3.0;
 double last=s;
 double direction=0;
 BOOL first=TRUE;
 while(1){
  double hquad=HORNERQUAD(s,c1,c2_2,c3_3);
  double delta;
  int signDelta;
  if(hquad!=0){
   s = s - HORNERCUBIC(s,c0,c1,c2,c3) / hquad;
  } else {
   return s;
  }
  delta=s-last;
  signDelta=(delta<0) ? -1 : (delta==0 ? 0 : 1);
  if(signDelta==0){
   return s;
  }
  if(first){
   last=s;
   direction=signDelta;
   first=FALSE;
  } else {
   if(signDelta!=direction){
    return s+(last-s)/2.0;
   }
   last=s;
   direction=signDelta;
  }
 }
 return 0;
}

int SolveCubic(double c[4], double s[3])
{
int     i, num, j;
double  sub,
        A, B, C,
        sq_A, p, q,
        cb_p, D;
double c0=c[0];
double c1=c[1];
double c2=c[2];
double c3=c[3];
double c2_2;
double c3_3;
double third=1.0/3.0;
if(abs(c3)<=1e-9)
 return SolveQuadratic(c,s);
// normalize the equation:x ^ 3 + Ax ^ 2 + Bx  + C = 0
A = c2 / c3;
B = c1 / c3;
C = c0 / c3;
// substitute x = y - A / 3 to eliminate the quadric term: x^3 + px + q = 0
sq_A = A * A;
p = third * (-third * sq_A + B);
q = 0.5 * (2.0/27.0 * A *sq_A - third * A * B + C);
// use Cardano's formula
cb_p = p * p * p;
D = q * q + cb_p;
if (abs(D)<=1e-9) //Epsilon carefully selected for stability
    {
    if (abs(q)<=1e-9)
        {
        // one triple solution
        s[0] = 0.0;
        num = 1;
        }
    else
        {
        // one single and one double solution
        double u = CubeRoot(-q);
        s[0] = 2.0 * u;
        s[1] = - u;
        num = 2;
        }
    }
else
    if (D < 0.0)
        {
        // casus irreductibilis: three real solutions
        double phi = 1.0/3.0 * acos(-q / sqrt(-cb_p));
        double t = 2.0 * sqrt(-p);
        s[0] = t * cos(phi);
        s[1] = -t * cos(phi + M_PI / 3.0);
        s[2] = -t * cos(phi - M_PI / 3.0);
        num = 3;
        }
    else
        {
        // one real solution
        double sqrt_D = sqrt(D);
        double u = CubeRoot(sqrt_D + abs(q));
        if (q > 0.0)
            s[0] = - u + p / u ;
        else
            s[0] = u - p / u;
        num = 1;
        }
// resubstitute
sub = A*third;
for (i = 0; i < num; i++)
    s[i] -= sub;
// Newton iterations to stabilize roots
for(j=0;j<4;j++){
 for(i=0;i<num;i++){
   s[i]=NewtonFixup(s[i],c0,c1,c2,c3);
 }
}
for(i=0;i<num;i++){
 if(abs(s[i])<1e-8)s[i]=0.0;
 if(abs(s[i]-1.0)<1e-8)s[i]=1.0;
}
return num;
}

//    Quadratic curves
double QuadControlPolygonLength(double *qc){
 return PointDistance(qc[0],qc[1],qc[2],qc[3])
       +PointDistance(qc[2],qc[3],qc[4],qc[5]);
}
double QuadChordLength(double *qc){
 return PointDistance(qc[0],qc[1],qc[4],qc[5]);
}
double QuadCalcLength(double *qc, double tolerance){
 double ret=0.0,pl,cl;
 double *cs=malloc(sizeof(double)*6);
 int stack=0,stacklen=6;
 CopyMemory(cs,qc,sizeof(double)*6);
 while(1){
  cl=QuadChordLength(&cs[stack]);
  pl=QuadControlPolygonLength(&cs[stack]);
  if(pl-cl<=tolerance){
   ret+=(2.0*pl+cl)/3.0;
   if(stack==0)
    break;
   stack-=6;
  }
  else {
   if(stack+16>stacklen){
    stacklen+=16+48;
    cs=realloc(cs,stacklen*sizeof(double));
   }
   QuadSubdivide(&cs[stack],&cs[stack+6],&cs[stack],0.5);
   stack+=6;
  }
 }
 free(cs);
 return ret;
}

//    Cubic curves
double CubicControlPolygonLength(double *qc){
 return PointDistance(qc[0],qc[1],qc[2],qc[3])
+PointDistance(qc[2],qc[3],qc[4],qc[5])
+PointDistance(qc[4],qc[5],qc[6],qc[7]);
}
double CubicChordLength(double *qc){
 return PointDistance(qc[0],qc[1],qc[6],qc[7]);
}
double CubicCalcLength(double *qc, double tolerance){
 double ret=0.0,pl,cl;
 double *cs=malloc(sizeof(double)*8);
 int stack=0,stacklen=8;
 CopyMemory(cs,qc,sizeof(double)*8);
 while(1){
  cl=CubicChordLength(&cs[stack]);
  pl=CubicControlPolygonLength(&cs[stack]);
  if(pl-cl<=tolerance){
   ret+=0.5*pl+0.5*cl;
   if(stack==0)
    break;
   stack-=8;
  }
  else {
   if(stack+16>stacklen){
    stacklen+=16+48;
    cs=realloc(cs,stacklen*sizeof(double));
   }
   CubicSubdivide(&cs[stack],&cs[stack+8],&cs[stack],0.5);
   stack+=8;
  }
 }
 free(cs);
 return ret;
}
double EllipseArea(double radiusX, double radiusY, MATRIX *xfm){
 double area=abs((double)(radiusX * radiusY))*M_PI;
 if (xfm)
  area*=MatrixDeterminant(xfm);
 return area;
}
static BOOL RectAccumulate(RECTDBL *prc, double x, double y, BOOL first){
 if(first){
  SetRect(prc,x,y,x,y);
 } else {
  RectAddPoint(prc,x,y);
 }
 return FALSE;
}
static int CubicDerivatives(double *zeros,double cur,double cp0,double cp1,double end){
 zeros[0]=(cp0-cur)*3.0;
 zeros[1]=(cp1-cp0-cp0+cur)*6.0;
 zeros[2]=(end+(cp0-cp1)*3.0-cur)*3.0;
 return SolveQuadratic(zeros,zeros);
}
int PathGetBounds(PATHITERATOR*pcb,LPVOID data,MATRIX*xfm,RECTDBL*prc){
 FULLPATHSEGMENT fps;
 FULLPATHCONTEXT fpc;
 LPVOID param;
 int i;
 BOOL first=TRUE;
 if(!prc)return -1;
 SetRect(prc,0,0,0,0);
 if(!pcb)return -1;
 param=pcb->Init(data);
 if(!param)return -1;
 FullPathContextInit(&fpc);
 while(FullPathIteratorNext(pcb,param,&fpc,&fps,xfm)>0){
  switch (fps.type){
   case PATH_MOVETO:
   case PATH_LINETO:
    first=RectAccumulate(prc,fps.coords[2],fps.coords[3],first);
    break;
   case PATH_QUADTO:{
    double *c=fps.coords;
    double zeros[2];
    first=RectAccumulate(prc,c[4],c[5],first);
    zeros[0]=-(c[2]+c[2]-c[0]-c[0])/((c[0]+c[4]-c[2]-c[2])*2);
    zeros[1]=-(c[3]+c[3]-c[1]-c[1])/((c[1]+c[5]-c[3]-c[3])*2);
    for(i=0;i<2;i++){
     if(zeros[i]>0.0 && zeros[i]<1.0){
      // Finding the extreme for zeros[i]
      double t=zeros[i];
      double u=(1.0-t);
      double usq=u*u;
      double tsq=t*t;
      double x=c[0]*usq+ 2.0*c[2]*t*u+ c[4]*tsq;
      double y=c[1]*usq+ 2.0*c[3]*t*u+ c[5]*tsq;
      first=RectAccumulate(prc,x,y,first);
     }
    }
    break;
   }
   case PATH_CUBICTO:{
    double *c=fps.coords;
    double zeros[6];
    int count;
    first=RectAccumulate(prc,c[6],c[7],first);
    count=CubicDerivatives(zeros,c[0],c[2],c[4],c[6]);
    count+=CubicDerivatives(&zeros[count],c[1],c[3],c[5],c[7]);
    for(i=0;i<count;i++){
     if(zeros[i]>0.0 && zeros[i]<1.0){
      // Finding the extreme for zeros[i]
      double t=zeros[i];
      double u=(1.0-t);
      double usq=u*u;
      double ucb=usq*u;
      double tsq=t*t;
      double tcb=tsq*t;
      double x = c[0]*ucb+ 3.0*c[2]*t*usq+ 3.0*c[4]*tsq*u+ c[6]*tcb;
      double y = c[1]*ucb+ 3.0*c[3]*t*usq+ 3.0*c[5]*tsq*u+ c[7]*tcb;
      first=RectAccumulate(prc,x,y,first);
     }
    }
    break;
   }
  }
 }
 pcb->Free(param);
 return (first) ? 0 : 1;
}

/*
int CubicInflectionPoints(double *c, double *solutions){
double ax=c[2]-c[0];
double bx=c[4]-c[2]-ax;
double cx=c[6]-c[4]-ax-bx-bx;
double ay=c[3]-c[1];
double by=c[5]-c[3]-ay;
double cy=c[7]-c[5]-ay-by-by;
double eqn[3];
double curve[8];
double splitpoint=0.0;
double solutions[2];
int solcount=0;
int roots,i;
eqn[0]=ax*by-ay*bx;
eqn[1]=ax*cy-ay*cx;
eqn[2]=bx*cy-by*cx;
roots=SolveQuadratic(eqn,eqn);
for(i=0;i<roots;i++){
 if(eqn[i]>0 && eqn[i]<1){
  solutions[solcount++]=eqn[i];
 }
}
if(solutions[0]>solutions[1]){
 double tmp=solutions[0];
 solutions[0]=solutions[1];
 solutions[1]=tmp;
}
return solcount;
}
*/

---