#include <iostream>
#include <CImg.h>
#include <complex>
//#include "CPDSI_functions.h"
using namespace std;
using namespace cimg_library;


CImg <double> Pasar_A_Grises(CImg <double> img)
{
	CImg<double>gris;
	if(img.dimv() == 3){
		gris=(img.get_channel(0)+img.get_channel(1)+img.get_channel(2))/3.0;
        }else{
		gris=img;
	}
	return gris;
}


CImg <double> Bordes(CImg <double> img){
   CImg <double> X(3,3,1,1), Y(3,3,1,1), XY(3,3,1,1), YX(3,3,1,1), imgX, imgY, imgXY, imgYX, img2(img.dimx(),img.dimy(),1,1,0);
   Y(0,0)=-1;
   Y(0,1)=-2;
   Y(0,2)=-1;
   Y(1,0)=0;
   Y(1,1)=0;
   Y(1,2)=0;	
   Y(2,1)=1;
   Y(2,0)=2;
   Y(2,2)=1;

   X(0,0)=-1;
   X(0,1)=0;
   X(0,2)=1;
   X(1,0)=-2;
   X(1,1)=0;
   X(1,2)=2;	
   X(2,0)=-1;
   X(2,1)=0;
   X(2,2)=1;

   imgX=img.get_convolve(X);
   imgY=img.get_convolve(Y);
   
   img2=imgX.abs()+imgY.abs(); //sumamos las dos imagenes utilizando el valor absoluto en vez de la suma de los cuadrados y la raiz
   img2.normalize(0,1);
   img2.threshold(0.2);	
   return img2;
}

CImg <double> Correlacion(CImg <double>imgA,CImg <double> imgB)
{
	complex <double> complejo(0.,0.), I(0.,1.), complejoFA(0.,0.), complejoFB(0.,0.);
	double valorFA0, valorFA1, valorFB0, valorFB1;

	int M=imgA.dimx()+imgB.dimx();
	int N=imgA.dimy()+imgB.dimy();
	imgB.rotate(180);
	imgA.resize(M,N,1,1,0);
	imgB.resize(M,N,1,1,0);
	CImgl <> FA = imgA.get_FFT();
	CImgl <> FB = imgB.get_FFT();
	
	CImg <double> real1(M,N,1,1,0);
	CImg <double> imag1(M,N,1,1,0);
	CImg <double> correlacion;
	for(int col=0; col<M;col++){
 	   for(int fil=0; fil<N; fil++){
	
		valorFA0 = FA[0](col,fil);
		valorFA1 = FA[1](col,fil);
		complejoFA = valorFA0 + I* valorFA1;
		
		valorFB0 = FB[0](col,fil);
		valorFB1 = FB[1](col,fil);
		complejoFB = valorFB0 + I* valorFB1;

		complejo = complejoFA * complejoFB;
			
	        real1(col,fil) = real(complejo);
		imag1(col,fil) = imag(complejo);
	   }
	}	  

	CImgl <double> corrF(real1,imag1);
	correlacion = corrF.get_FFT(true)[0];
	
	return correlacion;

}

CImg <double> Ajustar_Tono(CImg <double> cachitoA, CImg <double> cachitoB, CImg <double> imgB)
{
	const CImgStats statsA0(cachitoA.get_channel(0));
	const CImgStats statsB0(cachitoB.get_channel(0));
	double tono0 = statsA0.mean - statsB0.mean;
	
	const CImgStats statsA1(cachitoA.get_channel(1));
	const CImgStats statsB1(cachitoB.get_channel(1));
	double tono1 = statsA1.mean - statsB1.mean;
	
	const CImgStats statsA2(cachitoA.get_channel(2));
	const CImgStats statsB2(cachitoB.get_channel(2));
	double tono2 = statsA2.mean - statsB2.mean;

	for(int i=0; i<imgB.dimy();i++){
		for(int j=0; j<imgB.dimx(); j++){
			imgB(j,i,0,0)=imgB(j,i,0,0) + tono0;
			imgB(j,i,0,1)=imgB(j,i,0,1) + tono1;
			imgB(j,i,0,2)=imgB(j,i,0,2) + tono2;
		}
	}
	return imgB;

}

