curve.c

Go to the documentation of this file.

#include "curve.h"

double CurveGetX(CURVEDATA *curve, double t){
 switch(curve->u.type){
  case 0:
   return curve->u.pt.x0;
  case 1:
   return curve->u.line.x0+t*(curve->u.line.x1-curve->u.line.x0);
  case 2:{
   double x0=curve->u.quad.x0;
   double cx=curve->u.quad.x1;
   double x1=curve->u.quad.x2;
   x1=cx+(x1-cx)*t;
   x0=x0+(cx-x0)*t;
   cx=x0+(x1-x0)*t;
   return cx;
  }
  case 3:{
   double x0=curve->u.cubic.x0;
   double cx0=curve->u.cubic.x1;
   double cx1=curve->u.cubic.x2;
   double x1=curve->u.cubic.x3;
   x1=cx1+(x1-cx1)*t;
   x0=x0+(cx0-x0)*t;
   cx0=cx0+(cx1-cx0)*t;
   cx1=cx0+(x1-cx0)*t;
   cx0=x0+(cx0-x0)*t;
   return cx0+(cx1-cx0)*t;
  }
  default:return 0;
 }
}

double CurveGetY(CURVEDATA *curve, double t){
 switch(curve->u.type){
  case 0:
   return curve->u.pt.y0;
  case 1:
   return curve->u.line.y0+t*(curve->u.line.y1-curve->u.line.y0);
  case 2:{
   double y0=curve->u.quad.y0;
   double cy=curve->u.quad.y1;
   double y1=curve->u.quad.y2;
   y1=cy+(y1-cy)*t;
   y0=y0+(cy-y0)*t;
   cy=y0+(y1-y0)*t;
   return cy;
  }
  case 3:{
   double y0=curve->u.cubic.y0;
   double cy0=curve->u.cubic.y1;
   double cy1=curve->u.cubic.y2;
   double y1=curve->u.cubic.y3;
   y1=cy1+(y1-cy1)*t;
   y0=y0+(cy0-y0)*t;
   cy0=cy0+(cy1-cy0)*t;
   cy1=cy0+(y1-cy0)*t;
   cy0=y0+(cy0-y0)*t;
   return cy0+(cy1-cy0)*t;
  }
  default:return 0;
 }
}

double CurveGetX0(CURVEDATA *curve){
 switch(curve->u.type){
  case 0:
   return curve->u.pt.x0;
  case 1:
  case 2:
  case 3:
   return curve->u.line.x0;
  default:return 0;
 }
}

double CurveGetY0(CURVEDATA *curve){
 switch(curve->u.type){
  case 0:
   return curve->u.pt.y0;
  case 1:
   return curve->u.line.y0;
  case 2:
   return curve->u.line.y0;
  case 3:
   return curve->u.line.y0;
  default:return 0;
 }
}

double CurveGetX1(CURVEDATA *curve){
 switch(curve->u.type){
  case 0:
   return curve->u.pt.x0;
  case 1:
   return curve->u.line.x1;
  case 2:
   return curve->u.quad.x2;
  case 3:
   return curve->u.cubic.x3;
  default:return 0;
 }
}

double CurveGetY1(CURVEDATA *curve){
 switch(curve->u.type){
  case 0:
   return curve->u.pt.y0;
  case 1:
   return curve->u.line.y1;
  case 2:
   return curve->u.quad.y2;
  case 3:
   return curve->u.cubic.y3;
  default:return 0;
 }
}


void InitPointData(CURVEDATA *curve, double x, double y){
 LINEDATA *line=&curve->u.line;
 line->type=0;
 line->x0=x;
 line->y0=y;
 line->minX=line->maxX=x;
 line->minY=line->maxY=y;
}

void InitLineData(CURVEDATA *curve, double x0, double y0, double x1, double y1){
 LINEDATA *line=&curve->u.line;
 if(x0==x1 && y0==y1){
  InitPointData(curve,x0,y0);
 } else {
  line->type=1;
  line->x0=x0;
  line->x1=x1;
  line->y0=y0;
  line->y1=y1;
  line->minX=min(x0,x1);
  line->maxX=max(x0,x1);
  line->minY=min(y0,y1);
  line->maxY=max(y0,y1);
 }
}

