/* MPMorph - Amiga Morphing program */
/* Copyright (C) © 1993-97 Mark John Paddock */

/* This program is free software; you can redistribute it and/or modify */
/* it under the terms of the GNU General Public License as published by */
/* the Free Software Foundation; either version 2 of the License, or */
/* any later version. */

/* This program is distributed in the hope that it will be useful, */
/* but WITHOUT ANY WARRANTY; without even the implied warranty of */
/* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the */
/* GNU General Public License for more details. */

/* You should have received a copy of the GNU General Public License */
/* along with this program; if not, write to the Free Software */
/* Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA. */

/* mpaddock@cix.compulink.co.uk */

#include "Render.h"

#define CalcXYdiff(MyPoint) {\
	if (x < MyPoint->Cx) {\
  		MyPoint->xdiff = MyPoint->Cx - x;\
   }\
   else {\
   	MyPoint->xdiff = x - MyPoint->Cx;\
   }\
	if (y < MyPoint->Cy) {\
		MyPoint->ydiff = MyPoint->Cy - y;\
	}\
	else {\
		MyPoint->ydiff = y - MyPoint->Cy;\
  	}\
}

#define CalcCdiff(MyPoint) {\
	if (x < MyPoint->Cx) {\
  		MyPoint->xdiff = MyPoint->Cx - x;\
   }\
   else {\
   	MyPoint->xdiff = x - MyPoint->Cx;\
   }\
	if (y < MyPoint->Cy) {\
		MyPoint->ydiff = MyPoint->Cy - y;\
	}\
	else {\
		MyPoint->ydiff = y - MyPoint->Cy;\
  	}\
  	MyPoint->Cdiff = MyPoint->xdiff * MyPoint->xdiff +\
						  MyPoint->ydiff * MyPoint->ydiff;\
}

#define CalcCdiffQ(MyPoint) {\
  	MyPoint->Cdiff = MyPoint->xdiff * MyPoint->xdiff +\
						  MyPoint->ydiff * MyPoint->ydiff;\
}

#define CheckPoint(RetPoint) {\
  	if (AntiAlias) {\
	   if (RetPoint->xd < 0) {\
	   	RetPoint->xd = 0;\
	   }\
	   if (RetPoint->yd < 0) {\
	   	RetPoint->yd = 0;\
	   }\
	   if (RetPoint->x1d < 0) {\
	   	RetPoint->x1d = 0;\
	   }\
	   if (RetPoint->y1d < 0) {\
	   	RetPoint->y1d = 0;\
	   }\
	   if (RetPoint->xd > (width-1)) {\
	   	RetPoint->xd = width-1;\
	   }\
	   if (RetPoint->yd > (height-1)) {\
   		RetPoint->yd = height-1;\
	   }\
	   if (RetPoint->x1d > (width-1)) {\
	   	RetPoint->x1d = width-1;\
	   }\
	   if (RetPoint->y1d > (height-1)) {\
	   	RetPoint->y1d = height-1;\
	   }\
  	}\
  	else {\
	   if (RetPoint->x < 0) {\
	   	RetPoint->x = 0;\
	   }\
	   if (RetPoint->y < 0) {\
	   	RetPoint->y = 0;\
	   }\
	   if (RetPoint->x1 < 0) {\
	   	RetPoint->x1 = 0;\
	   }\
	   if (RetPoint->y1 < 0) {\
	   	RetPoint->y1 = 0;\
	   }\
	   if (RetPoint->x > (width-1)) {\
	   	RetPoint->x = width-1;\
	   }\
	   if (RetPoint->y > (height-1)) {\
   		RetPoint->y = height-1;\
	   }\
	   if (RetPoint->x1 > (width-1)) {\
	   	RetPoint->x1 = width-1;\
	   }\
	   if (RetPoint->y1 > (height-1)) {\
	   	RetPoint->y1 = height-1;\
	   }\
	}\
}

