#include <CImg.h>
#include <iostream>
#include <fstream>
#include <string>
#include <cmath>
#include <math.h>
#include <stdlib.h>
#include <time.h>
//#include "../CPDSI_functions/CPDSI_functions.h"
using namespace cimg_library;
using namespace std;



// ----- Sobel ---------------------------------------------------------------------------
void Sobel(CImg <double> imagen, CImg <double> &gradiente){
   CImg <float> hx (3,3,1,1,0), hy (3,3,1,1,0);
   
   hy(0,0)=-1;		hx(0,0)=-1;	
   hy(0,1)=-2;   	hx(1,0)=-2;
   hy(0,2)=-1;  	hx(2,0)=-1;
   hy(2,0)=1;	 	hx(0,2)=1;	
   hy(2,1)=2;    	hx(1,2)=2;
   hy(2,2)=1;    	hx(2,2)=1;

   imagen.channel(0);

   CImg <float> img_filx = imagen.get_convolve(hx).get_abs(); 
   CImg <float> img_fily = imagen.get_convolve(hy).get_abs();
   
// armamos el grandiente
   CImg <float> suma (img_filx.dimx(), img_filx.dimy(),1,1,0);
	
   cimg_forXY(imagen,x,y){
		suma(x,y)=img_filx(x,y) + img_fily(x,y); 
	   	(suma(x,y)<127)? suma(x,y)=0 :suma(x,y)=255;
    }
    gradiente=suma;
}

// ----- Busca los Vertices de la tapa del  CD   ------------------------------------------------
CImg <double> BuscaVertice(CImg <double> img){
  int x;
  CImg <double> resultado(img);
  resultado.fill(0);
  for (int w=0; w<img.dimy(); w++){
	if (!(w==0)){	
	x=0;
	for (int y=w; y>=0;y--){
	//cout<<" "<<x<<" "<<y<<" ";
		if (img(x,y)==255) { resultado(x,y)=255;
					//cout<<"encontro un punto"<<endl;
				     return(resultado);
					}
		x++;	
	}}
  }
}

//----Calcula los coeficientes de las imagenes--------------
void coefi(CImg <> A,CImg <> b,CImg <> &r){
//  Ax=b tenemos A y b para solucionar el sistema calculamos la inversa 
//  Ainv*b=x 

	CImg <> Ainv(A.dimx(),A.dimy(),1,1,0); 
	Ainv=invert(A,true);	 // el true es porque utiliza la factorización LU

	double aux=0.0;
	
	for(int y=0;y<Ainv.dimy();y++){ // recorre por fila
		for(int x=0;x<Ainv.dimx();x++){
			aux += Ainv(x,y)*b(0,x);		// y,x
		}
		r(0,y)=aux;
		aux=0.0;
	}
	//r.display();  // este vendria a ser x los coeficientes que necesitamos para la 
 }


