Curve representation by ICAS (Inner Centered Arcs)

#include "IcasBase.h"
#include <math.h>
//#include <stdio.h> //NULL
#ifndef TRACE_ON
#define NULL (void*)0
#endif

static double    g_t    = 0.0;
static int       g_reso = 0;//density of curve points
static int       g_deg  = 0;
static int       g_asiz = 0;
static double    g_grid_dist = 5.0;
static ICASPOINT g_a[ L_ICAS_MAX ];//working area
ICASPOINT        g_C[ L_ICAS_MAX ];//porting area

//Get Inner Center
int GetInC(ICASPOINT *S,  ICASPOINT *E,  ICASPOINT *C, 
		   ICASPOINT *SS, ICASPOINT *EE, ICASPOINT *I,
		   double *ds, double *de)
{
	double s,e,c;//length of edges
	double x1,y1,x2,y2,x3,y3,t,ix,iy;
	double sx,sy,ex,ey;

	//OutputDebugString("OK");//Windows.h

	x1 = (E->x - S->x);
	y1 = (E->y - S->y);
	c = sqrt(x1*x1 + y1*y1);//SE
	x2 = (C->x - S->x);
	y2 = (C->y - S->y);
	e = sqrt(x2*x2 + y2*y2);//SC
	x3 = (C->x - E->x);
	y3 = (C->y - E->y);
	s = sqrt(x3*x3 + y3*y3);//CE
	t = s + e + c;//surround length of triangle
	if (t <= EPS) {
		return -1;
	}
	ix = ((s*S->x+e*E->x+c*C->x)/t);
	iy = ((s*S->y+e*E->y+c*C->y)/t);
	I->x = ix;//Inner Center
	I->y = iy;

	s = (x1 *(iy-S->y) - y1 *(ix-S->x))/(x1 * y2 - x2 * y1);
	sx = S->x + s * x2;
	sy = S->y + s * y2;
	SS->x = sx;//Next Control point (front)
	SS->y = sy;

	e = (x1 *(iy-E->y) - y1 *(ix-E->x))/(x1 * y3 - x3 * y1);
	ex = E->x + e * x3;
	ey = E->y + e * y3;
	EE->x = ex;//Next Control point (rear)
	EE->y = ey;

	//length of IS and IE 
	if (NULL != ds ) {
		*ds = sqrt((ix-S->x)*(ix-S->x)+(iy-S->y)*(iy-S->y));
	}
	if (NULL != de ) {
		*de = sqrt((ix-E->x)*(ix-E->x)+(iy-E->y)*(iy-E->y));
	}
	return 0;
}
//Get Edge Center
int GetMidC(ICASPOINT *S,  ICASPOINT *E,  ICASPOINT *C, 
		    ICASPOINT *SS, ICASPOINT *EE, ICASPOINT *I,
		   double *ds, double *de)
{
	double s;//length of edges
	double Mx,My,x1,y1;

	Mx = (E->x + S->x)*0.5;
	My = (E->y + S->y)*0.5;
	I->x = Mx;//Mid Center
	I->y = My;
	x1 = (C->x - E->x)*0.5;
	y1 = (C->y - E->y)*0.5;
	SS->x = S->x - x1;//Next Control point (front)
	SS->y = S->y - y1;
	EE->x = E->x + x1;//Next Control point (rear)
	EE->y = E->y + y1;

	if ((ds != NULL) || (de != NULL)) {
		s = sqrt((Mx-S->x)*(Mx-S->x)+(My-S->y)*(My-S->y));
		//length of MS and ME
		if (NULL != ds ) {
			*ds = s;
		}
		if (NULL != de ) {
			*de = s;
		}
	}
	return 0;
}
//Biarc initiator
int BrcInC(ICASPOINT *S,  ICASPOINT *C,  ICASPOINT *E,
		   ICASPOINT *SS, ICASPOINT *EE, ICASPOINT *I,
		   int Typ, double *ds, double *de)
{
	int sts = 0;
	switch(Typ & L_SEG_TANGENT)
	{
	case L_LEFT_TAN:
		if (0 == g_deg) {
			if (GetInC(S, C, E, SS, EE, I, ds, de)<0) return -2;
			sts = L_LEFT_TAN;
		}
		break;
	case L_LEFT_TAN2:
		if (1 == g_deg) {
			if (GetMidC(S, E, C, SS, EE, I, ds, de)<0) return -3;
			sts = L_LEFT_TAN2;
		}
		break;
	case L_RIGHT_TAN:
		if (0 == g_deg) {
			if (GetInC(C, E, S, SS, EE, I, ds, de)<0) return -2;
			sts = L_RIGHT_TAN;
		}
		break;
	case L_RIGHT_TAN2:
		if (1 == g_deg) {
			if (GetMidC(E, S, C, EE, SS, I, de, ds)<0) return -3;
			sts = L_RIGHT_TAN2;
		}
		break;
	default://not tangent to data segment
		break;
	}
	if (sts) return 0;//

	switch(Typ & L_CURV_TANGENT)
	{
	default://
	case L_NORMAL:
	case L_BRIDGE:
		if (GetInC(S, E, C, SS, EE, I, ds, de)<0) return -4;
		break;
	case L_LINE:
	case L_IDENTDIR:
		I->x = S->x + (E->x - S->x)*g_t;
		I->y = S->y + (E->y - S->y)*g_t;
		return 1;//break;
	}
	return 0;
}

ICASPOINT * GetParametricPoint(ICASTRI *T, double tt)
{
	static ICASPOINT P;
	double s = 0.5;//from here to tt
	//Set Global parameters (g_deg, g_t)
	g_deg = 0;
	g_t   = tt;
	//get the point parametricaly, corresponding g_t [0.0-1.0)
	GetPointR(&(T->P[0]), &(T->P[1]), &(T->P[2]), T->Typ, &s, &P);

	return &P;
}