CImg <double> Unir(CImg <double> imgA, CImg <double> imgB, int desp_horizontal, int desp_vertical)
{
	CImg <double> cachitoA, cachitoB;
	int x1,x2,x3,y1,y2,y3;
	x1 = desp_horizontal;
	x2 = imgA.dimx()-1;
	x3 = desp_horizontal + imgB.dimx();
	
	if(desp_vertical>=0){
		y1=desp_vertical;
		y2=imgA.dimy()-1;
		if (imgA.dimy()>imgB.dimy() + desp_vertical){
			y3=imgA.dimy()-1;
		}else{
			y3=imgB.dimy() + desp_vertical;
		}
		
	}else{
		y1=abs(desp_vertical);
		y2=imgB.dimy()-1;
		if (imgB.dimy()>imgA.dimy() + abs(desp_vertical)){
			y3=imgB.dimy()-1;
		}else{
			y3=imgA.dimy() + abs(desp_vertical);
		}
		
	}
			
	CImg <double> unidas(x3 + 1, y3 + 1, 1, 3, 0);
	if(desp_vertical>=0){
		for(int i=0; i<imgA.dimy();i++){
			for(int j=0; j<imgA.dimx(); j++){
				for(int v=0;v<3;v++){
					unidas(j,i,0,v)=imgA(j,i,0,v);
				}
	 		}
		}
		cachitoA = imgA.get_crop(x1+10,y1+10,x2-10,y2-10);
		cachitoB = imgB.get_crop(10,10,x2-x1-10,y2-y1-10);
		imgB = Ajustar_Tono(cachitoA, cachitoB, imgB);

		//parte que se superpone
		for(int i=y1+10; i<=y2-10; i++){
			for(int j=x1+10; j<=x2-10; j++){
				for(int v=0;v<3;v++){
					unidas(j,i,0,v)=(imgB(j-x1,i-y1,0,v)+unidas(j,i,0,v))/2; //para corregir el tono de la imagenB, que quede mas parecido al de la imagenA
				}
			}
		}
		
		//parte que no se superpone a la derecha	
		for(int i=y1; i<=y2-10; i++){
			for(int j=x2+1-10; j<=x3; j++){
				for(int v=0;v<3;v++){
					unidas(j,i,0,v)=imgB(j-x1,i-y1,0,v);
				}
			}
		}

		//parte que no se superpone abajo
		for(int i=y2+1-10; i<=y3; i++){
			for(int j=x1; j<=x3;j++){
				for(int v=0;v<3;v++){
					unidas(j,i,0,v)=imgB(j-x1,i-y1,0,v);
				}
			}
		}
	}else{
		for(int i=0; i<imgA.dimy();i++){
			for(int j=0; j<imgA.dimx(); j++){
				for(int v=0;v<3;v++){
					unidas(j,i+y1,0,v)=imgA(j,i,0,v);	 
				}
	 		}
		}

		cachitoA = imgA.get_crop(x1+10,10,x2-10,y2-y1-10);
		cachitoB = imgB.get_crop(10,y1+10,x2-x1-10,y2-10);
		imgB = Ajustar_Tono(cachitoA, cachitoB, imgB);
			
		//parte que no se superpone de arriba
		for(int i=0; i<y1+10; i++){
			for(int j=x1; j<=x3;j++){
				for(int v=0;v<3;v++){
					unidas(j,i,0,v)=imgB(j-x1,i,0,v);
				}
			}
		}
		
		//parte que se superpone
		for(int i=y1+10; i<=y2-10; i++){
			for(int j=x1+10; j<=x2-10; j++){
				for(int v=0;v<3;v++){
					unidas(j,i,0,v)=(imgB(j-x1,i,0,v) + unidas(j,i,0,v))/2; //para corregir el tono de la imagenB, que quede mas parecido al de la imagenA
				}
			}
		}
		
		//parte que no se superpone a la derecha
		for(int i=y1; i<=y2; i++){
			for(int j=x2+1-10; j<=x3; j++){
				for(int v=0; v<3; v++){
					unidas(j,i,0,v)=imgB(j-x1,i,0,v);
				}
			}
		}
		

		
	}
	return unidas;



}


int main(){
	
	CImg<double> imgA("lata1.jpg"),imgB("lata2.jpg"), grisA, grisB, correlacion,unidas;
	int col, fil;
	const double negro[3]={0.0,0.0,0.0};
	
	///Pasamos las imágenes en color a tonos de gris
	grisA = Pasar_A_Grises(imgA);
	grisB = Pasar_A_Grises(imgB);
	
	///Realizamos la detección de bordes mediante máscaras de Sobel
	grisA = Bordes(grisA);
	grisB = Bordes(grisB);

	///Le sacamos a la imagen 10 pixeles de los bordes parra evitar bordes no deseados (los del contorno, en algunos casos)
	grisA.draw_rectangle(0,0,grisA.dimx()-1,9,negro);
	grisA.draw_rectangle(0,0,9,grisA.dimy()-1,negro);
	grisA.draw_rectangle(0,grisA.dimy()-11,grisA.dimx()-1,grisA.dimy()-1,negro);
	grisA.draw_rectangle(grisA.dimx()-11,0,grisA.dimx()-1,grisA.dimy()-1,negro);
	
	grisB.draw_rectangle(0,0,grisB.dimx()-1,9,negro);
	grisB.draw_rectangle(0,0,9,grisB.dimy()-1,negro);
	grisB.draw_rectangle(0,grisB.dimy()-11,grisB.dimx()-1,grisB.dimy()-1,negro);
	grisB.draw_rectangle(grisB.dimx()-11,0,grisB.dimx()-1,grisB.dimy()-1,negro);

	///Calculamos la correlación entre las imágenes de bordes
	correlacion = Correlacion(grisA, grisB);

	///Calculamos el punto donde se da el valor máximo de la correlación y calculamos los desplazamientos relativos entre las imágenes
	const CImgStats stats(correlacion);
	int max_x =stats.xmax;
	int max_y =stats.ymax; 
	int desp_h = max_x -grisB.dimx()+1; 
	int desp_v = max_y -grisB.dimy()+1; 

	///Unimos las dos imágenes con la información de los desplazamientos
	unidas = Unir(imgA,imgB,desp_h,desp_v);
	

	CImgDisplay ventana1(imgA, "imagen A");
	CImgDisplay ventana2(imgB, "imagen B");
	//CImgDisplay ventana3(grisA, "imagen A");
	//CImgDisplay ventana4(grisB, "imagen B");
	CImgDisplay ventana5(unidas, "unidas");
	//CImgDisplay ventana6(correlacion, "correlacion");
	
	unidas.normalize(0,255);
	unidas.save("unida.jpg");

	ventana1.wait();
	while ((!ventana1.is_closed)
		&& (ventana1.key != cimg::keyQ))
	{}
}