#version 3.8;
global_settings {
	assumed_gamma 1.0
}


#include "math.inc"

#declare E = 0.000001;

#declare Camera_Orthographic = true;		// <----################################################################################

#declare Camera_Position = <0, 0, -20>;		// front view
#declare Camera_Look_At  = <0, 0,  0> ;
#declare Fraction = 0.0001;					// functions as a zoom for the orthographic view: 4 zooms in 4x, 8 zooms in 8x, etc.

#declare Aspect = image_width/image_height;
// ###########################################
camera {
	#if (Camera_Orthographic = true)
		orthographic
		right     x*image_width/(Fraction)
		up 		y*image_height/(Fraction)
	#else
		right     x*image_width/image_height
		up y
	#end
	location  Camera_Position
	look_at   Camera_Look_At} 
// ###########################################
	
sky_sphere {pigment {rgb 1}}
light_source {<0, 0, -image_width/Fraction> color rgb 1}







#macro getNormal (p1, p2, p3, Normal)

	#local v1 = <p2[0]-p1[0], p2[1]-p1[1], p2[2]-p1[2]>;
	#local v2 = <p3[0]-p1[0], p3[1]-p1[1], p3[2]-p1[2]>;
	#local NormalX = v1[1]*v2[2] - v1[2]*v2[1];
	#local NormalY = v1[2]*v2[0] - v1[0]*v2[2];
	#local NormalZ = v1[0]*v2[1] - v1[1]*v2[0];
	#local Normal = <NormalX, NormalY, NormalZ>;
	
	Normal
	
#end

/*
void normalizeVector(#local vec[3]){
	#local mag=pow(vec[0]*vec[0]+vec[1]*vec[1]+vec[2]*vec[2],0.5);
	vec[0]=vec[0]/mag;
	vec[1]=vec[1]/mag;
	vec[2]=vec[2]/mag;
	return;
}

#local dot(#local x[3], #local y[3]){
	return x[0]*y[0]+x[1]*y[1]+x[2]*y[2];
}

void cross(#local x[3], #local y[3], #local out[3]){
	out[0]=x[1]*y[2]-x[2]*y[1];
	out[1]=x[2]*y[0]-x[0]*y[2];
	out[2]=x[0]*y[1]-x[1]*y[0];
	return;
}


void vectorBetweenNodes(struct NODE* tail, struct NODE* head, #local out[3]){
	out[0]=head->x-tail->x;
	out[1]=head->y-tail->y;
	out[2]=head->z-tail->z;
	return;
}

#local norm(#local  type, #local x[3]){
	if (type==2)
		return pow(dot(x,x),0.5);
	pr#local f("OOPS! Bad type in norm(...).\n");
	return 0;
}

// angle between nodes in 3d space
#local getNodeAngle(struct NODE* node1, struct NODE* node2, struct NODE* node3, #local normal[3]){
	#local v1[3];
	#local v3[3];
	vectorBetweenNodes(node2,node1,v1);
	vectorBetweenNodes(node2,node3,v3);
	#local cp[3];
	#local tempNormal[3]={normal[0],normal[1],normal[2]};
	cross(v1,v3,cp);
	normalizeVector(tempNormal);
	#local ang=atan2(dot(cp,tempNormal),dot(v1,v3));
	return ang;
}

// angle between points in 3d space
#local getpoint Angle(#local a[3], #local b[3], #local c[3], #local normal[3]){
	#local v1[3]={b[0]-a[0],b[1]-a[1],b[2]-a[2]};
	#local v3[3]={c[0]-a[0],c[1]-a[1],c[2]-a[2]};
	#local cp[3];
	#local tempNormal[3]={normal[0],normal[1],normal[2]};
	cross(v1,v3,cp);
	normalizeVector(tempNormal);
	#local ang=atan2(dot(cp,tempNormal),dot(v1,v3));
	return ang;
}

// gives 3x3 determinant as follows for (A,B,C,D)
// |Dx-Ax Dy-Ay Dz-Az|
// |Dx-Bx Dy-By Dz-Bz|
// |Dx-Cx Dy-Cy Dz-Cz|
#local magicDeterminant3(struct NODE* A, struct NODE* B, struct NODE* C, struct NODE* D){
	#local out=0;

	out+=(D->x-A->x)*((D->y-B->y)*(D->z-C->z)-(D->z-B->z)*(D->y-C->y));
	out+=(D->y-A->y)*((D->z-B->z)*(D->x-C->x)-(D->x-B->x)*(D->z-C->z));
	out+=(D->z-A->z)*((D->x-B->x)*(D->y-C->y)-(D->y-B->y)*(D->x-C->x));
	
	return out;	
}

*/