void InitQuadData(CURVEDATA *curve, double *params){
 QUADDATA *line=&curve->u.quad;
 line->type=2;
 line->x0=params[0];
 line->y0=params[1];
 line->x1=params[2];
 line->y1=params[3];
 line->x2=params[4];
 line->y2=params[5];
 line->minX=min(line->x0,line->x1);
 line->minX=min(line->x2,line->minX);
 line->maxX=max(line->x0,line->x1);
 line->maxX=max(line->x2,line->maxX);
 line->minY=min(line->y0,line->y1);
 line->minY=min(line->y2,line->minY);
 line->maxY=max(line->y0,line->y1);
 line->maxY=max(line->y2,line->maxY);
}


void InitCubicData(CURVEDATA *curve, double *params){
 CUBICDATA *line=&curve->u.cubic;
  line->type=3;
  line->x0=params[0];
  line->y0=params[1];
  line->x1=params[2];
  line->y1=params[3];
  line->x2=params[4];
  line->y2=params[5];
  line->x3=params[6];
  line->y3=params[7];
  line->minX=min(line->x0,line->x1);
  line->minX=min(line->x2,line->minX);
  line->minX=min(line->x3,line->minX);
  line->maxX=max(line->x0,line->x1);
  line->maxX=max(line->x2,line->maxX);
  line->maxX=max(line->x3,line->maxX);
  line->minY=min(line->y0,line->y1);
  line->minY=min(line->y2,line->minY);
  line->minY=min(line->y3,line->minY);
  line->maxY=max(line->y0,line->y1);
  line->maxY=max(line->y2,line->maxY);
  line->maxY=max(line->y3,line->maxY);
}

void SplitCurve(
 CURVEDATA *curve,
 CURVEDATA *curveLeft,
 CURVEDATA *curveRight,
 double t
){
 switch(curve->u.type){
  case 0:{
   double x=curve->u.pt.x0;
   double y=curve->u.pt.y0;
   InitPointData(curveLeft,x,y);
   InitPointData(curveRight,x,y);
   break;
  }
  case 1:{
   double line1[4];
   double line2[4];
   LineSubdivide(&curve->u.line.x0,line1,line2,t);
   InitLineData(curveLeft,line1[0],line1[1],line1[2],line1[3]);
   InitLineData(curveRight,line2[0],line2[1],line2[2],line2[3]);
   break;
  }
  case 2:{
   double quad1[6];
   double quad2[6];
   QuadSubdivide(&curve->u.quad.x0,quad1,quad2,t);
   InitQuadData(curveLeft,quad1);
   InitQuadData(curveRight,quad2);
   break;
  }
  case 3:{
   double cubic1[8];
   double cubic2[8];
   CubicSubdivide(&curve->u.cubic.x0,cubic1,cubic2,t);
   InitCubicData(curveLeft,cubic1);
   InitCubicData(curveRight,cubic2);
   break;
  }
  default:{
   break;
  }
 }
}

CURVEDATA *ReduceCurve(
 CURVEDATA *curve,
 CURVEDATA *curveDst,
 double tMin,
 double tMax
){
 switch(curve->u.type){
  case 0:
   InitPointData(curveDst,curve->u.pt.x0,curve->u.pt.y0);
   return curve;
  case 1:{
   double dx0=curve->u.line.x0+tMin*(curve->u.line.x1-curve->u.line.x0);
   double dy0=curve->u.line.y0+tMin*(curve->u.line.y1-curve->u.line.y0);
   double dx1=curve->u.line.x0+tMax*(curve->u.line.x1-curve->u.line.x0);
   double dy1=curve->u.line.y0+tMax*(curve->u.line.y1-curve->u.line.y0);
   InitLineData(curveDst,dx0,dy0,dx1,dy1);
   return curveDst;
  }
  case 2:{
   double quad[8];
   ASSERT(tMin>=0 && tMin<=1.0);
   ASSERT(tMax>=0 && tMax<=1.0);
   QuadSubdivide(&curve->u.quad.x0,NULL,quad,tMin);
   tMax=(tMin==1.0) ? 0.0 : (tMax-tMin)/(1.0-tMin);
   QuadSubdivide(quad,quad,NULL,tMax);
   InitQuadData(curveDst,quad);
   return curveDst;
  }
  case 3:{
   double cubic[8];
   ASSERT(tMin>=0.0 && tMin<=1.0);
   ASSERT(tMax>=0.0 && tMax<=1.0);
   CubicSubdivide(&curve->u.cubic.x0,NULL,cubic,tMin);
   tMax=(tMin==1.0) ? 0.0 : (tMax-tMin)/(1.0-tMin);
   CubicSubdivide(cubic,cubic,NULL,tMax);
   InitCubicData(curveDst,cubic);
   return curveDst;
  }
  default:return curveDst;
 }
}