//get parametric curve point
//'R'ecurrsive
int GetPointR(ICASPOINT *S, ICASPOINT *C, ICASPOINT *E,
			  int Typ, double *s, ICASPOINT *P )
{
	static double ss = 0.5;//difference
	ICASPOINT SS,EE,I;
	//Get Inner Center
	if (1 == BrcInC(S, C, E, &SS, &EE, &I, Typ, NULL, NULL)) {
		*P = I;  //L_LINE | L_IDENTDIR
		return 0;
	}
	//match g_t?
	if (fabs(g_t - *s) < EPS) {
		*P = I;
		//TRACE("R(%d)\n",g_deg);
	}
	//go reccursive
	else {
		g_deg++;
		ss *= 0.5;
		//
		if (*s > g_t) {
			*s = *s - ss;
			GetPointR(S, &SS, &I, Typ & L_LEFTMASK,  s, P);
		}
		else {//if (g_t > *s) {
			*s = *s + ss;
			GetPointR(&I, &EE, E, Typ & L_RIGHTMASK, s, P);
		}
		//
		g_deg--;
		ss *= 2.0;
	}
	return 0;
}
//interpolate by precision
//'R'ecurrsive
int GetArrayR(ICASPOINT *S, ICASPOINT *C, ICASPOINT *E, int Typ )
{
	ICASPOINT SS,EE,I;
	double ds, de;
	//Get Inner Center
	if (1 == BrcInC(S, C, E, &SS, &EE, &I, Typ, &ds, &de)) return 0;//L_LINE

	//enough precision?
	if (g_grid_dist >= ds && g_grid_dist >= de) {
		if (g_asiz >= L_ICAS_MAX) return 0;
		g_a[ g_asiz ] = I;
		++g_asiz;
		//TRACE("ArrayR[%d]\n",g_asiz);
	}
	//go reccursive
	else {
		//g_deg++;

		//Left first
		if (g_grid_dist < ds)
		GetArrayR(S, &SS, &I, Typ & L_LEFTMASK);

		//Right last
		if (g_grid_dist < de)
		GetArrayR(&I, &EE, E, Typ & L_RIGHTMASK);

		//g_deg--;
	}
	return 0;
}
//interpolate by tesselate degree
//'R'ecurrsive 2
int GetArrayR2(ICASPOINT *S, ICASPOINT *C, ICASPOINT *E, int Typ, int dig )
{
	ICASPOINT SS,EE,I;
	//enough tesselation?
	if (dig >= g_reso) {
		return 0;
	}
	//Get Inner Center
	if (1 == BrcInC(S, C, E, &SS, &EE, &I, Typ, NULL, NULL)) return 0;//L_LINE

	//go reccursive
	g_deg++;

	//Left first
	GetArrayR2(S, &SS, &I, Typ & L_LEFTMASK, dig+1);

	if (g_asiz >= L_ICAS_MAX) {
		g_deg--;
		return 0;
	}
	g_a[ g_asiz ] = I;
	++g_asiz;
	//TRACE("ArrayR2<%d>\n",g_asiz);

	//Right last
	GetArrayR2(&I, &EE, E, Typ & L_RIGHTMASK, dig+1);

	g_deg--;
	return 0;
}

