#include <list>
#include <math.h>
#include <Magick++.h>
#include <fstream.h>
#include <stdlib.h>

#define max(a,b) (a>b?a:b)
#define min(a,b) (a<b?a:b)

#define FEQ 0.00001
using namespace std;
using namespace Magick;



bool traversed[512][512];

enum ANIMTYPE {LINEAR,PARAB};
int morphtype=LINEAR;

double fn(int ctr,int nframes,int type)
{


  if(type==LINEAR)
    return 1.0*(ctr-1)/(nframes-1);
  if(type==PARAB)
    return 1.0*(ctr-1)*(ctr-1)/(nframes-1);
  
  return 0;

}

bool feq(double a,double b)
{
  if(fabs(a-b)<=FEQ)
    return true;

  return false;

}

class leda_point
{
public:
 
  double x,y;

public:
  
  leda_point(double a,double b)
  {
    x=a; y=b;
  }

  leda_point()
  {
  }


  
  double length(leda_point a)
  {
    return sqrt((x-a.x)*(x-a.x)+(y-a.y)*(y-a.y));
  }

  double xcoord()
  {
    return x;
  }
  double ycoord()
  {
    return y;

  }
  

  bool operator ==(leda_point a)
  {
    if(feq(length(a),0))
      return true;
    return false;
    
  }
  
};

class CTriangle
{
public:
  
  leda_point p1,p2,p3;



  CTriangle()
  {
  }
  
  
  

  CTriangle(leda_point a,leda_point b,leda_point c)
  {
    p1=a;
    p2=b;
    p3=c;
  }
  

  void bayertoxy(double a,double b,double c,int& x,int& y)
  {

     x= (int)(a*p1.xcoord()+b*p2.xcoord()+c*p3.xcoord());
     
     y= (int)(a*p1.ycoord()+b*p2.ycoord()+c*p3.ycoord());
  }
  
  void coord(leda_point p,double& a,double& b,double& c)
  {

    double d=(p1.x-p3.x)*(p2.y-p3.y)-(p2.x-p3.x)*(p1.y-p3.y);
    d=-d;

    a=((p2.x-p3.x)*(p.y-p3.y)-(p.x-p3.x)*(p2.y-p3.y))/d;
    b=((p1.x-p3.x)*(p.y-p3.y)-(p.x-p3.x)*(p1.y-p3.y))/d;
    b=-b;
    c=1-a-b;
    
    
  }

 
    
  

  

  void return_xrange(int y,double& minx, double& maxx)
  {
    double m1,c1,m2,c2,m3,c3;



    minx=9999;
    maxx=0;
    double x;
    // first line

    if(y>=min(p1.y,p2.y) && y<=max(p1.y,p2.y))
      {
	
	x=return_x(p1,p2,y);

	// lines on line
	if(feq(x,-1))
	  {
	    minx=min(p1.x,p2.x);
	    maxx=max(p1.x,p2.x);
	    return;
	  }

	minx=min(minx,x);
	maxx=max(maxx,x);
	
      }


    if(y>=min(p2.y,p3.y) && y<=max(p2.y,p3.y))
      {
	
	x=return_x(p2,p3,y);

	// lines on line
	if(feq(x,-1))
	  {
	    minx=min(p2.x,p3.x);
	    maxx=max(p2.x,p3.x);
	    return;
	  }

	minx=min(minx,x);
	maxx=max(maxx,x);
	
      }

    if(y>=min(p1.y,p3.y) && y<=max(p1.y,p3.y))
      {
	
	x=return_x(p1,p3,y);

	// lines on line
	if(feq(x,-1))
	  {
	    minx=min(p1.x,p3.x);
	    maxx=max(p1.x,p3.x);

	    return;
	  }

	minx=min(minx,x);
	maxx=max(maxx,x);
	
      }

  }


  double return_x(leda_point p1,leda_point p2,double y)
  {
    double m1,c1;

    if(p1.x==p2.x)
      return p1.x;
    
    m1=(p1.y-p2.y)/(p1.x-p2.x);
    c1=p1.y-m1*p1.x;
    
    if(feq(m1,0)) return -1;
    
    return (y-c1)/m1;
  
  }


  
}; // end of triangle class