//-----Sist de Resolución de Ecuaciones--------------------
void SolverWarp(double &c1,double &c2,double &c3,double &c4,double &c5,double &c6,double &c7,double &c8,int coord[][2]){

	// x, y: columna, fila
	int x1=coord[0][0];  int y1=coord[0][1];
	int x2=coord[1][0];  int y2=coord[1][1];
	int x3=coord[2][0]; int y3=coord[2][1];
	int x4=coord[3][0]; int y4=coord[3][1];
	CImg <> AI(8,8,1,1,0);
	CImg <> b(1,8,1,1,0);
	CImg <> r(1,8,1,1,0);
	
	// Siempre buscamos "deformar" la tapa del CD a un cuadrado de 250*250 
	b(0,0) = y1 ; 
	b(0,1) = x1 ;
	b(0,2) = y2 ;
	b(0,3) = x2 ;
	b(0,4) = y3 ;
	b(0,5) = x3 ;
	b(0,6) = y4 ;
	b(0,7) = x4 ;
	
	x1=0;  y1=0;
	x2=0;  y2=249;
	x3=249; y3=249;
	x4=249; y4=0;

	AI(0,0)=y1;	//cord[0][0];
	AI(0,1)=x1;	//cord[0][1];
	AI(0,2)=x1*y1;	//cord[0][0]*cord[0][1];
	AI(0,3)=1;
	AI(0,4)=0;
	AI(0,5)=0;
	AI(0,6)=0;
	AI(0,7)=0;
	
	AI(1,0)=0;
	AI(1,1)=0;
	AI(1,2)=0;
	AI(1,3)=0;
	AI(1,4)=y1;	//cord[0][0];
	AI(1,5)=x1;	//cord[0][1];
	AI(1,6)=x1*y1;	//cord[0][0]*cord[0][1];
	AI(1,7)=1;
	 
	//---------------------------------
	
	// los valores de la dos siguientes filas son con el 2º vértice x2 e y2
	
	AI(2,0)=y2;	//cord[1][0];
	AI(2,1)=x2;	//cord[1][1];
	AI(2,2)=x2*y2;  //cord[1][0]*cord[1][1];
	AI(2,3)=1;
	AI(2,4)=0;
	AI(2,5)=0;
	AI(2,6)=0;
	AI(2,7)=0;
	
	AI(3,0)=0;
	AI(3,1)=0;
	AI(3,2)=0;
	AI(3,3)=0;
	AI(3,4)=y2;	//cord[1][0];
	AI(3,5)=x2; 	//cord[1][1];
	AI(3,6)=x2*y2;	//cord[1][0]*cord[1][1];
	AI(3,7)=1;
	
	//---------------------------------
	
	// los valores de la dos siguientes filas son con el 3º vértice x3 e y3
	
	AI(4,0)=y3;	//cord[2][0];
	AI(4,1)=x3;	//cord[2][1];
	AI(4,2)=x3*y3;	//cord[2][0]*cord[2][1];
	AI(4,3)=1;
	AI(4,4)=0;
	AI(4,5)=0;
	AI(4,6)=0;
	AI(4,7)=0;
	
	AI(5,0)=0;
	AI(5,1)=0;
	AI(5,2)=0;
	AI(5,3)=0;
	AI(5,4)=y3;	//cord[2][0];
	AI(5,5)=x3;	//cord[2][1];
	AI(5,6)=x3*y3;	//cord[2][0]*cord[2][1];
	AI(5,7)=1;
	
	//---------------------------------
	
	// los valores de la dos siguientes filas son con el 4º vértice x4 e y4
	
	AI(6,0)=y4;	//cord[3][0];
	AI(6,1)=x4;	//cord[3][1];
	AI(6,2)=x4*y4;	//cord[3][0]*cord[3][1];
	AI(6,3)=1;
	AI(6,4)=0;
	AI(6,5)=0;
	AI(6,6)=0;
	AI(6,7)=0;
	
	AI(7,0)=0;
	AI(7,1)=0;
	AI(7,2)=0;
	AI(7,3)=0;
	AI(7,4)=y4;	//cord[3][0];
	AI(7,5)=x4;	//cord[3][1];
	AI(7,6)=x4*y4;	//cord[3][0]*cord[3][1];
	AI(7,7)=1;
	
	AI.transpose();
	
	coefi(AI,b,r);
	c1=r(0,0);
	c2=r(0,1);
	c3=r(0,2);
	c4=r(0,3);
	c5=r(0,4);
	c6=r(0,5);
	c7=r(0,6);
	c8=r(0,7);
	
}//----fin  SolverWarp ----------------------------------
void SolverProy(double &c1,double &c2,double &c3,double &c4,double &c5,double &c6,double &c7,double &c8,int coord[][2]){
// x, y: columna, fila
	int x1=coord[0][0];  int y1=coord[0][1];
	int x2=coord[1][0];  int y2=coord[1][1];
	int x3=coord[2][0]; int y3=coord[2][1];
	int x4=coord[3][0]; int y4=coord[3][1];
	CImg <> AI(8,8,1,1,0);
	CImg <> b(1,8,1,1,0);
	CImg <> r(1,8,1,1,0);
	
	
	// Siempre buscamos "deformar" la tapa del CD a un cuadrado de 250*250 
	b(0,0) = y1 ; 
	b(0,1) = x1 ;
	b(0,2) = y2 ;
	b(0,3) = x2 ;
	b(0,4) = y3 ;
	b(0,5) = x3 ;
	b(0,6) = y4 ;
	b(0,7) = x4 ;
	
	x1=0;  y1=0;
	x2=0;  y2=249;
	x3=249; y3=249;
	x4=249; y4=0;

	AI(0,0)=x1;	//cord[0][0];
	AI(0,1)=0;	//cord[0][1];
	AI(0,2)=-x1*b(0,1);	//cord[0][0]*cord[0][1];
	AI(0,3)=y1;
	AI(0,4)=0;
	AI(0,5)=-y1*b(0,1);
	AI(0,6)=1;
	AI(0,7)=0;
	
	AI(1,0)=0;
	AI(1,1)=x1;
	AI(1,2)=-x1*b(0,0);
	AI(1,3)=0;
	AI(1,4)=y1;	//cord[0][0];
	AI(1,5)=-y1*b(0,0);	//cord[0][1];
	AI(1,6)=0;	//cord[0][0]*cord[0][1];
	AI(1,7)=1;
	 
	//---------------------------------
	
	// los valores de la dos siguientes filas son con el 2º vértice x2 e y2
	
	AI(2,0)=x2;	//cord[1][0];
	AI(2,1)=0;	//cord[1][1];
	AI(2,2)=-x2*b(0,3);  //cord[1][0]*cord[1][1];
	AI(2,3)=y2;
	AI(2,4)=0;
	AI(2,5)=-y2*b(0,3);
	AI(2,6)=1;
	AI(2,7)=0;
	
	AI(3,0)=0;
	AI(3,1)=x2;
	AI(3,2)=-x2*b(0,2);
	AI(3,3)=0;
	AI(3,4)=y2;	//cord[1][0];
	AI(3,5)=-y2*b(0,2); 	//cord[1][1];
	AI(3,6)=0;	//cord[1][0]*cord[1][1];
	AI(3,7)=1;
	
	//---------------------------------
	
	// los valores de la dos siguientes filas son con el 3º vértice x3 e y3
	
	AI(4,0)=x3;	//cord[2][0];
	AI(4,1)=0;	//cord[2][1];
	AI(4,2)=-x3*b(0,5);	//cord[2][0]*cord[2][1];
	AI(4,3)=y3;
	AI(4,4)=0;
	AI(4,5)=-y3*b(0,5);
	AI(4,6)=1;
	AI(4,7)=0;
	
	AI(5,0)=0;
	AI(5,1)=x3;
	AI(5,2)=-x3*b(0,4);
	AI(5,3)=0;
	AI(5,4)=y3;	//cord[2][0];
	AI(5,5)=-y3*b(0,4);	//cord[2][1];
	AI(5,6)=0;	//cord[2][0]*cord[2][1];
	AI(5,7)=1;
	
	//---------------------------------
	
	// los valores de la dos siguientes filas son con el 4º vértice x4 e y4
	
	AI(6,0)=x4;	//cord[3][0];
	AI(6,1)=0;	//cord[3][1];
	AI(6,2)=-x4*b(0,7);	//cord[3][0]*cord[3][1];
	AI(6,3)=y4;
	AI(6,4)=0;
	AI(6,5)=-y4*b(0,7);
	AI(6,6)=1;
	AI(6,7)=0;
	
	AI(7,0)=0;
	AI(7,1)=x4;
	AI(7,2)=-x4*b(0,6);
	AI(7,3)=0;
	AI(7,4)=y4;	//cord[3][0];
	AI(7,5)=-y4*b(0,6);	//cord[3][1];
	AI(7,6)=0;	//cord[3][0]*cord[3][1];
	AI(7,7)=1;
	
	AI.transpose();
	
	coefi(AI,b,r);
	c1=r(0,0);
	c2=r(0,1);
	c3=r(0,2);
	c4=r(0,3);
	c5=r(0,4);
	c6=r(0,5);
	c7=r(0,6);
	c8=r(0,7);
	
}//----fin  SolverProy ----------------------------------
//-------Warping----------------------------------
void Warping (CImg <double> &img, CImg <double> imgOrig,int coord[][2]){
// img en una imagen de 250 por 250 inicializada en 0
// imgOrig es la img orginal 
   
	double  c1,c2,c3,c4,c5,c6,c7,c8;

	SolverWarp(c1,c2,c3,c4,c5,c6,c7,c8,coord);
	//SolverProy(c1,c2,c3,c4,c5,c6,c7,c8,coord);
	CImg<double> img_0(250,250,1,3,0), imgpo(250,250,1,3,0); 
	ofstream arch;

	// implementamos el mapeo inverso:
	double xt,yt;

	cimg_forXY(img,x,y){  // recorremos la img vacía
	  yt=double(c1*y+c2*x+c3*x*y+c4);
	  xt=double(c5*y+c6*x+c7*x*y+c8);
    	   
 	  // yt=double((c1*(double)x+c4*(double)y+c7)/(c3*(double)x+c6*(double)y+1));
	  // xt=double((c2*(double)x+c5*(double)y+c8)/(c3*(double)x+c6*(double)y+1));
	   // aca es donde viene lo de la interpolación..

	   int ixt,iyt;
	   
	   ixt=int(xt);
	   iyt=int(yt);

           arch.open ("salida.txt", ios::app);
	   arch <<" el punto ixt iyt :"<<ixt<<"  "<<iyt<< "\n"; 
           
           if ((xt>0) & (xt<249) & (yt>0) & (yt<249)){	           	
	   	// interpolación de orden cero
	   	img_0(x,y,0,0)=imgOrig(ixt,iyt,0,0);
	   	img_0(x,y,0,1)=imgOrig(ixt,iyt,0,1);
	   	img_0(x,y,0,2)=imgOrig(ixt,iyt,0,2);
	   }

	   //interpolación de primero orden (bi-lineal)
	   imgpo(x,y,0,0)=imgOrig.get_channel(0).linear_pix4d(xt,yt);
	   imgpo(x,y,0,1)=imgOrig.get_channel(1).linear_pix4d(xt,yt);
	   imgpo(x,y,0,2)=imgOrig.get_channel(2).linear_pix4d(xt,yt);

	   //interpolación de segund orden (bi-cubica)
	   img(x,y,0,0)=imgOrig.get_channel(0).cubic_pix2d(xt,yt);
	   img(x,y,0,1)=imgOrig.get_channel(1).cubic_pix2d(xt,yt);
	   img(x,y,0,2)=imgOrig.get_channel(2).cubic_pix2d(xt,yt);
    	   
	   arch.close(); 
	// Tenemos que poner el if porque mirando en el archivo encontramos valores negativos --> por eso hacemos la conversión ixt=int(abs(xt));
	// Y también hay coordenadas que se pasan por ejemplo  el punto ixt iyt :260  233


}//del for
	img.display();
	   img_0.save("w_o_0.jpg");     // interpolación de orden cero
	   imgpo.save("w_p_o.jpg");    //interpolación de primero orden (bi-lineal)
	   img.save("w_s_o.jpg");      //interpolación de segund orden (bi-cubica)
  
}//del void
//-----------------MAIN--------------------------------------------------------------
int main(){
	
	//CImg <double> img("fotos/ejemplo6.jpg");
	//CImg <double> img("fotos/ejemplo8.JPG");
	CImg <double> img("fotos/ejemplo5.JPG");
	img.resize(250,250);  // Hacemos un resize para operar con menos cant de filas y columnas
	CImg <double> img1(img);

	// Descomposición HSI
	img.RGBtoHSI();
	CImg <double> patron = img.get_channel(2); 		// intensidad
	patron.normalize(0,255);

	CImg <double> patronS = img.get_channel(1); 		// saturación
	patronS.normalize(0,255);

	double control;
	if (patron(1,1)>127) control=1.0;
		else control=0.0;
	int cont=0;
	cimg_forXY(patron,x,y){
		if (patronS(x,y)<45) { patron(x,y)=0; patronS(x,y)=0; cont++; }
	}
	
	double porc=(double)cont/62500;
	int T=(int)(img1.mean()*abs((control-porc)));
	
	//--------------histograma  de I --------------------------

	CImg <double> h_i=patron.get_histogram(256,0,255);

	const unsigned char rojo[]  = { 255,0,0 },verde[] = { 0,255,0 }, azul[]  = { 0,0,255 };

	CImg <unsigned char> hi(300,300,1,3,255);
	hi.draw_graph(h_i,rojo, 2, 2500, 0 , 1 ); 
	hi.draw_axis(0,255,2500,0,azul,-60,-60,-1,-1,0.7f);
	//hi.display();
	

	//fin de histograma------------------------------------

	// --- Detección de Bordes ----------------------------
	
	CImg <double> grad(patron);
	//Sobel(patron,grad);
	
	CImg <double> I(patron),res(patron),aux(patron);	
	cimg_forXY(I,x,y){
		(I(x,y)>T)?I(x,y)=255:I(x,y)=0;
	}
	Sobel(I,grad);
	grad.display("gradiente");
	grad.save("gradiente.jpg");
	//grad.display("gradiente umbral");
	CImg <double> h(3,3,1,1,0);
	cimg_forXY(grad,x,y){
			h(0,0)=grad(x-1,y-1);
			h(0,1)=grad(x-1,y,0);
			h(0,2)=grad(x-1,y+1);
			h(1,0)=grad(x,y-1);
			h(1,1)=grad(x,y);
			h(1,2)=grad(x,y+1);
			h(2,0)=grad(x+1,y-1);
			h(2,1)=grad(x+1,y);
			h(2,2)=grad(x+1,y+1);
			grad(x,y)=h.median();
		}
	//--------Buscamos los vertices-----------------------
	grad.save("gradientefil.jpg");
	int cord[4][2];

	res=BuscaVertice(grad);

	CImg <double> estad = res.get_stats();
        float omax = estad(0,5);
        int y=(int) (omax/res.dimx()); 
        int x=(int) omax-(y*res.dimx());
        cord[0][0]=x; cord[0][1]=y;
	// cout<<" "<<x<<" "<<y<<" ";
	
	for (int i=1; i<=3;i++){
	grad.rotate(90);     //rotamos la imagen para buscar los puntos siguientes
	//grad.display("imagen rotada");
	aux=BuscaVertice(grad);
	aux.rotate(-90*i);

	CImg <double> estad = aux.get_stats();
        float omax = estad(0,5);
        int y=(int)(omax/aux.dimx()); 
        int x=(int)omax-(y*aux.dimx());
        cord[i][0]=x; cord[i][1]=y;
        //cout<<" "<<x<<" "<<y<<" ";

	res+=aux;            //sumamos el punto encontrado a los puntos
	// res.display("vertice");
	}
	grad.rotate(90);
	//fin de buscar puntos------------------------------
	res.save("vertices.jpg");
	// ----- Ampliación del Warping

        //-------- Warping y tiempos de ejecución--------------------
	time_t inicio, fin;
	struct tm *tpoInicio, *tpoFin;
	// Empezamos a medir el tiempo
	inicio=time(NULL);
	CImg <double> final(250,250,1,3,0);

	Warping(final,img1,cord);
	
	fin=time(NULL);
	cout<<"Tiempo trancurrido (seg) del Warping: " << difftime(fin,inicio)<< endl;
	
	final.display("salio?");

	// -- Visualización--------------------------------
	
	CImgList<double> visu(3);
	
	visu[0]=patron;
	visu[1]=grad;
	visu[2]=final;
		
	CImgDisplay disp (visu, "imagen cd,  gradiente y vertices",0);
	
	while (!disp.is_closed){
	// si hay un resize del display, la img se redimensiona  
        if (disp.is_resized) disp.resize(disp.window_dimx(),disp.window_dimx()/2).display(visu);	
		}
	
	
	
}