/*------------------------------------------------------------------------

  File        : render3d.cpp

  Description : A 3D object software-renderer

  Author      : David Tschumperlé
  
  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
  (at your option) 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.
  
  ----------------------------------------------------------------------*/

#include "../CImg.h"
// The lines below are not necessary in your own code, it simply allows 
// the source compilation with compilers that do not respect the C++ standart.
#if ( defined(_MSC_VER) && _MSC_VER<=1200 ) || defined(__DMC__)
#define std
#endif
using namespace cimg_library;

CImgl<> *points_ptr;
int primitive_sorting(const void* A,const void* B);

int main(int argc,char **argv) {

	// Read command line arguments
	//----------------------------
	cimg_usage("A real-time 3D object renderer");
	const char *pname = cimg_option("-p",(char*)NULL,"Points datafile");
	const char *fname = cimg_option("-f",(char*)NULL,"Primitives datafile");
	const char *tname = cimg_option("-t",(char*)NULL,"Texture datafile");
	const char *geom  = cimg_option("-g","640x480","Display geometry");
	const bool fullsc = cimg_option("-fullscreen",false,"Fullscreen mode");
	const char *bgcol = cimg_option("-bgcol","100,100,100","Background color");
	const char *lgcol = cimg_option("-lgcol","255,255,255","Light color");
	const char *fccol = cimg_option("-fccol","255,255,255","Faces color");
	const char *lncol = cimg_option("-lncol","255,255,255","Lines color");
	const char *ptcol = cimg_option("-ptcol","255,255,255","Points color");
	const char *bgimg = cimg_option("-bgimg",(char*)NULL,"Background image");
	const char *lgimg = cimg_option("-lgimg",(char*)NULL,"Light image");
	const char *tximg = cimg_option("-tximg",(char*)NULL,"Texture image");
	const unsigned int specul = cimg_option("-specul",50,"Light specularity (when tximg is not defined)");
	unsigned int draw_face    = cimg_option("-df",2,"Draw primitives with faces (0=off,1=flat,2=textured)");
	bool draw_line    = cimg_option("-dl",false,"Draw primitives with lines (on/off)");
	bool draw_corner  = cimg_option("-dc",false,"Draw primitives with corners (on/off)");
	const bool obj_center = cimg_option("-center",true,"Center and scale object");
	const bool add_bbox   = cimg_option("-bbox",false,"Display bounding box");
	const char *params0   = cimg_option("-params","0,0,0,2","Default parameters value (rx,ry,rz,zoom)");
	const char *dparams0  = cimg_option("-dparams","0,0,0,0","Parameters variation (drx,dry,drz,dzoom)");
	const char *file_o    = cimg_option("-o",(char*)NULL,"Output rendered images");
	const char *file_iseq = cimg_option("-iseq",(char*)NULL,"Input motion sequence");
	const char *file_oseq = cimg_option("-oseq",(char*)NULL,"Output motion sequence");
	const unsigned int nb_iter = cimg_option("-iter",-1,"Number of rendering frames");
	const bool helpflag = cimg_option("-h",false,"Display this help");
	if (helpflag) {
		fprintf(stderr,"\nWhen the render window is visible, the following keys can be pressed :\n\n");
		fprintf(stderr,"\tF1 : Swap face rendering (none,flat,textured)\n");
		fprintf(stderr,"\tF2 : Enable/Disable primitives outlines\n");
		fprintf(stderr,"\tF3 : Enable/Disable primitives corners\n");
	}

	// Init default parameters colors and images (texture,light,background)
	//---------------------------------------------------------------------

	// Parameters and parameters variations
	float alpha=0,beta=0,gamma=0,zoom=1,zoomc=zoom;
        std::sscanf(params0,"%g%*c%g%*c%g%*c%g",&alpha,&beta,&gamma,&zoom);
	alpha*=(float)cimg::PI/180.0f; beta*=(float)cimg::PI/180.0f; gamma*=(float)cimg::PI/180.0f;
	CImg<> params = CImg<>(3,3).fill(1,0,0,0,(float)std::cos(alpha),(float)std::sin(alpha),0,-(float)std::sin(alpha),(float)std::cos(alpha))*
		CImg<>(3,3).fill((float)std::cos(beta),0,-(float)std::sin(beta),0,1,0,(float)std::sin(beta),0,(float)std::cos(beta))*
		CImg<>(3,3).fill((float)std::cos(gamma),-(float)std::sin(gamma),0,(float)std::sin(gamma),(float)std::cos(gamma),0,0,0,1),
		paramsc(params);

	float dalpha=0,dbeta=0,dgamma=0,dzoom=0;
	std::sscanf(dparams0,"%g%*c%g%*c%g%*c%g",&dalpha,&dbeta,&dgamma,&dzoom);
	dalpha*=(float)cimg::PI/180.0f; dbeta*=(float)cimg::PI/180.0f; dgamma*=(float)cimg::PI/180.0f;
	CImg<> dparams = CImg<>(3,3).fill(1,0,0,0,(float)std::cos(dalpha),(float)std::sin(dalpha),0,-(float)std::sin(dalpha),(float)std::cos(dalpha))*
		CImg<>(3,3).fill((float)std::cos(dbeta),0,-(float)std::sin(dbeta),0,1,0,(float)std::sin(dbeta),0,(float)std::cos(dbeta))*
		CImg<>(3,3).fill((float)std::cos(dgamma),-(float)std::sin(dgamma),0,(float)std::sin(dgamma),(float)std::cos(dgamma),0,0,0,1);
	
	// Colors and geometry
	unsigned char background_color[3],light_color[3],point_color[3],line_color[3],face_color[3];
	unsigned int r,g,b;
	r=255,g=255,b=255; std::sscanf(bgcol,"%d%*c%d%*c%d",&r,&g,&b); background_color[0]=r; background_color[1]=g; background_color[2]=b;
	r=255,g=255,b=255; std::sscanf(lgcol,"%d%*c%d%*c%d",&r,&g,&b); light_color[0]=r; light_color[1]=g; light_color[2]=b;
	r=255,g=255,b=255; std::sscanf(fccol,"%d%*c%d%*c%d",&r,&g,&b); face_color[0]=r; face_color[1]=g; face_color[2]=b;
	r=255,g=255,b=255; std::sscanf(lncol,"%d%*c%d%*c%d",&r,&g,&b); line_color[0]=r; line_color[1]=g; line_color[2]=b;
	r=255,g=255,b=255; std::sscanf(ptcol,"%d%*c%d%*c%d",&r,&g,&b); point_color[0]=r; point_color[1]=g; point_color[2]=b;

	// Background image
	unsigned int width=640,height=480;
	std::sscanf(geom,"%d%*c%d",&width,&height);
	CImg<unsigned char> background = bgimg?CImg<unsigned char>(bgimg):CImg<unsigned char>(width,height,1,3);
	if (!bgimg) {	
		memset(background.ptr(0,0,0,0),background_color[0],background.size()/3);
		memset(background.ptr(0,0,0,1),background_color[1],background.size()/3);
		memset(background.ptr(0,0,0,2),background_color[2],background.size()/3); 
	}

	// Light image
	CImg<unsigned char> light = lgimg?CImg<unsigned char>(lgimg):CImg<unsigned char>(256,256,1,3);
	if (!lgimg) light.draw_gaussian(light.width/2.0f,light.height/2.0f,specul,light_color);

	// Texture image
	CImg<unsigned char> texture = tximg?CImg<unsigned char>(tximg):CImg<unsigned char>(256,1,1,3);
	if (!tximg) cimg_mapX(texture,x) { texture(x,0,0,0) = x; texture(x,0,0,2) = x/2; texture(x,0,0,1) = x/3; }

	// Input-output motion sequences
	CImgl<> iseq, oseq;
	if (file_iseq) iseq = CImgl<>(file_iseq);

	// Display window
	CImgDisplay disp(background,pname?cimg::basename(pname):"Demo 3D object",0,3+(fullsc?4:0)+16);
	
	// Init the 3D object
	//-------------------
	CImgl<> points, mappings;
	CImgl<unsigned int> primitives;

	if (pname) {
		// load 3d object on the disk
		points = CImgl<>(pname);
		if (fname) primitives = CImgl<unsigned int>(fname);
		else cimgl_map(points,l) primitives.insert(CImg<unsigned int>(1,1).fill(l));
		if (tname) mappings = CImgl<>(tname); else mappings = CImgl<>(points.size);
	}
	else {
          cimg::warn(true,"No input 3D object specified. ** Running demo mode. Please try '%s -h' to get rendering options **.",
                         cimg::basename(argv[0]));
		// Create the example 3D object (a torus and a textured cube)
		points.insert(CImg<>(1,3).fill(-30,-30,-30));
		points.insert(CImg<>(1,3).fill( 30,-30,-30));
		points.insert(CImg<>(1,3).fill( 30, 30,-30));
		points.insert(CImg<>(1,3).fill(-30, 30,-30));
		points.insert(CImg<>(1,3).fill(-30,-30, 30));
		points.insert(CImg<>(1,3).fill( 30,-30, 30));
		points.insert(CImg<>(1,3).fill( 30, 30, 30));
		points.insert(CImg<>(1,3).fill(-30, 30, 30));
		mappings.insert(CImg<>(1,2).fill(0,0));
		mappings.insert(CImg<>(1,2).fill(1,0));
		mappings.insert(CImg<>(1,2).fill(1,1));
		mappings.insert(CImg<>(1,2).fill(0,1));
		mappings.insert(CImg<>(1,2).fill(0,0));
		mappings.insert(CImg<>(1,2).fill(1,0));
		mappings.insert(CImg<>(1,2).fill(1,1));
		mappings.insert(CImg<>(1,2).fill(0,1));
		primitives.insert(CImg<unsigned int>(1,3).fill(0,2,1));
		primitives.insert(CImg<unsigned int>(1,3).fill(0,3,2));
		primitives.insert(CImg<unsigned int>(1,3).fill(4,5,6));
		primitives.insert(CImg<unsigned int>(1,3).fill(4,6,7));
		primitives.insert(CImg<unsigned int>(1,3).fill(0,4,7));
		primitives.insert(CImg<unsigned int>(1,3).fill(0,7,3));
		primitives.insert(CImg<unsigned int>(1,3).fill(1,2,5));
		primitives.insert(CImg<unsigned int>(1,3).fill(5,2,6));
		primitives.insert(CImg<unsigned int>(1,3).fill(0,1,4));
		primitives.insert(CImg<unsigned int>(1,3).fill(4,1,5));
		primitives.insert(CImg<unsigned int>(1,3).fill(3,7,2));
		primitives.insert(CImg<unsigned int>(1,3).fill(2,7,6));
		const float	a = 80, b = 25;
		const unsigned int na = 24, nb = 12;
		for (unsigned int v=0; v<na; v++)
			for (unsigned int u=0; u<nb; u++) {
				const float 
					alpha = (float)(u*2*cimg::PI/nb),
					beta = (float)(v*2*cimg::PI/na),
					x = (float)((a+b*std::cos(alpha))*std::cos(beta)),
					y = (float)((a+b*std::cos(alpha))*std::sin(beta)),
					z = (float)(b*std::sin(alpha));
				points.insert(CImg<>(1,3).fill(x,y,z));
				mappings.insert(CImg<>());
			}
		for (unsigned int vv=0; vv<na; vv++)
			for (unsigned int uu=0; uu<nb; uu++) {
				const int nv = (vv+1)%na, nu = (uu+1)%nb;
				primitives.insert(CImg<unsigned int>(1,3).fill(8+nb*vv+nu,8+nb*vv+uu,8+nb*nv+uu));
				primitives.insert(CImg<unsigned int>(1,3).fill(8+nb*vv+nu,8+nb*nv+uu,8+nb*nv+nu));
			}
		// Create the demo background
		CImg<> background0(background,false);
		cimg_mapXY(background0,x,y) {
                  background0(x,y,0) = (float)std::fabs(x*std::cos(16.0*y/background.height)+y*std::sin(9.0*x/background.width))/6;
                  background0(x,y,1) = (float)std::fabs(x*std::sin(18.0*y/background.height)-y*std::cos(11.0*x/background.width))/3;
                  background0(x,y,2) = (float)std::fabs(x*std::cos(13.0*y/background.height)-y*std::sin(8.0*x/background.width));
		}
		background = background0.normalize(0,255);
		texture = CImg<unsigned char>("img/logo.ppm");
		light = CImg<unsigned char>("img/light3d.ppm");
		draw_face = 2;
		draw_corner = true; point_color[0] = point_color[1] = point_color[2] = 255;
		draw_line = true; line_color[0] = 0; line_color[1] = 100; line_color[2] = 255;
    dparams = CImg<>(3,3).fill(1,0,0,0,(float)std::cos(0.05),(float)std::sin(0.05),0,-(float)std::sin(0.05),(float)std::cos(0.05));
	}
	
	// Center the object and compute the Bounding box
	if (obj_center) {
		CImg<> orig = CImg<>(1,3).fill(0);
		{ cimgl_map(points,l) orig+=points[l]; orig/=points.size; }
		{ cimgl_map(points,l) points[l]-=orig; }
		double rmax = 0;
		cimgl_map(points,l) { const double r = points[l].norm(); if (r>rmax) rmax=r; }
		rmax/=100; { cimgl_map(points,l) points[l]/=rmax; }
	}
	
	if (add_bbox) {
		float 
			xmin=points[0](0),xmax=xmin,
			ymin=points[0](1),ymax=ymin,
			zmin=points[0](2),zmax=zmin;
		cimgl_map(points,l) {
			const CImg<> &pt = points[l];
			if (pt(0)<xmin) xmin=pt(0);
			if (pt(0)>xmax) xmax=pt(0);
			if (pt(1)<ymin) ymin=pt(1);
			if (pt(1)>ymax) ymax=pt(1);
			if (pt(2)<zmin) zmin=pt(2);
			if (pt(2)>zmax) zmax=pt(2);
		}
		primitives.insert(CImg<unsigned int>(1,2).fill(points.size+0,points.size+1));
		primitives.insert(CImg<unsigned int>(1,2).fill(points.size+1,points.size+2));
		primitives.insert(CImg<unsigned int>(1,2).fill(points.size+2,points.size+3));
		primitives.insert(CImg<unsigned int>(1,2).fill(points.size+3,points.size+0));
		primitives.insert(CImg<unsigned int>(1,2).fill(points.size+4,points.size+5));
		primitives.insert(CImg<unsigned int>(1,2).fill(points.size+5,points.size+6));
		primitives.insert(CImg<unsigned int>(1,2).fill(points.size+6,points.size+7));
		primitives.insert(CImg<unsigned int>(1,2).fill(points.size+7,points.size+4));
		primitives.insert(CImg<unsigned int>(1,2).fill(points.size+0,points.size+4));
		primitives.insert(CImg<unsigned int>(1,2).fill(points.size+1,points.size+5));
		primitives.insert(CImg<unsigned int>(1,2).fill(points.size+3,points.size+7));
		primitives.insert(CImg<unsigned int>(1,2).fill(points.size+2,points.size+6));
		points.insert(CImg<>(1,3).fill(xmin,ymin,zmin));
		points.insert(CImg<>(1,3).fill(xmax,ymin,zmin));
		points.insert(CImg<>(1,3).fill(xmax,ymax,zmin));
		points.insert(CImg<>(1,3).fill(xmin,ymax,zmin));
		points.insert(CImg<>(1,3).fill(xmin,ymin,zmax));
		points.insert(CImg<>(1,3).fill(xmax,ymin,zmax));
		points.insert(CImg<>(1,3).fill(xmax,ymax,zmax));
		points.insert(CImg<>(1,3).fill(xmin,ymax,zmax));
	}

	// Compute points normals
	CImgl<> normals(points.size,CImg<>(1,3).fill(0));
	cimgl_map(primitives,l) if (primitives[l].height==3) {
		const unsigned int i0 = primitives[l](0), i1 = primitives[l](1), i2 = primitives[l](2);
		const CImg<> &p0 = points[i0], &p1 = points[i1], &p2 = points[i2];
		CImg<> p0p1(p1), p0p2(p2), p1p2(p2);
		p0p1-=p0;	p0p2-=p0;	p1p2-=p1;
		CImg<> p0p1p(p0p1);
		normals[i0] += p0p1.cross(p0p2);
		normals[i1] += p0p1p.cross(p1p2);
		normals[i2] += p0p2.cross(p1p2);
	}
	{ cimgl_map(normals,l) normals[l]/=normals[l].norm()+1e-8;}

	// Create transformed mapping coordinates to speed up rendering    
	CImgl<int> imappings;
	{ cimgl_map(mappings,l) switch (mappings[l].height) {
	case 1: imappings.insert(CImg<int>(1,2).fill((int)(mappings[l](0)*(texture.width-1)),0)); break;
	case 2: imappings.insert(CImg<int>(1,2).fill((int)(mappings[l](0)*(texture.width-1)),(int)(mappings[l](1)*(texture.height-1)))); break;
	default: imappings.insert(CImg<int>()); break;
	}
	}

	fprintf(stderr,"\n> Object initialized : %d vertices, %d primitives\n",points.size,primitives.size);

	//-------------------------
	// Start of the main loop
	//-------------------------

	CImgl<> rpoints(points), ppoints(points), rnormals(normals);
	bool stopflag = false;
	int mousexc=-1,mouseyc=-1;
	CImg<unsigned char> img,img0(background);
	const unsigned char white[3] = {255,255,255};
	
	for (unsigned int iter=0; iter<nb_iter && !disp.closed && !stopflag; iter++) {
		if (file_o) { fprintf(stderr,"> Rendering frame %d           \r",iter); fflush(stderr); }
                disp.resize();
                img0 = background.get_resize(disp);
		img=img0;
                
		// Handle input-output motion sequences
		//-------------------------------------
		if (file_oseq) {
                  img.draw_text(2,2,white,NULL,1,"Sequence frame %d",oseq.size);
                  oseq.insert(params);	
		}
		if (iseq.size) { params = iseq[0]; iseq.remove(0); }

		// Compute the rotated points/normals and the projected points
		//-------------------------------------------------------------
		const float focale = 500;
      
		cimgl_map(points,l)	{
			const float 
				X  = points[l](0),  Y  = points[l](1),  Z  = points[l](2),
				Xn = normals[l](0), Yn = normals[l](1), Zn = normals[l](2);
			CImg<> &rpt = rpoints[l], &ppt = ppoints[l], &npt = rnormals[l];
			rpt(0) = params(0,0)*X + params(1,0)*Y + params(2,0)*Z;
			rpt(1) = params(0,1)*X + params(1,1)*Y + params(2,1)*Z;
			rpt(2) = params(0,2)*X + params(1,2)*Y + params(2,2)*Z;

			npt(0) = params(0,0)*Xn + params(1,0)*Yn + params(2,0)*Zn;
			npt(1) = params(0,1)*Xn + params(1,1)*Yn + params(2,1)*Zn;
			npt(2) = params(0,2)*Xn + params(1,2)*Yn + params(2,2)*Zn;
			
			const float pZ = ppt(2) = rpt(2) + focale;
			ppt(0) = zoom*focale*rpt(0)/pZ+img.width/2;
			ppt(1) = zoom*focale*rpt(1)/pZ+img.height/2;
		}
        
		// Compute and sort visible primitives
		//-------------------------------------
		CImgl<unsigned int> visible_primitives;
		{ cimgl_map(primitives,l) {
			switch (primitives[l].height) {
			case 3: // when the primitive is a triangle
				{
					const unsigned int i0 = primitives[l](0), i1 = primitives[l](1), i2 = primitives[l](2);
					const CImg<> &pp0 = ppoints[i0];
					CImg<> pp1 = ppoints[i1], pp2 = ppoints[i2];			
					pp1-=pp0;
					pp2-=pp0;
					if (pp1[0]*pp2[1]-pp2[0]*pp1[1]<0) visible_primitives.insert(primitives[l]);
				}
				break;
			case 2: // when the primitive is a line or a point
			case 1:
				visible_primitives.insert(primitives[l]);
				break;
			}
		}}
		points_ptr = &ppoints;
                std::qsort((void*)(visible_primitives.data),visible_primitives.size,sizeof(CImg<unsigned int>),primitive_sorting);
		{ cimgl_map(rnormals,l) {
			const CImg<> &npt = rnormals[l];
			npt(0)=(npt(0)+1)*light.width/2;
			npt(1)=(npt(1)+1)*light.height/2;
		}}

		// Draw the object on the screen
		//-------------------------------
       
		{ cimgl_map(visible_primitives,l) {

			//-------------------------
			// rendering with textures
			//-------------------------
			if (draw_face==2) switch(visible_primitives[l].height) {
			case 3: // face
                          {
                            const unsigned int i0 = visible_primitives[l](0), i1 = visible_primitives[l](1), i2 = visible_primitives[l](2);
                            const CImg<>
                              &p0  = ppoints[i0], &p1  = ppoints[i1], &p2  = ppoints[i2],
                              &np0 = rnormals[i0], &np1 = rnormals[i1], &np2 = rnormals[i2];
                            const CImg<int> &m0  = imappings[i0], &m1 = imappings[i1], &m2 = imappings[i2];
                            if (!m0.height || !m1.height || !m2.height) // no mapping (only light)                              
                              img.draw_triangle((int)p0(0),(int)p0(1),(int)p1(0),(int)p1(1),(int)p2(0),(int)p2(1),
                                                light,(int)np0(0),(int)np0(1),(int)np1(0),(int)np1(1),(int)np2(0),(int)np2(1));
                            else { // texture mapping + light
                              img.
                                draw_triangle((int)p0(0),(int)p0(1),(int)p1(0),(int)p1(1),(int)p2(0),(int)p2(1),
                                              texture,m0(0),m0(1),m1(0),m1(1),m2(0),m2(1)).
                                draw_triangle((int)p0(0),(int)p0(1),(int)p1(0),(int)p1(1),(int)p2(0),(int)p2(1),
                                              light,(int)np0(0),(int)np0(1),(int)np1(0),(int)np1(1),(int)np2(0),(int)np2(1),0.5f);
                            }
                          }
                          break;
                          
			case 2: // line
				{
					const unsigned int i0 = visible_primitives[l](0), i1 = visible_primitives[l](1);
					const CImg<> &p0 = ppoints[i0], &p1 = ppoints[i1];
					const CImg<int> &m0 = imappings[i0], &m1 = imappings[i1];
					if (!m0.height || !m1.height) img.draw_line((int)p0(0),(int)p0(1),(int)p1(0),(int)p1(1),line_color);
					else img.draw_line((int)p0(0),(int)p0(1),(int)p1(0),(int)p1(1),texture,m0(0),m0(1),m1(0),m1(1));
				}
				break;

			case 1: // point
				{
					const unsigned int i0 = visible_primitives[l](0);
					const CImg<> &p0 = ppoints[i0];
					const CImg<int> &m0 = imappings[i0];
					const int X = (int)p0(0), Y = (int)p0(1);
					if (X>=0 && X<(int)img.width && Y>=0 && Y<(int)img.height)
						if (!m0.height) { img(X,Y,0) = point_color[0]; img(X,Y,1) = point_color[1]; img(X,Y,2) = point_color[2]; }
						else { img(X,Y,0) = texture(m0(0),m0(1),0,0); img(X,Y,1) = texture(m0(0),m0(1),0,1); img(X,Y,2) = texture(m0(0),m0(1),0,2); }
				}
				break;    
			}

			//-----------------------------
			// rendering with flat shading
			//-----------------------------
			if (draw_face==1) switch(visible_primitives[l].height) {
			case 3: // face
				{
					const unsigned int i0 = visible_primitives[l](0), i1 = visible_primitives[l](1), i2 = visible_primitives[l](2);
					const CImg<>
						&p0  = ppoints[i0], &p1  = ppoints[i1], &p2  = ppoints[i2],
						&np0 = rnormals[i0], &np1 = rnormals[i1], &np2 = rnormals[i2];
					float cz = -(np0(2)+np1(2)+np2(2))/3;                   
					if (cz<0) cz=0;
					unsigned char color[3];
					color[0] = (unsigned char)(face_color[0]*cz);
					color[1] = (unsigned char)(face_color[1]*cz);
					color[2] = (unsigned char)(face_color[2]*cz);
					img.draw_triangle((int)p0(0),(int)p0(1),(int)p1(0),(int)p1(1),(int)p2(0),(int)p2(1),color);
				}
				break;

			case 2: // line
				{
					const unsigned int i0 = visible_primitives[l](0), i1 = visible_primitives[l](1);
					const CImg<> &p0 = ppoints[i0], &p1 = ppoints[i1];
					img.draw_line((int)p0(0),(int)p0(1),(int)p1(0),(int)p1(1),line_color);                    
				}
				break;

			case 1: // point
				{
					const unsigned int i0 = visible_primitives[l](0);
					const CImg<> &p0 = ppoints[i0];                    
					const int X = (int)p0(0), Y = (int)p0(1);
					if (X>=0 && X<(int)img.width && Y>=0 && Y<(int)img.height) {
						img(X,Y,0) = point_color[0]; img(X,Y,1) = point_color[1]; img(X,Y,2) = point_color[2]; 
					}                        
				}
				break;    
			}

			//-------------------------
			// rendering with outlines
			//-------------------------
			if (draw_line) switch(visible_primitives[l].height) {
			case 3: // face
				{
					const unsigned int i0 = visible_primitives[l](0), i1 = visible_primitives[l](1), i2 = visible_primitives[l](2);
					const CImg<> &p0  = ppoints[i0], &p1  = ppoints[i1], &p2  = ppoints[i2];
					img.draw_line((int)p0(0),(int)p0(1),(int)p1(0),(int)p1(1),line_color); 
					img.draw_line((int)p0(0),(int)p0(1),(int)p2(0),(int)p2(1),line_color);
					img.draw_line((int)p2(0),(int)p2(1),(int)p1(0),(int)p1(1),line_color);
				}
				break;

			case 2: // line
				{
					const unsigned int i0 = visible_primitives[l](0), i1 = visible_primitives[l](1);
					const CImg<> &p0 = ppoints[i0], &p1 = ppoints[i1];
					img.draw_line((int)p0(0),(int)p0(1),(int)p1(0),(int)p1(1),line_color);                    
				}
				break;

			case 1: // point
				{
					const unsigned int i0 = visible_primitives[l](0);
					const CImg<> &p0 = ppoints[i0];                    
					const int X = (int)p0(0), Y = (int)p0(1);
					if (X>=0 && X<(int)img.width && Y>=0 && Y<(int)img.height) {
						img(X,Y,0) = point_color[0]; img(X,Y,1) = point_color[1]; img(X,Y,2) = point_color[2]; 
					}                        
				}
				break;    
			}

			//------------------------------
			// rendering with corner points
			//------------------------------
			if (draw_corner) switch(visible_primitives[l].height) {
			case 3: // face
				{
					const unsigned int i0 = visible_primitives[l](0), i1 = visible_primitives[l](1), i2 = visible_primitives[l](2);
					const CImg<> &p0  = ppoints[i0], &p1  = ppoints[i1], &p2  = ppoints[i2];
					int X = (int)p0(0), Y = (int)p0(1);
					if (X>=0 && X<(int)img.width && Y>=0 && Y<(int)img.height) { img(X,Y,0) = point_color[0]; img(X,Y,1) = point_color[1]; img(X,Y,2) = point_color[2]; }
					X = (int)p1(0); Y = (int)p1(1);
					if (X>=0 && X<(int)img.width && Y>=0 && Y<(int)img.height) { img(X,Y,0) = point_color[0]; img(X,Y,1) = point_color[1]; img(X,Y,2) = point_color[2]; }
					X = (int)p2(0); Y = (int)p2(1);
					if (X>=0 && X<(int)img.width && Y>=0 && Y<(int)img.height) { img(X,Y,0) = point_color[0]; img(X,Y,1) = point_color[1]; img(X,Y,2) = point_color[2]; }
				}
				break;

			case 2: // line
				{
					const unsigned int i0 = visible_primitives[l](0), i1 = visible_primitives[l](1);
					const CImg<> &p0 = ppoints[i0], &p1 = ppoints[i1];
					int X = (int)p0(0), Y = (int)p0(1);
					if (X>=0 && X<(int)img.width && Y>=0 && Y<(int)img.height) { img(X,Y,0) = point_color[0]; img(X,Y,1) = point_color[1]; img(X,Y,2) = point_color[2]; }
					X = (int)p1(0); Y = (int)p1(1);
					if (X>=0 && X<(int)img.width && Y>=0 && Y<(int)img.height) { img(X,Y,0) = point_color[0]; img(X,Y,1) = point_color[1]; img(X,Y,2) = point_color[2]; }
				}
				break;

			case 1: // point
				{
					const unsigned int i0 = visible_primitives[l](0);
					const CImg<> &p0 = ppoints[i0];                    
					const int X = (int)p0(0), Y = (int)p0(1);
					if (X>=0 && X<(int)img.width && Y>=0 && Y<(int)img.height) { img(X,Y,0) = point_color[0]; img(X,Y,1) = point_color[1]; img(X,Y,2) = point_color[2]; }
				}
				break;    
			}
		}}

		// Handle mouse and keys for 3D animation
		//----------------------------------------
		if (disp.key) {
			switch(disp.key) {
			case cimg::keyQ:  stopflag = true; break;
			case cimg::keyF1: draw_face=(draw_face+1)%3; break;
			case cimg::keyF2: draw_line=!draw_line; break;
			case cimg::keyF3: draw_corner=!draw_corner; break;
			}
			disp.key=0;
		}

		if (!disp.button) {
			params *= dparams;
			zoom += dzoom;
			if (mousexc>=0 || mouseyc>=0) mousexc=mouseyc=-1;
		}
		else {
			if (mousexc<0 && mouseyc<0) {
				mousexc = disp.mousex;
				mouseyc = disp.mousey;
				paramsc = params;
				zoomc = zoom;
			}
			float u = (float)disp.mousex - mousexc,
				v = (float)disp.mousey - mouseyc,
                          n = (float)std::sqrt(u*u+v*v);               

			if (disp.button&1 && n>=1) {
				u/=n; v/=n;
				n*=(float)(4*cimg::PI/std::sqrt((double)(disp.width*disp.width+disp.height*disp.height)));
				const CImg<> P = CImg<>(3,3).fill(u,-v,0,v,u,0,0,0,1),
					R = CImg<>(3,3).fill((float)std::cos(n),0,-(float)std::sin(n),0,1,0,(float)std::sin(n),0,(float)std::cos(n));
				params = P*R*P.get_transpose()*paramsc;
			}
			if (disp.button&2) zoom = zoomc+v*5/disp.height;
		}

		// End of the main loop
		disp.display(img).wait(20);
		if (file_o) img.save(file_o,iter);
	
	}

	if (file_oseq) { fprintf(stderr,"> Saving motion sequence '%s'\n",file_oseq); oseq.save(file_oseq); }
  return 0;
}

// function that compares the z between two primitives
//-----------------------------------------------------
int primitive_sorting(const void* A,const void* B) {
	const CImg<unsigned int>& a = *((CImg<unsigned int>*)A), b = *((CImg<unsigned int>*)B);
	const CImgl<>& points = *points_ptr;
	const float
		za = a.height==3?(points[a(0)](2)+points[a(1)](2)+points[a(2)](2))/3:(a.height==2?(points[a(0)](2)+points[a(1)](2))/2:(a.height==1?points[a(0)](2):0)),
		zb = b.height==3?(points[b(0)](2)+points[b(1)](2)+points[b(2)](2))/3:(b.height==2?(points[b(0)](2)+points[b(1)](2))/2:(b.height==1?points[b(0)](2):0));
	if (zb>za) return 1;
	return -1;
}