/* Point is on a point... */
#define ReturnAPoint(Pointa,Pointx) {\
  	if (AntiAlias) {\
  		Pointx->xd = Pointa->x;\
  		Pointx->yd = Pointa->y;\
  		Pointx->x1d = Pointa->x1;\
  		Pointx->y1d = Pointa->y1;\
  	}\
  	else {\
  		Pointx->x = Pointa->x;\
  		Pointx->y = Pointa->y;\
  		Pointx->x1 = Pointa->x1;\
  		Pointx->y1 = Pointa->y1;\
  	}\
}

/* Point is on a line... */
#define Return2Points(a,c,p,x,y) {\
	if (AntiAlias) {\
		double AB,cs,as;\
		AB = (as = sqrt((double)a->Cdiff)) + (cs = sqrt((double)c->Cdiff));\
		p->xd = x +\
				((a->x - a->Cx) * cs +\
				 (c->x - c->Cx) * as) /\
				 AB;\
		p->x1d = x +\
				((a->x1 - a->Cx) * cs +\
				 (c->x1 - c->Cx) * as) /\
				AB;\
		p->yd = y +\
				((a->y - a->Cy) * cs +\
				 (c->y - c->Cy) * as) /\
				AB;\
		p->y1d = y +\
				((a->y1 - a->Cy) * cs +\
				 (c->y1 - c->Cy) * as) /\
				AB;\
	}\
	else {\
		LONG AB;\
		AB = a->Cdiff + c->Cdiff;\
		p->x = x +\
				((a->x - a->Cx) * c->Cdiff +\
				 (c->x - c->Cx) * a->Cdiff) /\
				 AB;\
		p->x1 = x +\
				((a->x1 - a->Cx) * c->Cdiff +\
				 (c->x1 - c->Cx) * a->Cdiff) /\
				AB;\
		p->y = y +\
				((a->y - a->Cy) * c->Cdiff +\
				 (c->y - c->Cy) * a->Cdiff) /\
				AB;\
		p->y1 = y +\
				((a->y1 - a->Cy) * c->Cdiff +\
				 (c->y1 - c->Cy) * a->Cdiff) /\
				AB;\
	}\
	CheckPoint(p);\
}

/* Caclulate the distances between various points		*/
	/* ab is length x,y to line ab */
	/* ab1 is length x,y to intersect of ab */
	/* ab2 is length intersect of ab to bb */
#define CalcDiffs(a,bb,ab1,ab2,ab,x,y) {\
	if (a->Cx == bb->Cx) {\
		ab  = a->xdiff;\
		ab1 = a->ydiff;\
		ab2 = bb->ydiff;\
	}\
	else {\
		if (a->Cy == bb->Cy) {\
			ab1 = a->xdiff;\
			ab2 = bb->xdiff;\
			ab  = a->ydiff;\
		}\
		else {\
			double m;\
			double c;		/* y = mx+c line a to bb */\
			double n;\
			double d;		/* y = nx+d line perp to y=mx+c thru x,y */\
			double xsolve,ysolve;\
			m = (double)(bb->Cy - a->Cy) / (double)(bb->Cx - a->Cx);\
			c = bb->Cy - (bb->Cx * m);\
			n = (-1) / m;\
			d = y - n * x;\
			xsolve = (d - c)/(m - n);\
			ysolve = m * xsolve + c;\
			ab1 = sqrt(((a->Cy - ysolve) * (a->Cy - ysolve) + (a->Cx - xsolve) * (a->Cx - xsolve)));\
			ab2 = sqrt(((bb->Cy - ysolve) * (bb->Cy - ysolve) + (bb->Cx - xsolve) * (bb->Cx - xsolve)));\
			ab  = sqrt(((y - ysolve) * (y - ysolve) + (x - xsolve) * (x - xsolve)));\
		}\
	}\
}

/* Used to determine if a point lies on a line,
 * or in triangle of 3 points
 */