//Inner Center-ed Curve inside Triangle
int IcasCurve(ICASTRI *Tri, int CurvTyp, int Reso, double dReso)
{
	int i,Csiz = 0;

	if (CurvTyp & L_CLOSE) {//CurvTyp & 0x1 : close
	}
	else {//CurvTyp==0 : open
		if (Tri->Typ & L_BRIDGE) return 0;//elliminate BRIDGE part
	}
	//
	if (Tri->Typ & L_IDENTDIR) Reso = 0;//to be flat line
	//
	g_asiz = 0;//max L_ICAS_MAX
	g_deg  = 0;//tesselate degree
	//apply resolution
	switch( Reso )
	{
	case 0://line segment
	case 1://single inner center
	case 2://three inner centers
	case 3://seven inner centers
	case 4://if you want more tesselation, set bigger 'g_reso' and call GetArrayR2().
		g_reso = Reso;
		GetArrayR2(&(Tri->P[0]), &(Tri->P[1]), &(Tri->P[2]), Tri->Typ, 0);//tesselate degree
		break;
	case 5://if you want more resolution, set smaller 'g_grid_dist' and call GetArrayR().
		g_grid_dist = dReso;
		GetArrayR(&(Tri->P[0]), &(Tri->P[1]), &(Tri->P[2]), Tri->Typ);//resolution 
		break;
	}
	//copy to porting area
	g_C[0]        = Tri->P[0];//first
	for(i=0; i<g_asiz; ++i) {
		g_C[i+1]  = g_a[i];
	}
	g_C[g_asiz+1] = Tri->P[2];//last
	Csiz = g_asiz+2;
	//TRACE("Csiz=%d\n",Csiz);
	return Csiz;
}
//move Control point = Q
int GetSideTri(ICASTRI *Tri, ICASPOINT *Q, int Idx)
{
	int a = 0, b = 2;
	double P1ax,P1ay,Pbax,Pbay,PbQx,PbQy,s,d;
	if (Tri->Typ & L_IDENTDIR) return 0;
	if (Idx!=0) {
		a = 2; b = 0;
	}
	P1ax = Tri->P[1].x - Tri->P[a].x;
	P1ay = Tri->P[1].y - Tri->P[a].y;
	Pbax = Tri->P[b].x - Tri->P[a].x;
	Pbay = Tri->P[b].y - Tri->P[a].y;
	PbQx = Tri->P[b].x - Q->x;
	PbQy = Tri->P[b].y - Q->y;
	d = P1ax * PbQy - P1ay * PbQx;
	if (fabs(d) > EPS) {
		s = (PbQy * Pbax - PbQx * Pbay)/d;
		Tri->P[1].x = Tri->P[a].x + s * P1ax;
		Tri->P[1].y = Tri->P[a].y + s * P1ay;
	}
	return 0;
}
//move Data point = Q
int GetSideTri2(ICASTRI *Tri, ICASPOINT *Q, int Idx)
{
	int a = 0, b = 2;
	double P1ax,P1ay,P1bx,P1by,PaQx,PaQy,s,d;
	if (Tri->Typ & L_IDENTDIR) return 0;
	if (Idx!=0) {
		a = 2; b = 0;
	}
	P1ax = Tri->P[1].x - Tri->P[a].x;
	P1ay = Tri->P[1].y - Tri->P[a].y;
	P1bx = Tri->P[1].x - Tri->P[b].x;
	P1by = Tri->P[1].y - Tri->P[b].y;
	PaQx = Tri->P[a].x - Q->x;
	PaQy = Tri->P[a].y - Q->y;
	d = P1ax * P1by - P1ay * P1bx;
	if (fabs(d) > EPS) {
		s = (P1bx * PaQy - PaQx * P1by)/d;
		Tri->P[1].x = Tri->P[a].x + s * P1ax;
		Tri->P[1].y = Tri->P[a].y + s * P1ay;
		Tri->P[b]   = *Q;
	}
	return 0;
}
//Terminal control point //isosceles triangle 
int GetCntlTerm(ICASPOINT *S, ICASPOINT *E, ICASPOINT *EE, ICASPOINT *P )
{
	double Sx,Sy,Ux,Uy,Vx,Vy,s,t,u,v;
	Ux = (E->x - S->x);
	Uy = (E->y - S->y);
	Vx = (EE->x - E->x);
	Vy = (EE->y - E->y);
	u  = sqrt(Ux*Ux + Uy*Uy);
	v  = sqrt(Vx*Vx + Vy*Vy);

	if ( u < EPS || v < EPS ) return -1;
	//unit vector
	Ux = Ux / u;
	Uy = Uy / u;
	Vx = Vx / v;
	Vy = Vy / v;
	//vector of control line
	Vx = Ux + Vx;
	Vy = Uy + Vy;
	//perpendicular vector to SE
	Ux = E->y - S->y;
	Uy = S->x - E->x;
	//constant vector
	Sx = (S->x - E->x)*0.5;
	Sy = (S->y - E->y)*0.5;

	t = Ux * Vy - Vx * Uy;
	if ( fabs(t) < EPS ) return -2;
	s = (Ux * Sy - Uy * Sx)/t;

	//Terminal control point
	P->x = E->x + s*Vx;
	P->y = E->y + s*Vy;
	return 0;
}
//
int GetCross(ICASPOINT *D1, ICASPOINT *D2,
			 ICASPOINT *C1, ICASPOINT *C2, ICASPOINT *P)
{
	double D21x,D21y,C21x,C21y,DCx,DCy,t,d;
	D21x = D1->x - D2->x;
	D21y = D1->y - D2->y;
	C21x = C1->x - C2->x;
	C21y = C1->y - C2->y;
	DCx  = D2->x - C2->x;
	DCy  = D2->y - C2->y;
	d = D21x * C21y - C21x * D21y;
	if (fabs(d) > EPS) {
		t = (C21x * DCy - C21y * DCx)/d;
		P->x = D2->x + t*D21x;
		P->y = D2->y + t*D21y;
	}
	return 0;
}
//Hit test
BOOL CSClip(double *sx, double *sy,double *x, double *y, 
		    double *l,  double *r, double *t, double *b )
{
	double s,x0,x1,y0,y1,xd,yd;
	//Cohen-Sutherland Clipping algorithm
	if      (*x < *l && *sx < *l) return FALSE;//B0100
	else if (*x > *r && *sx > *r) return FALSE;//B0001
	if      (*y > *t && *sy > *t) return FALSE;//B1000
	else if (*y < *b && *sy < *b) return FALSE;//B0010
	//
	x0 = *l - *sx;
	x1 = *r - *sx;
	y0 = *b - *sy;
	y1 = *t - *sy;
	xd = *x - *sx;
	yd = *y - *sy;
	//
	s = *sx + y0 * xd / yd;
	if ( *l < s && s < *r ) {
		return TRUE;
	}
	s = *sx + y1 * xd / yd;
	if ( *l < s && s < *r ) {
		return TRUE;
	}

	s = *sy + x0 * yd / xd;
	if ( *b < s && s < *t ) {
		return TRUE;
	}
	s = *sy + x1 * yd / xd;
	if ( *b < s && s < *t ) {
		return TRUE;
	}
	return FALSE;
}
BOOL TriInsideChk(ICASPOINT P[], ICASPOINT *Q)
{
	double P1x,P1y,P10x,P10y,P12x,P12y,QPx,QPy,u,v,d;
	P1x  = P[1].x ;
	P1y  = P[1].y ;
	P10x = P[0].x - P1x;
	P10y = P[0].y - P1y;
	P12x = P[2].x - P1x;
	P12y = P[2].y - P1y;

	d = P10x * P12y - P10y * P12x;
	if (fabs(d) < EPS) return FALSE;
	QPx  = Q->x - P1x;
	QPy  = Q->y - P1y;

	u = (P12y * QPx - P12x * QPy) / d;
	if (0.0 < u && u < 1.0)
	{
		v = (P10x * QPy - P10y * QPx) / d;
		if (0.0 < v && v < 1.0)
		{
			if ( u+v < 1.0)
			{
				return TRUE;
			}
		}
	}
	return FALSE;
}
//Get ControlPoints
int GetCntl3(ICASPOINT *S, ICASPOINT *C, ICASPOINT *E,
		     ICASPOINT *C1,ICASPOINT *C2,ICASPOINT *C3 )
{
	double Ux,Uy;
	Ux = (C->x - E->x);
	Uy = (C->y - E->y);
	if (Ux*Ux + Uy*Uy < EPS) return -1;
	C1->x = S->x + Ux;
	C1->y = S->y + Uy;
	C3->x = S->x - Ux;
	C3->y = S->y - Uy;

	Ux = (E->x - S->x);
	Uy = (E->y - S->y);
	if (Ux*Ux + Uy*Uy < EPS) return -2;
	C2->x = C->x + Ux;
	C2->y = C->y + Uy;
	return 0;
}
void Point2Work(ICASWORK *pW, ICASPOINT *pP)
{
	pW->typ = pP->typ;
	pW->x   = pP->x;
	pW->y   = pP->y;
	pW->l   = 0.0;
	pW->ex  = 0.0;
	pW->ey  = 0.0;
	pW->vx  = 0.0;
	pW->vy  = 0.0;
	pW->sts = 0;//not ready
}
//complete E members from E and S values
int Work2Vector(ICASWORK *pS, ICASWORK *pE)
{
	if (0 == pE->sts) {//not ready
		pE->vx = pE->x - pS->x;
		pE->vy = pE->y - pS->y;
		pE->l  = sqrt((pE->vx * pE->vx) + (pE->vy * pE->vy));
		pE->ex = pE->vx / pE->l;//unit vector x
		pE->ey = pE->vy / pE->l;//unit vector y
		pE->sts = 1;//Done
	}
	return pE->sts;
}
int ChkCntl(ICASWORK *SS, ICASWORK *S,
			ICASWORK *E,  ICASWORK *EE,
		    ICASPOINT *C1, ICASPOINT *M )
{
	if (( S->typ & L_MULTIPLE)&&( E->typ & L_MULTIPLE)) {
		if (fabs(S->x - E->x) < EPS && fabs(S->y - E->y) < EPS) {//S=E  prevent store
			//TRACE("ChkCntl:S==E\n");
			return -3;//S=E  prevent to store
		}
	}
	//SE middle point
	M->x = (S->x + E->x)*0.5;
	M->y = (S->y + E->y)*0.5;

	if ((SS->typ & L_MULTIPLE)&&( S->typ & L_MULTIPLE)) {//SS=S
		if (fabs(SS->x - S->x) < EPS && fabs(SS->y - S->y) < EPS) {
			if (( E->typ & L_MULTIPLE)&&(EE->typ & L_MULTIPLE)) {//E=EE
				if (fabs(E->x - EE->x) < EPS && fabs(E->y - EE->y) < EPS) {
					//TRACE("ChkCntl:SS==S && E==EE\n");
					*C1 = *M;//C1 = SE middle point
					return -4;//L_MLTSTART | L_MLTEND
				}
			}
			else {
				//TRACE("ChkCntl:SS==S\n");
				//GetCntlTerm(S,E,EE,C1);
				GetCntlTerm((ICASPOINT*)S,(ICASPOINT*)E,(ICASPOINT*)EE,C1);
				return -5;//L_MLTSTART
			}
		}
	}
	else if (( E->typ & L_MULTIPLE)&&(EE->typ & L_MULTIPLE)) {//E=EE
		if (fabs(E->x - EE->x) < EPS && fabs(E->y - EE->y) < EPS) {
			//TRACE("ChkCntl:E=EE\n");
			//GetCntlTerm(E,S,SS,C1);
			GetCntlTerm((ICASPOINT*)E,(ICASPOINT*)S,(ICASPOINT*)SS,C1);
			return -6;//L_MLTEND
		}
	}
	return 0;
}
//calculate control point
int GetCntl(ICASWORK *SS, ICASWORK *S,
			ICASWORK *E,  ICASWORK *EE,
		    ICASPOINT *C1, ICASPOINT *C2, ICASPOINT *M )
{
	double Vx,Vy,Sx,Sy,Ex,Ey,s,t,u;

	int sts = ChkCntl(SS, S, E, EE, C1, M );//set M
	if (0 > sts) return sts;//MULTIPLE check valid

	Work2Vector(SS,S);//A
	Work2Vector(S,E); //B
	Work2Vector(E,EE);//C

	//vector product
	//S-SS(A) vs S-E(B)
	s = S->vx * E->vy  - S->vy * E->vx;
	if (fabs(s) < EPS) {//M = SE middle point
		return -1;//S-SS and S-E are the identical line
	}
	//S-E(B)  vs E-EE(C)
	t = E->vx * EE->vy - E->vy * EE->vx;
	if (fabs(t) < EPS) {//M = SE middle point
		return -2;//S-E and E-EE are the identical line
	}
	//E-EE(C) vs S-SS(A)
	u = S->vx * EE->vy - S->vy * EE->vx;

	Sx = S->ex + E->ex; //A.ex + B.ex;//unit vector
	Sy = S->ey + E->ey; //A.ey + B.ey;
	Ex = E->ex + EE->ex;//B.ex + C.ex;
	Ey = E->ey + EE->ey;//B.ey + C.ey;
	Vx = E->vx;
	Vy = E->vy;

	if (L_PARALLEL == SS->typ) {//point SS on the identical line
		Sx = S->ex; //A.ex;
		Sy = S->ey; //A.ey;
	}
	if (L_PARALLEL == EE->typ) {//point EE on the identical line
		Ex = EE->ex;//C.ex;
		Ey = EE->ey;//C.ey;
	}
	//identical turnnig direction
	if (s*t > 0)
	{
		if (fabs(u) < EPS && (L_PARALLEL == SS->typ && L_PARALLEL == EE->typ)) {
			double l = E->l * 0.5;
			C1->x = S->x + S->ex * l; //A.ex * l;
			C1->y = S->y + S->ey * l; //A.ey * l;
			C2->x = E->x - EE->ex * l;//C.ex * l;
			C2->y = E->y - EE->ey * l;//C.ey * l;
			//boxy
			M->x = (C1->x + C2->x)*0.5;//C1C2 middle point
			M->y = (C1->y + C2->y)*0.5;
			return 2;
		}
		s  = (Ey*Vx - Ex*Vy)/(Sx*Ey - Ex*Sy);
		C1->x = S->x + (Sx*s);
		C1->y = S->y + (Sy*s);
		return 1;
	}
	//vector SS-S || vector E-EE (|| = parallel)
	else if (fabs(u) < EPS) {
		double l = E->l * 0.5;
		C1->x = S->x + S->ex * l; //Ax * l;
		C1->y = S->y + S->ey * l; //Ay * l;
		C2->x = E->x - EE->ex * l;//Cx * l;
		C2->y = E->y - EE->ey * l;//Cy * l;
		//Z fold//M = SE middle point
		return 2;
	}
	//Z fold
	else {
		double Ax,Ay,Cx,Cy;
		double C12x,C12y,w,v;
		double l = E->l * 0.5;
		w = sqrt(Sx*Sx + Sy*Sy);
		v = sqrt(Ex*Ex + Ey*Ey);
		Ax = Sx / w;
		Ay = Sy / w;
		Cx = Ex / v;
		Cy = Ey / v;
		C1->x = S->x + Ax * l;
		C1->y = S->y + Ay * l;
		C2->x = E->x - Cx * l;
		C2->y = E->y - Cy * l;
		C12x = C2->x - C1->x;
		C12y = C2->y - C1->y;
		u = Vx * C12y - Vy * C12x;

		if (fabs(u) < EPS) {//fail safe//M = SE middle point
		}
		//cross point SE and C1C2
		else {
			s = ((S->y - C1->y)*C12x - C12y*(S->x - C1->x))/u;
			M->x = S->x + (Vx*s);
			M->y = S->y + (Vy*s);
		}
		return 2;
	}
}
int ChkParallel(ICASPOINT Seg[], int SegCnt)
{
	int i, cnt = 0;
	int siz = SegCnt;
	int typ = L_NORMAL;
	double s,t;
	ICASPOINT U,V;
	if (siz < 3) return 0;
	--siz;//last -1
	U.x = Seg[0].x - Seg[siz].x;
	U.y = Seg[0].y - Seg[siz].y;
	for(i=0; i<siz; ++i)
	{
		V.x = Seg[i+1].x - Seg[i].x;
		V.y = Seg[i+1].y - Seg[i].y;
		s = V.x * V.x + V.y * V.y;//square length
		t = U.x * V.y - U.y * V.x;//vector product
		//identical check
		if (s < EPS) {//fabs(s) < EPS) {
			typ = L_MULTIPLE;
			Seg[i].typ = L_MULTIPLE;//same point also next too
			//TRACE("ChkParallel[%d]:MULTIPLE\n",i);
		}
		//parallel check
		else if (typ != L_MULTIPLE && fabs(t) < EPS) {
			typ = L_NORMAL;
			Seg[i].typ = L_PARALLEL;//on the unique line
			++cnt;
			//TRACE("ChkParallel[%d]:PARALLEL\n",i);
		}
		else {
			//TRACE("ChkParallel[%d]:%d\n",i,typ);
			Seg[i].typ = typ;//for previous L_MULTIPLE
			typ = L_NORMAL;
		}
		//
		U = V;
	}
	V.x = Seg[siz].x - Seg[0].x;
	V.y = Seg[siz].y - Seg[0].y;
	s = V.x * V.x + V.y * V.y;
	t = U.x * V.y - U.y * V.x;
	if (s < EPS) {//fabs(s) < EPS) {
		Seg[siz].typ = L_MULTIPLE;//same point 
		//TRACE("ChkParallelEND[%d]:MULTIPLE\n",siz);
	}
	else if (typ != L_MULTIPLE && fabs(t) < EPS) {
		Seg[siz].typ = L_PARALLEL;//on the line between Pre/Post 
		++cnt;
		//TRACE("ChkParallelEND[%d]:PARALLEL\n",siz);
	}
	else {
		//TRACE("ChkParallelEND[%d]:%d\n",siz,typ);
		Seg[siz].typ = typ;
	}
	return cnt;
}
int get_quad( ICASPOINT *O, ICASPOINT *P )//O:origin P:segment data
{
	if ( P->x > O->x + EPS ) {
		if ( P->y  < O->y - EPS ) return( 4 );
		else                      return( 1 );
	}
	else if ( P->x < O->x - EPS ) {
		if ( P->y  > O->y + EPS ) return( 2 );
		else                      return( 3 );
	}
	else {
		if (P->y  > O->y )        return( 2 );
		else                      return( 4 );
	}
}
//
BOOL ChkArea(ICASPOINT Seg[], int SegCnt, ICASPOINT *O)
{
    int i, j, q_ang, vs_quad, ve_quad, q_dif;
    double x_mid;
	int cnt = SegCnt;

	q_ang = 0;
	//Which quad is Seg[LAST] in?
	j = cnt-1;
	vs_quad = get_quad(O,&Seg[j]);

	for(i = 0; i < cnt; ++i)
	{
		//vertex on vertex
		if (fabs(Seg[i].x - O->x)<EPS && fabs(Seg[i].y - O->y)<EPS) {
		    return TRUE;//V_ON_V;
		}
		//Which quad is m_Seg[i] in?
		ve_quad = get_quad(O,&Seg[i]);

		//get quad transit
		q_dif = ve_quad - vs_quad;
		switch(q_dif) {
		case 2: case -2:
			if (fabs(Seg[j].y - Seg[i].y) < EPS) {
				x_mid = O->x;
			}
			else {
				x_mid = Seg[j].x - (((Seg[j].y-O->y)*(Seg[i].x-Seg[j].x))/(Seg[i].y-Seg[j].y));
			}
			if ( x_mid > O->x + EPS ) {
				q_dif *= -1;
			}
			else if ( x_mid < O->x - EPS ) {
			}
			else {//vertex on edge
				return TRUE;//V_ON_E;
			}
			break;
		case 3:
			q_dif = -1; break;
		case -3:
			q_dif = 1; break;
		}
		q_ang += q_dif;
		vs_quad = ve_quad;
		j = i;
	}
	q_ang = abs(q_ang);
	if (q_ang > 0 && 0 == q_ang%4) {//IN
		return TRUE;//V_IN;
	}
	else {//OUT
		return FALSE;//V_OUT;
	}
}