#macro PlanarPointWithinTriangle (P, V1, V2, V3)
					
	#local AB = array {V2[0][0] - V1[0][0], V2[0][1] - V1[0][1]};
	#local BC = array {V3[0][0] - V2[0][0], V3[0][1] - V2[0][1]};
	#local CA = array {V1[0][0] - V3[0][0], V1[0][1] - V3[0][1]};
	#local AP = array { P[0][0] - V1[0][0],  P[0][1] - V1[0][1]};
	#local BP = array { P[0][0] - V2[0][0],  P[0][1] - V2[0][1]};
	#local CP = array { P[0][0] - V3[0][0],  P[0][1] - V3[0][1]};

	#local N1 = array {AB[1], -AB[0]};
	#local N2 = array {BC[1], -BC[0]};
	#local N3 = array {CA[1], -CA[0]};
	
	#local S1 = AP[0] * N1[0] + AP[1] * N1[1];
	#local S2 = BP[0] * N2[0] + BP[1] * N2[1];
	#local S3 = CP[0] * N3[0] + CP[1] * N3[1];

	#local Tolerance = 0.0001;

	#if ( (S1<0 & S2<0 & S3<0) | (S1<Tolerance & S2<0 & S3<0) | (S2<Tolerance & S1<0 & S3<0) | (S3<Tolerance & S1<0& S2<0))
		#local Result = 1;
		//#debug "point  is inside triangle \n"
	#else
		#local Result = 0;
		//#debug "point  is outside triangle \n"
	#end
	
	Result
	

#end

#macro Reference2DSubarray (Array, Index)
	#local Size = dimension_size (Array, 2);
	#local Temp = array[1][Size];
	
	#for (i, 0, Size-1)
		//#debug concat ("i = ", str (i, 0, 0), "\n")
		#local Temp[0][i] = Array[Index][i];
	#end
	Temp
#end