static int FindCubicZeros(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;
 int num=SolveQuadratic(zeros,zeros);
 int ret=0;
 for(int i=0;i<num;i++){
  double t=zeros[i];
  if(t>0&&t<1){
   zeros[ret]=t;
   ret++;
  }
 }
 return ret;
}

static int SortDoubles(const double *a, const double *b){
 if(a==b||*a==*b)return 0;
 if(*a<*b)return -1;
 return 1;
}

int FindInflectionPoints(CURVEDATA *curve, double *inflections){
 int count=0;
 switch(curve->u.type){
  default:
   return count;// Set to 0
  case 2:{
   double zero;
   double p0,p1,p2;
   p0=curve->u.quad.x0;
   p1=curve->u.quad.x1;
   p2=curve->u.quad.x2;
   // Will return infinity if division results in 0
   zero=-(p1+p1-p0-p0)/(2.0*(p0-p1-p1+p2));
   if(zero>0&&zero<1){
    inflections[count++]=zero;
   }
   p0=curve->u.quad.y0;
   p1=curve->u.quad.y1;
   p2=curve->u.quad.y2;
   // Will return infinity if division results in 0
   zero=-(p1+p1-p0-p0)/(2.0*(p0-p1-p1+p2));
   if(zero>0&&zero<1){
    inflections[count++]=zero;
   }
   if(count==2){
    double inf1=inflections[0];
    double inf2=inflections[1];
    if(inf1>inf2){
     inflections[0]=inf2;
     inflections[1]=inf1;
    }
   }
   return count;
  }
  case 3:{
   double zeros[3];
   double p0,p1,p2,p3;
   int num,i;
   p0=curve->u.cubic.x0;
   p1=curve->u.cubic.x1;
   p2=curve->u.cubic.x2;
   p3=curve->u.cubic.x3;
   num=FindCubicZeros(zeros,p0,p1,p2,p3);
   for(i=0;i<num;i++){
    inflections[count++]=zeros[i];
   }
   p0=curve->u.cubic.y0;
   p1=curve->u.cubic.y1;
   p2=curve->u.cubic.y2;
   p3=curve->u.cubic.y3;
   num=FindCubicZeros(zeros,p0,p1,p2,p3);
   for(i=0;i<num;i++){
    inflections[count++]=zeros[i];
   }
   qsort(inflections,count,sizeof(double),SortDoubles);
   return count;
  }
 }
}

int YDirection(CURVEDATA *curve, BOOL yaxis, double y, double t){
 switch(curve->u.type){
  case 0:
   return 0;
  case 1:{
   LINEDATA *line=&curve->u.line;
   double y0=(yaxis) ? line->y0 : line->x0;
   double y1=(yaxis) ? line->y1 : line->x1;
   return (y1<y0) ? -1 : 1;
  }
  case 2:{
   QUADDATA *quad=&curve->u.quad;
   double y0=(yaxis) ? quad->y0 : quad->x0;
   double y2=(yaxis) ? quad->y2 : quad->x2;
   return (y2<y0) ? -1 : 1;
  }
  case 3:{
   CUBICDATA *cubic=&curve->u.cubic;
   double y0=(yaxis) ? cubic->y0 : cubic->x0;
   double y3=(yaxis) ? cubic->y3 : cubic->x3;
   return (y3<y0) ? -1 : 1;
  }
  default: return 0;
 }
}