//add new Triangle
int SetInC(ICASTRI Tri[], int TriCnt, 
		   ICASPOINT *S, ICASPOINT *E, ICASPOINT *C,
		   int IcasTyp)
{
	if (NULL == Tri || TriCnt >= L_ICAS_MAX) return 0;
	Tri[TriCnt].Typ = IcasTyp;
	Tri[TriCnt].P[0] = *S;
	Tri[TriCnt].P[0].typ = L_NORMAL;
	Tri[TriCnt].P[1] = *C;
	Tri[TriCnt].P[1].typ = L_NORMAL;
	Tri[TriCnt].P[2] = *E;
	Tri[TriCnt].P[2].typ = L_NORMAL;
	++TriCnt;
	return TriCnt;
}
//dispatch by number of Control points
int SetTriangle(ICASTRI Tri[], int TriCnt,
				int cnt, int addflg,
				ICASPOINT *S,  ICASPOINT *E,
				ICASPOINT *C1, ICASPOINT *C2, ICASPOINT *M )
{
	if (cnt == 1) {
		TriCnt = SetInC(Tri, TriCnt, S, E, C1, (L_NORMAL | addflg));
	}
	else if (cnt == 2) {
		TriCnt = SetInC(Tri, TriCnt, S, M, C1, (L_NORMAL | L_AUTO_R | addflg));
		TriCnt = SetInC(Tri, TriCnt, M, E, C2, (L_NORMAL | L_AUTO_L | addflg));
	}
	//identical line
	else if (cnt == -1 || cnt == -2) {
		TriCnt = SetInC(Tri, TriCnt, S, E, M, (L_IDENTDIR | addflg));
	}
	//multiple points(start & end)
	else if (cnt == -4) {
		TriCnt = SetInC(Tri, TriCnt, S, E, C1, (L_MLTSTART | L_MLTEND | L_LINE | addflg));
	}
	//multiple point(start)
	else if (cnt == -5) {
		TriCnt = SetInC(Tri, TriCnt, S, E, C1, (L_MLTSTART | addflg));
	}
	//multiple point(end)
	else if (cnt == -6) {
		TriCnt = SetInC(Tri, TriCnt, S, E, C1, (L_MLTEND | addflg));
	}
	return TriCnt;
}
//make Control points
BOOL IcasSetCntl(int Type, ICASPOINT Seg[], int *pSegCnt, ICASTRI Tri[], int *pTriCnt)
{
	int siz = *pSegCnt, cnt, tsiz = 0, i;
	ICASPOINT C1, C2, M;
	ICASWORK  SS, S, E, EE, E2, E3;
	//
	if (siz < 3) return FALSE;
	//
	*pTriCnt = 0;
	//
	ChkParallel(Seg, siz);
	//
	//Start part
	Point2Work(&SS, &(Seg[siz-1]));//Point2Work(ICASPOINT *pP, ICASWORK *pW)
	Point2Work(&S,  &(Seg[0]));
	Point2Work(&E,  &(Seg[1]));
	Point2Work(&EE, &(Seg[2]));

	if (Type) {//Close
		cnt = GetCntl( &SS, &S, &E, &EE, &C1, &C2, &M );
	}
	else {//Open
		if ( 0 == GetCntlTerm( (ICASPOINT*)&S, (ICASPOINT*)&E, (ICASPOINT*)&EE, &C1 ) ) cnt = 1;
		else return FALSE;
	}
	tsiz = SetTriangle(Tri, tsiz, cnt, L_START, (ICASPOINT*)&S, (ICASPOINT*)&E, &C1, &C2, &M );
	//Main part
	for(i=3; i<siz; i++) {
		SS = S;
		S  = E;
		E  = EE;
		Point2Work(&EE, &(Seg[i]));
		//
		cnt  = GetCntl( &SS, &S, &E, &EE, &C1, &C2, &M );
		tsiz = SetTriangle(Tri, tsiz, cnt, L_MAIN, (ICASPOINT*)&S, (ICASPOINT*)&E, &C1, &C2, &M );
	}
	if (Type) {//Close
		//End part
		Point2Work(&E2, &(Seg[0]));
		cnt  = GetCntl( &S, &E, &EE, &E2, &C1, &C2, &M );
		tsiz = SetTriangle(Tri, tsiz, cnt, L_END, (ICASPOINT*)&E, (ICASPOINT*)&EE, &C1, &C2, &M );
		//Bridge part
		Point2Work(&E3, &(Seg[1]));
		cnt  = GetCntl( &E, &EE, &E2, &E3, &C1, &C2, &M );
		tsiz = SetTriangle(Tri, tsiz, cnt, L_BRIDGE, (ICASPOINT*)&EE, (ICASPOINT*)&E2, &C1, &C2, &M );
	}
	else {//Open
		//End part
		if ( 0 == GetCntlTerm( (ICASPOINT*)&EE, (ICASPOINT*)&E, (ICASPOINT*)&S, &C1 ) ) cnt = 1;
		else return FALSE;
		tsiz = SetTriangle(Tri, tsiz, cnt, L_END, (ICASPOINT*)&E, (ICASPOINT*)&EE, &C1, &C2, &M );
	}
	//
	*pTriCnt = tsiz;
	return TRUE;
}
BOOL IcasOpen(ICASPOINT Seg[], int *pSegCnt, ICASTRI Tri[], int *pTriCnt)
{
	int siz = *pSegCnt, tsiz = *pTriCnt, cnt, i;
	ICASPOINT C1;
	ICASWORK  S, E, EE;//,SS;
	//
	if (siz < 3) return FALSE;
	//
	if (!(Tri[tsiz-1].Typ & L_BRIDGE)) return FALSE; //already open, perhaps...
	//
	if (   fabs(Tri[0].P[0].x - Tri[tsiz-1].P[2].x) > EPS
		|| fabs(Tri[0].P[0].y - Tri[tsiz-1].P[2].y) > EPS ) return FALSE;

	//Start part
	Point2Work(&S,  &(Seg[0]));//Point2Work(ICASPOINT *pP, ICASWORK *pW)
	Point2Work(&E,  &(Seg[1]));
	Point2Work(&EE, &(Seg[2]));

	cnt = GetCntlTerm( (ICASPOINT*)&S, (ICASPOINT*)&E, (ICASPOINT*)&EE, &C1 );
	if (0 == cnt) {
		if (Tri[0].Typ & L_AUTO_R) {
			//delete Tri[1]
			for(i=1; i<tsiz; ++i)
			{
				Tri[i] = Tri[i+1];
			}
			tsiz--;
		}
		Tri[0].P[0] = Seg[0];
		Tri[0].P[1] = C1;
		Tri[0].P[2] = Seg[1];
		Tri[0].Typ  = L_NORMAL;
	}
	else return FALSE;

	//End part
	Point2Work(&S,  &(Seg[siz-3]));
	Point2Work(&E,  &(Seg[siz-2]));
	Point2Work(&EE, &(Seg[siz-1]));//last

	cnt = GetCntlTerm( (ICASPOINT*)&S, (ICASPOINT*)&E, (ICASPOINT*)&EE, &C1 );
	if (0 == cnt) {
		cnt = tsiz -3;//delete max
		for(i=tsiz-1; i>=cnt; --i) {
			if (Tri[i].Typ & (L_AUTO_R | L_AUTO_L)) {
				tsiz--;
			}
		}
		Tri[tsiz-1].P[0] = Seg[siz-2];
		Tri[tsiz-1].P[1] = C1;
		Tri[tsiz-1].P[2] = Seg[siz-1];
		Tri[tsiz-1].Typ  = L_NORMAL;
	}
	else return FALSE;

	*pTriCnt = tsiz;
	return TRUE;
}
BOOL IcasClose(ICASPOINT Seg[], int *pSegCnt, ICASTRI Tri[], int *pTriCnt)
{
	int siz = *pSegCnt, tsiz = *pTriCnt, cnt, i;
	ICASPOINT C1, C2, M;
	ICASWORK  SS, S, E, EE, E2, E3;
	//
	if (siz < 3) return FALSE;
	//
	if (Tri[tsiz-1].Typ & L_BRIDGE) return FALSE; //already close, perhaps...
	//
	if (   fabs(Tri[0].P[0].x - Tri[tsiz-1].P[2].x) < EPS
		&& fabs(Tri[0].P[0].y - Tri[tsiz-1].P[2].y) < EPS ) return FALSE;

	//Start part
	Point2Work(&SS, &(Seg[siz-1]));//Point2Work(ICASPOINT *pP, ICASWORK *pW)
	Point2Work(&S,  &(Seg[0]));
	Point2Work(&E,  &(Seg[1]));
	Point2Work(&EE, &(Seg[2]));

	cnt = GetCntl( &SS, &S, &E, &EE, &C1, &C2, &M );
	if (1 == cnt) {
		Tri[0].P[1] = C1;
		Tri[0].Typ  = L_NORMAL;
	}
	else if (2 == cnt) {
		//Insert
		//make space //to be m_Tri[i] == m_Tri[i+1]
		for(i=tsiz; i>0; --i)
		{
			Tri[i] = Tri[i-1];
		}
		tsiz++;
		//
		Tri[0].P[0] = Seg[0];
		Tri[0].P[1] = C1;
		Tri[0].P[2] = M;
		Tri[0].Typ  = L_NORMAL | L_AUTO_R;
		Tri[1].P[0] = M;
		Tri[1].P[1] = C2;
		Tri[1].P[2] = Seg[1];
		Tri[1].Typ  = L_NORMAL | L_AUTO_L;
		//tsiz = SetTriangle(Tri, tsiz, cnt, L_START, (ICASPOINT*)&S, (ICASPOINT*)&E, &C1, &C2, &M );
	}
	else return FALSE;

	//End part
	Point2Work(&S,  &(Seg[siz-3]));
	Point2Work(&E,  &(Seg[siz-2]));
	Point2Work(&EE, &(Seg[siz-1]));//last
	Point2Work(&E2, &(Seg[0]));

	cnt  = GetCntl( &S, &E, &EE, &E2, &C1, &C2, &M );
	if (1 == cnt) {
		Tri[tsiz-1].P[0] = Seg[siz-2];
		Tri[tsiz-1].P[1] = C1;
		Tri[tsiz-1].P[2] = Seg[siz-1];//last
		Tri[tsiz-1].Typ  = L_NORMAL;
	}
	else if (2 == cnt) {
		tsiz = SetTriangle(Tri, tsiz, cnt, L_END, (ICASPOINT*)&E, (ICASPOINT*)&EE, &C1, &C2, &M );
	}
	else return FALSE;

	//Bridge part
	Point2Work(&E3, &(Seg[1]));
	cnt  = GetCntl( &E, &EE, &E2, &E3, &C1, &C2, &M );
	tsiz = SetTriangle(Tri, tsiz, cnt, L_BRIDGE, (ICASPOINT*)&EE, (ICASPOINT*)&E2, &C1, &C2, &M );

	*pTriCnt = tsiz;
	return TRUE;
}
BOOL MoveRecord(ICASTRI Tri[], int TriCnt,
				ICASPOINT Seg[], int SegCnt, ICASPOINT *D)
{
	int i;
	for(i=0; i<SegCnt; ++i) {
		Seg[i].x = Seg[i].x + D->x;
		Seg[i].y = Seg[i].y + D->y;
	}
	for(i=0; i<TriCnt; ++i) {
		Tri[i].P[0].x = Tri[i].P[0].x + D->x;
		Tri[i].P[0].y = Tri[i].P[0].y + D->y;
		Tri[i].P[1].x = Tri[i].P[1].x + D->x;
		Tri[i].P[1].y = Tri[i].P[1].y + D->y;
		Tri[i].P[2].x = Tri[i].P[2].x + D->x;
		Tri[i].P[2].y = Tri[i].P[2].y + D->y;
	}
	return TRUE;
}