/* Draws line between a and b and sees if line from x,0 to x,y intersects */
/* returns 1 if intersects, 0 if not intersect */
/* Version 3.3 returns 4 if on a line */
/* This used to be a function - now very complex define! */
#define intersect(Pointa,Pointb,x,y) (\
	(Pointb->Cx == Pointa->Cx) ?								/* Vertical a to b */\
		((Pointb->Cx == x) &&									/* x,y on a to b	3.3 fix!!!! */\
		 (((Pointb->Cy > y) && (Pointa->Cy < y)) ||\
		  ((Pointb->Cy < y) && (Pointa->Cy > y)))) ?\
			4															/* Between a and b */\
		:\
			0\
	:\
		(((Pointb->Cy > y) && (Pointa->Cy > y)) ||\
		 ((Pointb->Cx > x) && (Pointa->Cx > x)) ||\
		 ((Pointb->Cx < x) && (Pointa->Cx < x))) ?		/* Can not possibly cross */\
			0\
		:\
			((Pointb->Cy < y) && (Pointa->Cy < y)) ?\
				1\
			:\
				(Pointa->Cx == x) ?\
					(Pointa->Cy > y) ?\
						0\
					:\
						1\
				:\
					(Pointb->Cx == x) ?\
						(Pointb->Cy > y) ?\
							0\
						:\
							1\
					:\
						(Pointb->Cy == Pointa->Cy) ?\
							(Pointa->Cy > y) ?\
								0\
							:\
								1\
						:\
							((temp = (((Pointa->Cx < x)?(x - Pointa->Cx):(Pointa->Cx - x)) * Pointb->Cy +\
								  ((Pointb->Cx < x)?(x - Pointb->Cx):(Pointb->Cx - x)) * Pointa->Cy)) >\
								  (temp1 = ((Pointb->Cx < Pointa->Cx)?(Pointa->Cx - Pointb->Cx):(Pointb->Cx - Pointa->Cx)) * y)) ?\
								0\
							:\
								(temp == temp1) ?\
									4\
								:\
									1\
)

/* Given 3 points (a,b,c) with a point (p) in the triangle
 * determine the coordinates of point (p) on the first and last frame
 */