int Triangulate(CTriangle *arr_Triangles,leda_point *arr_points,int no_of_points);
leda_point interpolate(leda_point p1,leda_point p2,double lambda);
double interpolate(double p1,double p2,double lambda);
void find_fpoints();
void find_triangles(leda_point* fpoints,CTriangle* triangles);
void morph(double t);
void filename(int ctr,char*s);


//___Make morphing code

double spixels[512][512];
double mpixels[512][512];
double epixels[512][512];

leda_point sfpoints[100];
leda_point mfpoints[100];
leda_point efpoints[100];


int tsfpoints[200][3]; // contains the index values for the triangles


CTriangle str[200];
CTriangle mtr[200];
CTriangle etr[200];
char filenames[50][15];
char outputfile[20];

int nfpoints,ntr;
int rows,cols;
int nframes;




  
  
void main(int argc,char* argv[])
{

  char startfile[80],endfile[80],featurefile[80],middlefile[80];
  Image source,target;


  int i;

  if(argc<4)
    {
      cout << "\nMorph 1.0\nSyntax : morph startfile endfile no_of_frames featuredatafile[optional] outputfile[optional] delay" << endl << flush;
      exit(0);

    }
  
     

    strcpy(startfile,argv[1]);
    strcpy(endfile,argv[2]);

    nframes=atoi(argv[3]);
    if(nframes<2 || nframes>30)
      {

	cout << "Problem with no of frames. exiting...";
	exit(0);
      }

    //__generate feature file
    if(argc<5)
      {
	cout << "Calling Point Mapper...";

	char buffer[80];
	strcpy(buffer,"java pointMapper ");
	strcat(buffer,startfile);
	strcat(buffer," ");
	strcat(buffer,endfile);
	
	system(buffer);
	
	strcpy(featurefile,"points.dat");
	
      }
    else
      strcpy(featurefile,argv[4]);


    if(argc>=6)
      {
	strcpy(outputfile,argv[5]);
      }
   //___Read the images and put them in array as grayscale
    
   
  source = Image(startfile);
  cols = source.baseColumns();  
  rows = source.baseRows();
  source.classType(DirectClass);
  source.quantizeColorSpace(GRAYColorspace);

  target = Image(endfile);
  int tcols=target.baseColumns();  
  int trows=target.baseRows();
  target.classType(DirectClass);
  target.quantizeColorSpace(GRAYColorspace);

  if(tcols!=cols || trows!=rows)
    {
      cout << "Images don't match in size" << flush;
      exit(0);

    }

  
  Pixels viewsrc(source);
  PixelPacket *pixelsrc = viewsrc.get(0,0,cols,rows);

  Pixels viewtgt(target);
  PixelPacket *pixeltgt = viewtgt.get(0,0,cols,rows);
  
  for(int x=0;x<rows;x++)
    for(int y=0;y<cols;y++)
      {
	spixels[x][y] = ((ColorGray)(*pixelsrc)).shade();	
	pixelsrc++;
	epixels[x][y] = ((ColorGray)(*pixeltgt)).shade();
	pixeltgt++;
	

	traversed[x][y]=false;
	
      }
  pixelsrc-= cols*rows;

  //___Read the feature points

    ifstream inFile(featurefile);
    int xs,ys,xt,yt;
    
    inFile>>nfpoints;

    if(nfpoints<5)
      {
	cout << "Number of feature points < 5.Please specify more...";
	exit(0);
      }
    
    nfpoints+= 4;
  
    
    efpoints[0]=sfpoints[0] = leda_point(0,0);
    efpoints[1]=sfpoints[1] = leda_point(0,rows);
    efpoints[2]=sfpoints[2] = leda_point(cols,0);
    efpoints[3]=sfpoints[3] = leda_point(cols,rows);
    
    
    for( i=4;i<nfpoints;i++)
      {
	inFile>>xs>>ys>>xt>>yt;
	
	sfpoints[i] = leda_point(xs,ys);
	efpoints[i] = leda_point(xt,yt);      
      }

    inFile.close();
    
    //___Now morph the pictures__________

    double t;
    double inc=1.0/nframes;
    double par=inc;
    
    int ctr = 1;
    
    //__find Traingles and feature points for triangles for start image

    
  
    ntr=Triangulate(str,sfpoints,nfpoints);

    /*    cout << "no of triangles:  " << ntr;

    for(i=0;i<ntr;i++)
      {
	cout << endl << tsfpoints[i][0] << " " << tsfpoints[i][1] << " " << tsfpoints[i][2] << " " ;
	
      }

    cout << "no of f points : " << nfpoints << flush;
    
  for(i=0;i<nfpoints;i++)
      {
	cout << endl << sfpoints[i].xcoord() << " " << 	 sfpoints[i].ycoord() << " " << 	efpoints[i].xcoord() << " " << efpoints[i].ycoord() << " " ;
      }
  
    

    cout << flush;
    */
    //find_fpoints(); // not needed anymore

    //__generate triangles for the source image
    find_triangles(sfpoints,str);    
     //___Triangulate the end image
     find_triangles(efpoints,etr);
     
    //___write first file as source file
     filename(ctr,middlefile);

     
     source.write(middlefile);
    cout<<"\ncreating frame "<<ctr<<"......done" << flush;

    ctr++;
    
    while(ctr<nframes)
      {


	//par=0.5;
    
 
	//	cout << "\ncalling morph function..." << flush;

	cout<<"\ncreating frame no. "<<ctr << flush;
	
	morph(fn(ctr,nframes,morphtype));
    
	//cout << "\nmorph fn. called..." << flush;
    
	// write morphed image
	Image average(Geometry(cols,rows),"white");
	average.classType(DirectClass);
	Pixels newpix(average);
	PixelPacket *pixval=newpix.get(0,0,cols,rows);
    
	for(int i=0;i<rows;i++)
	  for(int j=0;j<cols;j++)
	    {
	      (*pixval) = ColorGray(mpixels[i][j]);
	      pixval++;
	      traversed[i][j]=false;
	    }
    
	newpix.sync();
	filename(ctr,middlefile);

	//cout << "writing image to : " << middlefile << endl;
    
	average.write(middlefile);
	//cout<< "\n......done" << flush;
	ctr++;
	//par += inc;
	cout << "\ndone..." <<flush;
      }

    //__write dest file as final image
    filename(ctr,middlefile);
    target.write(middlefile);
    cout<<"\ncreating frame no. "<<ctr << flush;





    if(argc>=6)
      {
	cout << "\nWriting animation  to : " << outputfile  << flush;
	
	
	
	char bigbuffer[1000];
	strcpy(bigbuffer,"convert -loop 0 -delay ");
	if(argc==7)
	  strcat(bigbuffer,argv[6]);
	else strcat(bigbuffer,"20");

	strcat(bigbuffer," ");
	
	for(i=0;i<nframes;i++)
	  {
	    strcat(bigbuffer,filenames[i]);
	    strcat(bigbuffer," ");
	    
	  }
	strcat(bigbuffer,outputfile);
	
	system(bigbuffer);
      }
    
 }