BOOL RotateRecord(ICASTRI Tri[], int TriCnt,
				  ICASPOINT Seg[], int SegCnt, int iD, 
				  ICASPOINT *D, int init)
{
	static ICASPOINT Mat;
	ICASPOINT u, w;
	ICASPOINT o = Seg[iD];//rotate origin 
	ICASPOINT v = *D;
	double d;
	int i,j;
	v.x -= o.x;//
	v.y -= o.y;
	d = sqrt(v.x*v.x + v.y*v.y);
	if (fabs(d) < EPS ) return FALSE;
	v.x = v.x / d;//unit vector
	v.y = v.y / d;
	if (init) {
		Mat.x = v.x;
		Mat.y = v.y;
		//TRACE("0:[%d]%lf,%lf",m_iD,v.x,v.y);
		return TRUE;
	}
	else {
		u.x = Mat.x * v.x + Mat.y * v.y;
		u.y = Mat.x * v.y - Mat.y * v.x;
		//TRACE("[%d]%lf,%lf\n",m_iD,u.x,u.y);
	}
	for(i=0; i<SegCnt; ++i) {
		if (i == iD) continue;
		w.x = Seg[i].x - o.x;
		w.y = Seg[i].y - o.y;
		Seg[i].x = u.x * w.x - u.y * w.y + o.x;
		Seg[i].y = u.y * w.x + u.x * w.y + o.y;
	}
	for(i=0; i<TriCnt; ++i) {
		for(j=0; j<3; ++j) {
			w.x = Tri[i].P[j].x - o.x;
			w.y = Tri[i].P[j].y - o.y;
			if(fabs(w.x) < EPS && fabs(w.y) < EPS ) continue;
			Tri[i].P[j].x = u.x * w.x - u.y * w.y + o.x;
			Tri[i].P[j].y = u.y * w.x + u.x * w.y + o.y;
		}
	}
	return TRUE;
}