#define Triangle(p,a,bb,cc,x,y) {\
	double AB1,AB2,BC1,BC2,CA1,CA2;\
	double AB,BC,CA;\
	double ABBCCA2;\
	double ABCA,BCCA,ABBC;\
	double BC1BC2,AB1AB2,CA1CA2;\
	CalcDiffs(a,bb,AB1,AB2,AB,x,y);\
	CalcDiffs(bb,cc,BC1,BC2,BC,x,y);\
	CalcDiffs(cc,a,CA1,CA2,CA,x,y);\
	ABBCCA2 = 2 * (AB + BC + CA);\
	ABCA = AB+CA;\
	BCCA = BC+CA;\
	ABBC = AB+BC;\
	BC1BC2 = BC1+BC2;\
	AB1AB2 = AB1+AB2;\
	CA1CA2 = CA1+CA2;\
	if (AntiAlias) {\
		p->xd = x +\
					(ABCA*(BC1*(cc->x - cc->Cx) + BC2*(bb->x - bb->Cx)) / BC1BC2 +\
					 BCCA*(AB1*(bb->x - bb->Cx) + AB2*(a->x - a->Cx)) / AB1AB2 +\
					 ABBC*(CA1*(a->x - a->Cx) + CA2*(cc->x - cc->Cx)) / CA1CA2) /\
					ABBCCA2;\
		p->x1d = x +\
					(ABCA*(BC1*(cc->x1 - cc->Cx) + BC2*(bb->x1 - bb->Cx)) / BC1BC2 +\
					 BCCA*(AB1*(bb->x1 - bb->Cx) + AB2*(a->x1 - a->Cx)) / AB1AB2 +\
					 ABBC*(CA1*(a->x1 - a->Cx) + CA2*(cc->x1 - cc->Cx)) / CA1CA2) /\
					ABBCCA2;\
		p->y1d = y +\
					(ABCA*(BC1*(cc->y1 - cc->Cy) + BC2*(bb->y1 - bb->Cy)) / BC1BC2 +\
					 BCCA*(AB1*(bb->y1 - bb->Cy) + AB2*(a->y1 - a->Cy)) / AB1AB2 +\
					 ABBC*(CA1*(a->y1 - a->Cy) + CA2*(cc->y1 - cc->Cy)) / CA1CA2) /\
					ABBCCA2;\
		p->yd = y +\
					(ABCA*(BC1*(cc->y - cc->Cy) + BC2*(bb->y - bb->Cy)) / BC1BC2 +\
					 BCCA*(AB1*(bb->y - bb->Cy) + AB2*(a->y - a->Cy)) / AB1AB2 +\
					 ABBC*(CA1*(a->y - a->Cy) + CA2*(cc->y - cc->Cy)) / CA1CA2) /\
					ABBCCA2;\
	}\
	else {\
		p->x = x +\
					(ABCA*(BC1*(cc->x - cc->Cx) + BC2*(bb->x - bb->Cx)) / BC1BC2 +\
					 BCCA*(AB1*(bb->x - bb->Cx) + AB2*(a->x - a->Cx)) / AB1AB2 +\
					 ABBC*(CA1*(a->x - a->Cx) + CA2*(cc->x - cc->Cx)) / CA1CA2) /\
					ABBCCA2;\
		p->x1 = x +\
					(ABCA*(BC1*(cc->x1 - cc->Cx) + BC2*(bb->x1 - bb->Cx)) / BC1BC2 +\
					 BCCA*(AB1*(bb->x1 - bb->Cx) + AB2*(a->x1 - a->Cx)) / AB1AB2 +\
					 ABBC*(CA1*(a->x1 - a->Cx) + CA2*(cc->x1 - cc->Cx)) / CA1CA2) /\
					ABBCCA2;\
		p->y1 = y +\
					(ABCA*(BC1*(cc->y1 - cc->Cy) + BC2*(bb->y1 - bb->Cy)) / BC1BC2 +\
					 BCCA*(AB1*(bb->y1 - bb->Cy) + AB2*(a->y1 - a->Cy)) / AB1AB2 +\
					 ABBC*(CA1*(a->y1 - a->Cy) + CA2*(cc->y1 - cc->Cy)) / CA1CA2) /\
					ABBCCA2;\
		p->y = y +\
					(ABCA*(BC1*(cc->y - cc->Cy) + BC2*(bb->y - bb->Cy)) / BC1BC2 +\
					 BCCA*(AB1*(bb->y - bb->Cy) + AB2*(a->y - a->Cy)) / AB1AB2 +\
					 ABBC*(CA1*(a->y - a->Cy) + CA2*(cc->y - cc->Cy)) / CA1CA2) /\
					ABBCCA2;\
	}\
}

/* Given 3 points (a,b,c) with a point (p) in the triangle
 * determine the coordinates of point (p) on the first and last frame
 * Integer version
 * Completely different to the real version
 */
#define Triangle_I(p,a,bb,cc,x,y) {\
	LONG ABBCCA2;\
	LONG bc,ac,ab;\
	ABBCCA2 = 2 * ((ab = a->Cdiff + bb->Cdiff) + cc->Cdiff);\
	p->x = x +\
				((a->x - a->Cx) * (bc = (bb->Cdiff + cc->Cdiff)) +\
				 (bb->x - bb->Cx) * (ac = (a->Cdiff + cc->Cdiff)) +\
				 (cc->x - cc->Cx) * (ab)) /\
				 ABBCCA2;\
	p->x1 = x +\
				((a->x1 - a->Cx) * (bc) +\
				 (bb->x1 - bb->Cx) * (ac) +\
				 (cc->x1 - cc->Cx) * (ab)) /\
				ABBCCA2;\
	p->y = y +\
				((a->y - a->Cy) * (bc) +\
				 (bb->y - bb->Cy) * (ac) +\
				 (cc->y - cc->Cy) * (ab)) /\
				ABBCCA2;\
	p->y1 = y +\
				((a->y1 - a->Cy) * (bc) +\
				 (bb->y1 - bb->Cy) * (ac) +\
				 (cc->y1 - cc->Cy) * (ab)) /\
				ABBCCA2;\
}

/* The following is the alogrithm stuff!!!!! */
/* and is not particularly well commented! */