int Triangulate(CTriangle *arr_Triangles,leda_point *arr_points,int no_of_points)
{

  // generate a file for the feature points

  int i,count;
  fstream ffile("ftemp",ios::out);
  ffile << 2 <<endl;
  ffile << no_of_points << endl;
  for(i=0;i<no_of_points;i++)
      ffile << arr_points[i].x << " " << arr_points[i].y << endl;
  ffile.flush();
  ffile.close();

    // now we have the feature point file

  system("qdelaunay i QJ Pp < ftemp > ttemp");

  fstream tfile("ttemp",ios::in);
  tfile >> count;
  for(i=0;i<count;i++)
       tfile >> tsfpoints[i][0] >>   tsfpoints[i][1] >>  tsfpoints[i][2];
  tfile.close();

  // got the triangles ready
  return count;
}




inline leda_point interpolate(leda_point p1,leda_point p2,double lambda)
  // gives average of two points weighted by lambda
{
  int x,y;
  
  x = interpolate(p1.xcoord(),p2.xcoord(),lambda);
  y = interpolate(p1.ycoord(),p2.ycoord(),lambda);

  return leda_point(x,y);
}

inline double interpolate(double p1,double p2,double lambda)
{
  return (1.0-lambda)*p1+lambda*p2;
}


//___ find traiangle for any arbit. feature point array
void find_triangles(leda_point* fpoints,CTriangle* triangles)
{

  int i,j;

  leda_point p1,p2,p3;
  for(i=0;i<ntr;i++)
        triangles[i]=CTriangle(fpoints[tsfpoints[i][0]],fpoints[tsfpoints[i][1]],fpoints[tsfpoints[i][2]]);

}