BOOL MoveVirtualPoint(ICASTRI Tri[], int TriCnt,
					  int iT, int iTR,
					  ICASPOINT *C, ICASPOINT *P, ICASPOINT *Q,
					  ICASPOINT *D, int init)
{
	static ICASPOINT A,B,AB,CD,AC,BD;
	ICASPOINT M;
	double d,s,t;

	if (init) {
		A = Tri[iT ].P[0];
		B = Tri[iTR].P[2];
		AB.x = B.x - A.x;
		AB.y = B.y - A.y;
		d = sqrt(AB.x * AB.x + AB.y * AB.y);
		if (fabs(d) > EPS) {
			AB.x = AB.x / d;//unit vector
			AB.y = AB.y / d;
		}
		else return FALSE;

		CD.x = Tri[iTR].P[1].x - Tri[iT].P[1].x;
		CD.y = Tri[iTR].P[1].y - Tri[iT].P[1].y;
		AC.x = Tri[iT ].P[1].x - A.x;
		AC.y = Tri[iT ].P[1].y - A.y;
		BD.x = Tri[iTR].P[1].x - B.x;
		BD.y = Tri[iTR].P[1].y - B.y;
		return TRUE;
	}
	M = A;
	d = AB.x * (D->x - M.x) + AB.y * (D->y - M.y);//Inner product
	//New Auto = M
	M.x += d * AB.x;//scale unit vector
	M.y += d * AB.y;
	M.typ = L_NORMAL;
	*C = M;

	d = AC.y * CD.x - AC.x * CD.y;
	if (fabs(d) < EPS) return FALSE;
	s = ((M.y - A.y)*CD.x - (M.x - A.x)*CD.y)/d;

	d = BD.y * CD.x - BD.x * CD.y;
	if (fabs(d) < EPS) return FALSE;
	t = ((M.y - B.y)*CD.x - (M.x - B.x)*CD.y)/d;

	P->x = A.x + s * AC.x;
	P->y = A.y + s * AC.y;
	Q->x = B.x + t * BD.x;
	Q->y = B.y + t * BD.y;
	P->typ = L_NORMAL;
	Q->typ = L_NORMAL;
	return TRUE;
}