/* version 3.3 works differntly with points on lines/points */

/* Given a point, returns the coordinates in the first and last image	*/
void __regargs
FindPoint(struct MyPoint *RetPoint,BOOL Everything,WORD x,WORD y) {
	struct MyPoint *MyPoint,*MyPoint1,*MyPoint2,*MyPoint3;
	struct MyPoint **p,**p1,**p2;
	UWORD NumPoints;
	ULONG temp,temp1;	/* Used by intersect macro */

	if (!(Mode & MODE_ONCE) || !Everything) {
		/* Normal calculation - or pre-calculation for pre-calculate */
		{
			/* Not search mode */
			if ((Mode & MODE_DELAU) && (Mode & MODE_STAT)) {
				p = Points;
				*p = &BigPoint;		// 4.7
				for (MyPoint = (struct MyPoint *)PointList.lh_Head;
						MyPoint->MyNode.mln_Succ;
		   			MyPoint = (struct MyPoint *)MyPoint->MyNode.mln_Succ) {
		   		CalcXYdiff(MyPoint);
					if ((MyPoint->xdiff + MyPoint->ydiff) < (*p)->Cdiff) {
						CalcCdiffQ(MyPoint);
						if (Everything && !MyPoint->Cdiff) {
							/* On a point and not pre-calculating */
							ReturnAPoint(MyPoint,RetPoint);
							return;
						}
						if (MyPoint->Cdiff < (*p)->Cdiff) {
							*p = MyPoint;
						}
					}
				}
				
			}
			else {
				if (!Depth) {
					Points[0] = &BigPoint;
					Points[1] = &BigPoint;
					Points[2] = &BigPoint;
					MyPoint1 = Points[0];
					MyPoint2 = Points[1];
					MyPoint3 = Points[2];
					for (MyPoint = (struct MyPoint *)PointList.lh_Head;
							MyPoint->MyNode.mln_Succ;
			   			MyPoint = (struct MyPoint *)MyPoint->MyNode.mln_Succ) {
		   			CalcXYdiff(MyPoint);
						if ((MyPoint->xdiff + MyPoint->ydiff) < MyPoint3->Cdiff) {
							CalcCdiffQ(MyPoint);
							if (Everything && !MyPoint->Cdiff) {
								/* On a point and not pre-calculating */
								ReturnAPoint(MyPoint,RetPoint);
								return;
							}
							if (MyPoint->Cdiff < MyPoint3->Cdiff) {
								if (MyPoint->Cdiff < MyPoint1->Cdiff) {
									MyPoint3 = Points[2] = Points[1];
									MyPoint2 = Points[1] = Points[0];
									MyPoint1 = Points[0] = MyPoint;
								}
								else {
									if (MyPoint->Cdiff < MyPoint2->Cdiff) {
										MyPoint3 = Points[2] = Points[1];
										MyPoint2 = Points[1] = MyPoint;
									}
									else {
										MyPoint3 = Points[2] = MyPoint;
									}
								}
							}
						}
					}
				}
				else {
					UWORD i;
					NumPoints = Depth+3;
					p = Points;
					for (i = 0;
						  i < NumPoints;
						  ++i) {
						*(p++) = &BigPoint;
					}
					/* Set up differences */
					p2 = &(Points[Depth+2]);
					for (MyPoint = (struct MyPoint *)PointList.lh_Head;
							MyPoint->MyNode.mln_Succ;
			   			MyPoint = (struct MyPoint *)MyPoint->MyNode.mln_Succ) {
			   		CalcXYdiff(MyPoint);
						if ((MyPoint->xdiff + MyPoint->ydiff) < (*p2)->Cdiff) {
							CalcCdiffQ(MyPoint);
							if (Everything && !MyPoint->Cdiff) {
								/* On a point and not pre-calculating */
								ReturnAPoint(MyPoint,RetPoint);
								return;
							}
							if (MyPoint->Cdiff < (*p2)->Cdiff) {
								p = Points;
								for (i = 0;
									  (i < NumPoints);
									  ++i,++p) {
									if ((*p)->Cdiff > MyPoint->Cdiff) {
										if (Depth+2-i) {
											memmove(p+1,p,(Depth+2-i)*sizeof(struct MyPoint *));
										}
										*p = MyPoint;
										break;
									}
								}
							}
						}
					}
				}
			}
		}
		if (!Everything) {
			if (!(Mode & MODE_DELAU)) {
				Points[Depth+3] = (struct MyPoint *)NumPoints;
			}
			return;
		}
	}
	else {
		/* things have been pre-calculated */
		if (Mode & MODE_DELAU) {
			if (((*Points)->Cx == x) &&
				 ((*Points)->Cy == y)) {
				ReturnAPoint((*Points),RetPoint);
				return;
			}
		}
		else {
			UWORD i;
			p = Points;
			NumPoints = (ULONG)Points[Depth+3];
			for (i = 0;
					((i < Depth+3) && *p);
					++i,++p) {
				MyPoint = *p;
				CalcCdiff(MyPoint);
				if (!MyPoint->Cdiff) {
					ReturnAPoint(MyPoint,RetPoint);
					return;
				}
			}
		}
	}
   /* otherwise try and find a triangle with the other points */
	if (!Depth && !(Mode & (MODE_STAT | MODE_DELAU))) {
		/* Special case, only have 3 points */
		MyPoint = Points[0];
		MyPoint1 = Points[1];
		MyPoint2 = Points[2];
	}
	else {
	   if (Mode & MODE_DELAU) {
			struct Triangle *T;
	   	MyPoint = Points[0];
			for (T = (struct Triangle *)MyPoint->TList.lh_Head;
					T->TNode.mln_Succ;
	   			T = (struct Triangle *)T->TNode.mln_Succ) {
	   		if (!((x < T->MinX) || (x > T->MaxX) || (y < T->MinY) || (y > T->MaxY))) {
					UWORD i1,i2,i3;
					if ((i1 = intersect(MyPoint,T->Point1,x,y)) == 4) {
						CalcCdiff(MyPoint);
						CalcCdiff((T->Point1));
						Return2Points(MyPoint,(T->Point1),(RetPoint),x,y);
						return;
					}
			   	if (i1 + (i2 = intersect(T->Point1,T->Point2,x,y)) < 2) {
						if ((i1 + i2 + (i3 = intersect(T->Point2,MyPoint,x,y))) == 1) {
							MyPoint1 = T->Point1;
							MyPoint2 = T->Point2;
							if (Integer) {
								CalcCdiff(MyPoint);
								CalcCdiff(MyPoint1);
								CalcCdiff(MyPoint2);
							}
							else {
								CalcXYdiff(MyPoint);
								CalcXYdiff(MyPoint1);
								CalcXYdiff(MyPoint2);
							}
							goto FoundPoints;
						}
						else {
							if (i3 == 4) {
								CalcCdiff(MyPoint);
								CalcCdiff((T->Point2));
								Return2Points((T->Point2),MyPoint,RetPoint,x,y);
								return;
							}
						}
					}
					else {
						if (i2 == 4) {
							CalcCdiff((T->Point1));
							CalcCdiff((T->Point2));
							Return2Points((T->Point1),(T->Point2),RetPoint,x,y);
							return;
						}
					}
				}
	   	}
	   	if (!(Mode & MODE_STAT)) {
	   		if (Integer) {
		   		CalcCdiff(MyPoint);
		   		CalcCdiff((Points[1]));
	   			CalcCdiff((Points[2]));
	   		}
	   		else {
		   		CalcXYdiff(MyPoint);
		   		CalcXYdiff((Points[1]));
	   			CalcXYdiff((Points[2]));
	   		}
	   	}
	   }
	   else {
		  	if (Mode & MODE_CLOSEST) {
				UWORD KeepCount = 0;
		  		UWORD i;
				LONG checkdiff;
				checkdiff = 0x7FFFFFFF;
				for (p = Points, i=0;
					  i<(NumPoints-2);
					  ++i,++p) {
					if ((*p)->Cdiff < checkdiff) {
						UWORD j;
						for (p1 = p+1, j=i;
							  j<(NumPoints-2);
							  ++j,++p1) {
							if (((*p)->Cdiff + (*p1)->Cdiff) < checkdiff) {
								UWORD i1,k;
								if ((i1 = intersect((*p),(*p1),x,y)) == 4) {
									Return2Points((*p),(*p1),RetPoint,x,y);
									return;
								}
								for (p2 = p1+1, k=j;
									  k<(NumPoints-2);
									  ++k,++p2) {
									/* loop thru 3 point combinations finding smallest triangle */
									if (((*p)->Cdiff + (*p1)->Cdiff + (*p2)->Cdiff) < checkdiff) {
										UWORD i2;
								   	if (i1 + (i2 = intersect((*p1),(*p2),x,y)) < 2) {
								   		UWORD i3;
								   		if ((i1 + i2 + (i3 = intersect((*p2),(*p),x,y))) == 1) {
									   		KeepCount = 1;
									   		MyPoint = *p;
												MyPoint1 = *p1;
												MyPoint2 = *p2;
												checkdiff = (*p)->Cdiff + (*p1)->Cdiff + (*p2)->Cdiff;
											}
											else {
												if (i3 == 4) {
													Return2Points((*p2),(*p),RetPoint,x,y);
													return;
												}
											}
										}
										else {
											if (i2 == 4) {
												Return2Points((*p1),(*p2),RetPoint,x,y);
												return;
											}
										}
									}
						   	}
					   	}
				   	}
			   	}
		  		}
		  		if (KeepCount) {
		  			goto FoundPoints;
		  		}
		  	}
		  	else {
		  		UWORD i;
				for (p = Points, i=0;
					  (i<(NumPoints-2));
					  ++i,++p) {
					UWORD j;
					for (p1 = p+1, j=i;
						  (j<(NumPoints-2));
						  ++j,++p1) {
						UWORD i1,k;
						i1 = intersect((*p),(*p1),x,y);
						if (i1 == 4) {
							Return2Points((*p),(*p1),RetPoint,x,y);
							return;
						}
						for (p2 = p1+1, k=j;
							  (k<(NumPoints-2));
							  ++k,++p2) {
							UWORD i2;
							/* loop thru 3 point combinations finding smallest triangle */
					   	if (i1 + (i2 = intersect((*p1),(*p2),x,y)) < 2) {
					   		UWORD i3;
					   		if ((i1 + i2 + (i3 = intersect((*p2),(*p),x,y))) == 1) {
						   		MyPoint = *p;
						   		MyPoint1 = *p1;
						   		MyPoint2 = *p2;
						   		goto FoundPoints;
						   	}
						   	else {
						   		if (i3 == 4) {
						   			Return2Points((*p2),(*p),RetPoint,x,y);
						   			return;
						   		}
						   	}
					   	}
					   	else {
					   		if (i2 == 4) {
					   			Return2Points((*p1),(*p2),RetPoint,x,y);
					   			return;
					   		}
					   	}
				   	}
			   	}
		  		}
			}
		}
	  	/* If still no triangle then assume stationary */
	  	if (Mode & MODE_STAT) {
	   	if (AntiAlias) {
				RetPoint->xd = x;
				RetPoint->x1d = x;
				RetPoint->yd = y;
				RetPoint->y1d = y;
				return;
			}
	  		RetPoint->x = x;
  			RetPoint->y = y;
  			RetPoint->x1 = x;
	  		RetPoint->y1 = y;
	  		return;
	  	}
  		MyPoint = Points[0];
  		MyPoint1 = Points[1];
  		MyPoint2 = Points[2];
	}
FoundPoints:
	/* Otherwise do the triangle calculation */
	if (Integer) {
		Triangle_I(RetPoint,MyPoint,MyPoint1,MyPoint2,x,y);
	}
	else {
		Triangle(RetPoint,MyPoint,MyPoint1,MyPoint2,x,y);
	}
   /* Check the point is actually on the image */
   CheckPoint(RetPoint);
}