#macro DelaunayTriangulate (points, numpoints)
	#local numpoints = numpoints - 3;
	
	#debug "Finding min & max \n"
	// find min and max boundaries of point  cloud
	#local xmin = 0;
	#local xmax = 0;
	#local ymin = 0;
	#local ymax = 0;
	
	#for (i, 0, numpoints)
		#if (points[i][0] > xmax)
			#local xmax = points[i][0];
		#elseif (points[i][0] < xmin)
			#local xmin = points[i][0];
		#end
		
		#if (points[i][1] > ymax)
			#local ymax = points[i][1];
		#elseif (points[i][1] < ymin)
			#local ymin = points[i][1];
		#end
	#end // end for

	#debug "Remapping points \n"
	// remap everything (preserving the aspect ratio) to between (0,0)-(1,1)
	#local Height = ymax-ymin;
	#local Width  = xmax-xmin;
	#local d = Height; // d=largest dimension
	
	#if (Width > d)
		#local d = Width;
	#end
	
	#for (i, 0, numpoints)
		#local points[i][0] = (points[i][0] - xmin)/d;
		#local points[i][1] = (points[i][1] - ymin)/d;
	#end // end for

	#debug "Sorting points \n"
	// sort points by proximity
	#local NbinRows = ceil(pow(numpoints, 0.25));
	#local bins = array [numpoints+1];
	
	#for (i, 0, numpoints)
		#local  p = (points[i][1] * NbinRows*0.999); // bin row
		#local  q = (points[i][0] * NbinRows*0.999); // bin column
		#if (mod(p, 2))
			#local bins[i] = (p+1) * NbinRows - q;
		#else
			#local bins[i] = p * NbinRows + q + 1;
		#end
	#end
	
	#debug "Insertion sort \n"
	#for (i, 1, numpoints) // insertion sort
		#local key = bins[i];
		#local tempF = array {points[i][0], points[i][1]};
		#local  j = i - 1;
		//#debug "------------------------------------------------\n"
		//#debug concat ("      i = ", str (i, 0, 0), "\n")
		//#debug concat ("bins[i] = ", str (bins[i], 0, 0), "\n")
		//#debug concat ("bins[j] = ", str (bins[j], 0, 0), "\n")
		#if (j > 0) // wrapper
			#while ( j >= 0 & bins[j] > key)
				#local bins[j+1] = bins[j];
				#local points[j+1][0] = points[j][0];
				#local points[j+1][1] = points[j][1];
				#local j = j-1;
				
			//	#debug concat ("j = ", str (j, 0, 0), "\n")
			//	#debug concat ("bins[i] = ", str (bins[i], 0, 0), "\n")
				#if (j<0) #break #end
			#end // end while
		#end // end if (j > 0) // wrapper
		
		#local bins[j+1] = key;
		#local points[j+1][0] = tempF[0];
		#local points[j+1][1] = tempF[1];
	#end

	
	// add big triangle around our point  cloud
	//#local points = array [numpoints+3][2]; // nope - points array already exists, redeclaring it wipes it out.
	#local numpoints = numpoints + 1;
	#local points[numpoints][0] = -100;
	#local points[numpoints][1] = -100;
	#local points[numpoints+1][0] = 100;
	#local points[numpoints+1][1] = -100;
	#local points[numpoints+2][0] = 0;
	#local points[numpoints+2][1] = 100;
	#local numpoints = numpoints + 2;		//  Necessary?

	#debug "Creating data structures \n"
	// data structures required
	#local verts = array[3][3];
	#local verts[0][0] = numpoints-3;
	#local verts[0][1] = numpoints-2;
	#local verts[0][2] = numpoints-1;
	#local tris = array[3][3];
	#local tris[0][0] = -1;
	#local tris[0][1] = -1;
	#local tris[0][2] = -1;
	#declare nT = 1;
	#local triangleStack = array [numpoints-3]; // is this a big enough stack?
	#local tos = -1;
	// insert all points and triangulate one by one
	#local CNP = dimension_size (points, 1);
	#debug concat ("Current number of points = ", str (CNP, 0, 0), "\n")
	#debug concat ("numpoints = ", str (numpoints, 0, 0), "\n")
	
	#for (ii, 0, (numpoints-3))
		#debug "------------------------------------ \n"
		#debug concat ("ii = ", str (ii, 0, 0), "\n")
		// find triangle T which contains points[i]
		#local  j = nT-1; // last triangle created
		
		#debug concat ("nT = ", str (nT, 0, 0), "\n")
		#debug concat ("jstart = ", str (j, 0, 0), "\n")
			
			
		#local TestPoint = Reference2DSubarray (points, ii);
		#local _V1 = Reference2DSubarray (points, verts[j][0]);
		#local _V2 = Reference2DSubarray (points, verts[j][1]);
		#local _V3 = Reference2DSubarray (points, verts[j][2]);
		
		#local Exit = 0;
		#local Loop = 1;
		#while (Loop)
							
			//#if (PlanarPointWithinTriangle (points[ii], points[verts[j][0]], points[verts[j][1]], points[verts[j][2]]) )
			#if ( PlanarPointWithinTriangle (TestPoint, _V1, _V2, _V3) )
				#declare nT = nT + 2;
				#debug concat ("nT = ", str (nT, 0, 0), "\n")
				// delete triangle T and replace it with three sub-triangles touching P 
				//#local tris = (nT)*3; // nope - tris array already exists, redeclaring it wipes it out.
				//#local verts = (nT)*3; // nope - verts array already exists, redeclaring it wipes it out.
				// vertices of new triangles
				#local verts[nT-2][0] = ii;
				#local verts[nT-2][1] = verts[j][1];
				#local verts[nT-2][2] = verts[j][2];
				#local verts[nT-1][0] = ii;
				#local verts[nT-1][1] = verts[j][2];
				#local verts[nT-1][2] = verts[j][0];
				// update adjacencies of triangles surrounding the old triangle
				// fix adjacency of A
				#local  adj1 = tris[j][0];
				#local  adj2 = tris[j][1];
				#local  adj3 = tris[j][2];
				
				#if (adj1 >= 0)
					#for (m, 0, 3)
						#if (tris[adj1][m] = j)
							#local tris[adj1][m] = j;
							#break;
						#end // end if
					#end // end for
				#end // end if
				
				#if (adj2 >= 0)
					#for (m, 0, 3)
						#if (tris[adj2][m] = j)
							#local tris[adj2][m] = nT-2;
							#break;
						#end // end if
					#end // end for
				#end // end if
				
				#if (adj3 >= 0)
					#for (m, 0, 3)
						#if (tris[adj3][m] = j)
							#local tris[adj3][m] = nT-1;
							#break;
						#end // end if
					#end // end for
				#end // end if
				
				// adjacencies of new triangles	
				#local tris[nT-2][0] = j;
				#local tris[nT-2][1] = tris[j][1];
				#local tris[nT-2][2] = nT-1;
				#local tris[nT-1][0] = nT-2;
				#local tris[nT-1][1] = tris[j][2];
				#local tris[nT-1][2] = j;
				// replace v3 of containing triangle with P and rotate to v1
				#local verts[j][2] = verts[j][1];
				#local verts[j][1] = verts[j][0];
				#local verts[j][0] = ii;
				// replace 1st and 3rd adjacencies of containing triangle with new triangles
				#local tris[j][1] = tris[j][0];
				#local tris[j][2] = nT-2;
				#local tris[j][0] = nT-1;
				// place each triangle containing P onto a stack, if the edge opposite P has an adjacent triangle
				
				#if (tris[j][1] >= 0)
					#local triangleStack[tos+1] = j;
				#end // end if
				
				#if (tris[nT-2][1] >= 0)
					#local triangleStack[tos+1] = nT-2;
				#end // end if
				
				#if (tris[nT-1][1] >= 0)
					#local triangleStack[tos+1] = nT-1;
				#end // end if
				
				#while (tos >= 0) // looping thru the stack
					#local L = triangleStack [tos-1];
					#local v1 = array {points[verts[L][2]][0], points[verts[L][2]][1]};
					#local v2 = array {points[verts[L][1]][1], points[verts[L][1]][1]};
					#local  oppVert = -1;
					#local  oppVertID = -1;
					
					#for (k, 0, 3)
						#if ((verts[tris[L][1]][k] != verts[L][1]) & (verts[tris[L][1]][k] != verts[L][2]))
							#local oppVert = verts[tris[L][1]][k];
							#local oppVertID = k;
							#break;
						#end // end if
					#end // end for
					
				
					#local v3 = array {points[oppVert][0], points[oppVert][1]};
					#local P = array {points[ii][0], points[ii][1]};
					
					// check if P in circumcircle of triangle on top of stack
					#local cosa = ((v1[0] - v3[0])*(v2[0] - v3[0]) + (v1[1] - v3[1])*(v2[1] - v3[1]));
					#local cosb = ((v2[0] -  P[0])*(v1[0] -  P[0]) + (v2[1] -  P[1])*(v1[1] -  P[1]));
					#local sina = ((v1[0] - v3[0])*(v2[1] - v3[1]) - (v1[1] - v3[1])*(v2[0] - v3[0]));
					#local sinb = ((v2[0] -  P[0])*(v1[1] -  P[1]) - (v2[1] -  P[1])*(v1[0] -  P[0]));
					
					#if ( 	((cosa<0) & (cosb<0)) | 
							((-cosa*((v2[0]-P[0])*(v1[1]-P[1])-(v2[1]-P[1])*(v1[0]-P[0]))) > (cosb*((v1[0]-v3[0])*(v2[1]-v3[1])-(v1[1]-v3[1])*(v2[0]-v3[0]))))
					)

						// swap diagonal, and redo triangles L R A & C
						// initial state:
						#local  R = tris[L][1];
						#local  C = tris[L][2];
						#local  A = tris[R][mod ((oppVertID+2), 3)];
						// fix adjacency of A
						#if (A >= 0)
							#for (m, 0, 3)
								#if (tris[A][m]=R)
									#local tris[A][m]=L;
									#break;
								#end // end if
							#end // end for
						#end // end if
						// fix adjacency of C
						#if (C>=0)
							#for (m, 0, 3)
								#if (tris[C][m]=L)
									#local tris[C][m]=R;
									#break;
								#end // end if
							#end // end for
						#end // end if
						// fix vertices and adjacency of R
						#for (m, 0, 3)
							#if (verts[R][m]=oppVert)
								#local verts[R][mod ((m+2), 3)]=ii;
								#break;
							#end // end if
						#end // end for
					
						#for (m, 0, 3)
							#if (tris[R][m]=L)
								#local tris[R][m]=C;
								#break;
							#end // end if
						#end // end for
						
						#for (m, 0, 3)
							#if (tris[R][m]=A)
								#local tris[R][m]=L;
								#break;
							#end // end if
						#end // end for
						
						#for (m, 0, 3)
							#if (verts[R][0] != ii)
								#local  temp1      = verts[R][0];
								#local  temp2      = tris[R][0];
								#local verts[R][0] = verts[R][1];
								#local verts[R][1] = verts[R][2];
								#local verts[R][2] = temp1;
								#local tris[R][0]  = tris[R][1];
								#local tris[R][1]  = tris[R][2];
								#local tris[R][2]  = temp2;
							#end // end if
						#end // end for
						
						// fix vertices and adjacency of L
						#local verts[L][2]=oppVert;
						
						#for (m, 0, 3)
							#if (tris[L][m]=C)
								#local tris[L][m]=R;
								#break;
							#end // end if
						#end // end for
						
						#for (m, 0, 3)
							#if (tris[L][m]=R)
								#local tris[L][m]=A;
								#break;
							#end // end if
						#end // end for
						
						// add L and R to stack if they have triangles opposite P;
						#if (tris[L][1]>=0)
							#local triangleStack[tos+1]=L;
						#end // end if
						
						#if (tris[R][1]>=0)
							#local triangleStack[tos+1]=R;
						#end // end if
					#end
				#end // end if
				#break;
			#end // end while (tos >= 0)
			
			#debug concat ("jend = ", str (j, 0, 0), "\n")
			//#if (j = 0) #local Exit = 1; #end   // ####################################### HACK #######################################
			// adjust j in the direction of target point  ii
			#local AB = array {points[verts[j][1]][0] - points[verts[j][0]][0], points[verts[j][1]][1] - points[verts[j][0]][1]};
			#local BC = array {points[verts[j][2]][0] - points[verts[j][1]][0], points[verts[j][2]][1] - points[verts[j][1]][1]};
			#local CA = array {points[verts[j][0]][0] - points[verts[j][2]][0], points[verts[j][0]][1] - points[verts[j][2]][1]};
			#local AP = array {points[ii][0] - points[verts[j][0]][0], points[ii][1] - points[verts[j][0]][1]};
			#local BP = array {points[ii][0] - points[verts[j][1]][0], points[ii][1] - points[verts[j][1]][1]};
			#local CP = array {points[ii][0] - points[verts[j][2]][0], points[ii][1] - points[verts[j][2]][1]};
			#local N1 = array {AB[1], -AB[0]};
			#local N2 = array {BC[1], -BC[0]};
			#local N3 = array {CA[1], -CA[0]};
			#local S1 = AP[0]*N1[0] + AP[1]*N1[1];
			#local S2 = BP[0]*N2[0] + BP[1]*N2[1];
			#local S3 = CP[0]*N3[0] + CP[1]*N3[1];
			
			// update j
			#if ((S1>0) & (S1>=S2) & (S1>=S3))
				#local j = tris[j][0];
				#debug concat ("j1 = ", str (j, 0, 0), "\n")
			#elseif ((S2>0) & (S2>=S1) & (S2>=S3))
				#local j = max (tris[j][1], 0);
				#debug concat ("j2 = ", str (j, 0, 0), "\n")
			#elseif ((S3>0) & (S3>=S1) & (S3>=S2))
				#local j = tris[j][2];
				#debug concat ("j3 = ", str (j, 0, 0), "\n")
			#end
			
			
			//#if (Exit) #local Loop = 0; #end
		#end // end while (1)
		//#debug concat ("nT = ", str (nT, 0, 0), "\n")
	#end // end #for (ii, 0, (numpoints-3))
	
	//#debug concat ("nT = ", str (nT, 0, 0), "\n")
	
	// count how many triangles we have that dont involve supertriangle vertices
	#local nT_final = nT;
	#local renumberAdj = array [nT];
	#local deadTris = array [nT];
	
	#for (i, 0, nT)		
		#if ( (verts[i][0] >= (numpoints-3)) | (verts[i][1] >= (numpoints-3)) | (verts[i][2] >= (numpoints-3)) )
			#local deadTris[i] = 1;
			#local renumberAdj[i] = nT-(nT_final-1);
		#else 
			#local renumberAdj[i] = nT-(nT_final);
		#end
	#end
	
	// delete any triangles that contain the supertriangle vertices
	#local  verts_final = array [nT_final][3]
	#local  tris_final  = array [nT_final][3]
	#local  index = 0;
	
	#for (i, 0, nT)		
		#if ( (verts[i][0] < (numpoints-3)) & (verts[i][1] < (numpoints-3)) & (verts[i][2] < (numpoints-3)) )
			#local verts_final[index][0]  = verts[i][0];
			#local verts_final[index][1]  = verts[i][1];
			#local verts_final[index][2]  = verts[i][2];
			#local tris_final[index][0]   = (1-deadTris[tris[i][0]]) * tris[i][0]-deadTris[tris[i][0]];
			#local tris_final[index][1]   = (1-deadTris[tris[i][1]]) * tris[i][1]-deadTris[tris[i][1]];
			#local tris_final[index+1][2] = (1-deadTris[tris[i][2]]) * tris[i][2]-deadTris[tris[i][2]];
		#end
	#end
	
	#for (i, 0, nT_final)
		#if (tris_final[i][0] >= 0)
			#local tris_final[i][0] = tris_final[i][0] - renumberAdj[tris_final[i][0]];
		#end
		
		#if (tris_final[i][1] >= 0)
			#local tris_final[i][1] = tris_final[i][1] - renumberAdj[tris_final[i][1]];
		#end
		
		#if (tris_final[i][2] >= 0)
			#local tris_final[i][2] = tris_final[i][2] - renumberAdj[tris_final[i][2]];
		#end
	#end

	// undo the mapping
	#local numpoints = numpoints - 3;
	
	#if (Width > d)
		#local d = Width;
	#end
	
	#for (i, 0, numpoints)
		#local points[i][0] = points[i][0] * d + xmin;
		#local points[i][1] = points[i][1] * d + ymin;
	#end

	#if (0)
		#debug ("Triangle Vertices:\n");
		#for (i, 0, nT_final)
			printf("[%d,%d,%d]\n",verts_final[i][0],verts_final[i][1],verts_final[i][2]);
			printf("Triangle Adjacencies:\n");
		#end
		#for (i, 0, nT_final)
			printf("[%d,%d,%d]\n",tris_final[i][0],tris_final[i][1],tris_final[i][2]);
			printf("point  Cloud:\n");
		#end
		#for (i, 0, numpoints)
			printf("%d:(%f,%f)\n",i,points[i][0],points[i][1]);
		#end
	#end // end if
#end // end macro




#declare Seed = seed (123);
#declare Points1 = array [13][2]
#for (i, 0, 9)
	#local Points1[i][0] = floor (rand(Seed)*100);
	#local Points1[i][1] = floor (rand(Seed)*100);
#end

#declare NP = dimension_size (Points1, 1)-1;

DelaunayTriangulate (Points1, NP)