void morph(double t)
{
  
  //find feature points of t frame

  int i,j,k;
  for(i=0;i<nfpoints;i++)
      mfpoints[i]=interpolate(sfpoints[i],efpoints[i],t);
  
  //form triangles of t frame
 
  
  find_triangles(mfpoints,mtr);
  
  //__for each triangle do morphing

  double minx,maxx,miny,maxy;

  double a,b,c;
  int xs,ys,xt,yt;

 
  for(i=0;i<ntr;i++)
    {
      /*
      minx=min(mtr[i].p1.xcoord(),mtr[i].p2.xcoord());
      minx=min(minx,mtr[i].p3.xcoord());
      maxx=max(mtr[i].p1.xcoord(),mtr[i].p2.xcoord());
      maxx=max(maxx,mtr[i].p3.xcoord());
      
      /*
      cout << "\ncoord of triangle : " << i;
      cout << "\n x1 : " << mtr[i].p1.xcoord() << " y1 : " << mtr[i].p1.ycoord();
      cout << "\n x2 : " << mtr[i].p2.xcoord() << " y2 : " << mtr[i].p2.ycoord();
      cout << "\n x3 : " << mtr[i].p3.xcoord() << " y3 : " << mtr[i].p3.ycoord();
      */
      
      
      miny=min(mtr[i].p1.ycoord(),mtr[i].p2.ycoord());
      miny=min(miny,mtr[i].p3.ycoord());
           
      maxy=max(mtr[i].p1.ycoord(),mtr[i].p2.ycoord());
      maxy=max(maxy,mtr[i].p3.ycoord());
      

      

      // scan conversion for triangles
      for(k=miny;k<=maxy;k++)
	{
	  mtr[i].return_xrange(k,minx,maxx);

	  //cout << "\nTriangle : " << i << "  ycord : " << k << "  minx : " << minx << "  maxx : " << maxx << flush;

	  for(j=minx;j<=maxx;j++)
	  {

	    //if(traversed[k][j])
	    // continue;
	    
	      mtr[i].coord(leda_point(j,k),a,b,c);
	    
	      /*	      if(a < 0 || a > 1 || b < 0 || b > 1 || c < 0 || c > 1)
		{
		cout << "\nsomething is wrong with the bayer coords..." << flush;
		cout << a << " " << b << "  " << c << " cord :  " << j << "," << k <<flush;
		}*/
		//got the coordinates
	    

	      /*    if(a<0) a=0;
	      if(a>1) a=1;
	      
	      if(b<0) b=0;
	      if(b>1) b=1;
	      
	      if(c<0) c=0;
	      if(c>1) c=1;
	      */
	     str[i].bayertoxy(a,b,c,xs,ys);
	    etr[i].bayertoxy(a,b,c,xt,yt);

	    double tt=interpolate(spixels[ys][xs],epixels[yt][xt],t);

	    if(tt>1)
	      cout << "\n%?&^^$3smth. is wrong with color values " << flush;
	    mpixels[k][j]=tt;

	  // if(feq(tt,0))
	  //    cout << "smth. is wrong here !!" << flush;

	    traversed[k][j]=true;
	  }
	}
      
    } // this triangle is done


}



void filename(int ctr,char*s)
{
  strcpy(s,"output_");
  int l=strlen(s);
  s[l]=64+ctr;

  s[l+1]=0;
  
  strcat(s,".jpg");

  strcpy(filenames[ctr-1],s);
  
  
}