int SolveCurve(CURVEDATA *curve, BOOL yaxis, double y, double *r){
 switch(curve->u.type){
  case 0:
   return 0;
  case 1:{
   LINEDATA *line=&curve->u.line;
   double eqn[2];
   double y0=(yaxis) ? line->y0 : line->x0;
   double y1=(yaxis) ? line->y1 : line->x1;
   eqn[1] = y1 - y0;
   eqn[0] = y0 - y;
   return SolveLinear(eqn,r);
  }
  case 2:{
   QUADDATA *quad=&curve->u.quad;
   double eqn[3];
   double y0=(yaxis) ? quad->y0 : quad->x0;
   double cy=(yaxis) ? quad->y1 : quad->x1;
   double y1=(yaxis) ? quad->y2 : quad->x2;
   eqn[2] = y1 + y0 - cy - cy;
   eqn[1] = cy + cy - y0 - y0;
   eqn[0] = y0 - y;
   return SolveQuadratic(eqn,r);
  }
  case 3:{
   CUBICDATA *cubic=&curve->u.cubic;
   double eqn[4];
   double y0=(yaxis) ? cubic->y0 : cubic->x0;
   double cy1=(yaxis) ? cubic->y1 : cubic->x1;
   double cy2=(yaxis) ? cubic->y2 : cubic->x2;
   double y1=(yaxis) ? cubic->y3 : cubic->x3;
   eqn[0] = y0 - y;
   eqn[1] = 3 * (cy1 - y0);
   eqn[2] = 3 * (cy2 + y0 - cy1 - cy1);
   eqn[3] = y1 + (cy1 - cy2) * 3.0 - y0;
   return SolveCubic(eqn,r);
  }
  default: return 0;
 }
}

BOOL IsHorizontalLine(CURVEDATA * curve){
 return (curve->u.c.minY==curve->u.c.maxY);
}

void DebugOutCurve(CURVEDATA *curve){
#ifdef DEBUG
 switch(curve->u.type){
  case 0:
   DebugOut("Curve[0]=%f,%f",
     (float)curve->u.line.x0,
     (float)curve->u.line.y0);
  case 1:
   DebugOut("Curve[1]=%f,%f,%f,%f",
     (float)curve->u.line.x0,
     (float)curve->u.line.y0,
     (float)curve->u.line.x1,
     (float)curve->u.line.y1);
  case 2:
   DebugOut("Curve[2]=%f,%f,%f,%f,%f,%f",
     (float)curve->u.quad.x0,
     (float)curve->u.quad.y0,
     (float)curve->u.quad.x1,
     (float)curve->u.quad.y1,
     (float)curve->u.quad.x2,
     (float)curve->u.quad.y2);
  case 3:
   DebugOut("Curve[3]=%f,%f,%f,%f,%f,%f,%f,%f",
     (float)curve->u.cubic.x0,
     (float)curve->u.cubic.y0,
     (float)curve->u.cubic.x1,
     (float)curve->u.cubic.y1,
     (float)curve->u.cubic.x2,
     (float)curve->u.cubic.y2,
     (float)curve->u.cubic.x3,
     (float)curve->u.cubic.y3);
 }
#endif
}

double CurveLength(CURVEDATA *curve){
 if(!curve)return 0;
 switch(curve->u.type){
  case 1:
   return PointDistance(
       curve->u.line.x0,curve->u.line.y0,
       curve->u.line.x1,curve->u.line.y1);
  case 2:
   return QuadCalcLength(&curve->u.quad.x0,1e-5);
  case 3:
   return CubicCalcLength(&curve->u.cubic.x0,1e-5);
  default:return 0;
 }
}

BOOL CurvesEqual(CURVEDATA *a, CURVEDATA *b){
 if(a==b)return TRUE;
 if(a->u.type!=b->u.type)return FALSE;
 switch(a->u.type){
  case 0:
   return (a->u.pt.x0==b->u.pt.x0 &&
           a->u.pt.y0==b->u.pt.y0);
  case 1:
   return (a->u.line.x0==b->u.line.x0 &&
           a->u.line.x1==b->u.line.x1 &&
           a->u.line.y0==b->u.line.y0 &&
           a->u.line.y1==b->u.line.y1
          );
  case 2:
   return (a->u.quad.x0==b->u.quad.x0 &&
           a->u.quad.y0==b->u.quad.y0 &&
           a->u.quad.x1==b->u.quad.x1 &&
           a->u.quad.y1==b->u.quad.y1 &&
           a->u.quad.x2==b->u.quad.x2 &&
           a->u.quad.y2==b->u.quad.y2
          );
  case 3:
   return (a->u.cubic.x0==b->u.cubic.x0 &&
           a->u.cubic.y0==b->u.cubic.y0 &&
           a->u.cubic.x1==b->u.cubic.x1 &&
           a->u.cubic.y1==b->u.cubic.y1 &&
           a->u.cubic.x2==b->u.cubic.x2 &&
           a->u.cubic.y2==b->u.cubic.y2 &&
           a->u.cubic.x3==b->u.cubic.x3 &&
           a->u.cubic.y3==b->u.cubic.y3
           );
  default:return TRUE;
 }
}

---