int ChkLinear(ICASPOINT Seg[], int iDL, int iDR, ICASPOINT *D)
{
	int sts = L_NORMAL;
	double Ux,Uy,d1,Vx,Vy,d2;
	Ux = Seg[iDL].x - D->x;
	Uy = Seg[iDL].y - D->y;
	d1 = Ux*Ux + Uy*Uy;
	if (fabs(d1) < EPS) {
		//TRACE("ChkLiner1:L_MULTIPLE\n");
		sts = L_MULTIPLE;
	}
	else {
		Vx = Seg[iDR].x - D->x;
		Vy = Seg[iDR].y - D->y;
		d2 = Vx*Vx + Vy*Vy;
		if (fabs(d2) < EPS) {
			//TRACE("ChkLiner2:L_MULTIPLE\n");
			sts = L_MULTIPLE;
		}
		else {
			d1 = Ux*Vy - Vx*Uy;
			if (fabs(d1) < EPS ) {
				//TRACE("ChkLiner3:L_PARALLEL\n");
				sts = L_PARALLEL;
			}
		}
	}
	return sts;
}

void CalcInsertCntl(ICASTRI Tri[], int iT,
					ICASPOINT *M, ICASPOINT *C)
{
	M->x = (Tri[iT].P[0].x + Tri[iT].P[2].x)*0.5;
	M->y = (Tri[iT].P[0].y + Tri[iT].P[2].y)*0.5;
	C->x = (M->x - Tri[iT].P[1].x)*0.5;
	C->y = (M->y - Tri[iT].P[1].y)*0.5;
	C->x += M->x;
	C->y += M->y;
	M->typ = L_NORMAL;
	C->typ = L_NORMAL;
}

void MidPoint(ICASPOINT Seg[], int iDS, int iDE,
			  ICASPOINT *M )
{
	M->x = (Seg[iDS].x + Seg[iDE].x)*0.5;
	M->y = (Seg[iDS].y + Seg[iDE].y)*0.5;
	M->typ = L_NORMAL;
}

/////////////////////////////////////////////////////////////////////////////