// This program is free software: you can use, modify and/or redistribute it
// under the terms of the simplified BSD Lcense. You should have received a
// copy of this license along this program. If not, see
// <http://www.opensource.org/licenses/bsd-license.html>.
//
// Copyright (C) 2017, Javier Sánchez Pérez <jsanchez@ulpgc.es>
// All rights reserved.


#include "crop_and_zoom.h"
#include "direct_method/transformation.h"
#include "sutherland_hodgman.h"


#include <stdio.h>
#include <stack>

/**
  *
  * Computes the inscribed rectangle by updating the corners
  * according to the projection of the image corners
  *
**/
void fast_rectangle(
  float *H,    //set of transformations to be cropped
  int nx,      //number of columns
  int ny,      //number of rows
  int nframes, //number of frames
  int nparams, //type of transformation
  float &x1,   //x coordinate of top-left corner
  float &y1,   //y coordinate of top-left corner
  float &x2,   //x coordinate of bottom-right corner
  float &y2    //y coordinate of bottom-right corner
)
{  
  float *H_1=new float[nparams];

  x1=y1=0; x2=nx-1; y2=ny-1;

  //find the inside common rectangle
  for(int i=0;i<nframes;i++)
  { 
    inverse_transform(&(H[i*nparams]), H_1, nparams);
    float xp, yp;

    //project the image corners and store the limits
    project(0,0,H_1,xp,yp,nparams);
    if(xp>x1) x1=xp;
    if(yp>y1) y1=yp;

    project(nx-1,0,H_1,xp,yp,nparams);
    if(xp<x2) x2=xp;
    if(yp>y1) y1=yp;

    project(0,ny-1,H_1,xp,yp,nparams);
    if(xp>x1) x1=xp;
    if(yp<y2) y2=yp;

    project(nx-1,ny-1,H_1,xp,yp,nparams);
    if(xp<x2) x2=xp;
    if(yp<y2) y2=yp;
  }

  delete []H_1;
}



/**
  *
  * Find the intersecting polygon to all the homographies 
  *
**/
poly find_intersecting_polygon(
  float *H,    //set of transformations to be cropped
  int nx,      //number of columns
  int ny,      //number of rows
  int nframes, //number of frames
  int nparams  //type of transformation
)
{
  float *H_1=new float[nparams];

  poly clipper=poly_new();
  clipper->len=4;
  clipper->alloc=0;
  clipper->v=new vec_t[4];
  clipper->v[0].x=0   ; clipper->v[0].y=0; 
  clipper->v[1].x=nx-1; clipper->v[1].y=0;
  clipper->v[2].x=nx-1; clipper->v[2].y=ny-1; 
  clipper->v[3].x=0   ; clipper->v[3].y=ny-1;

  //find the inside common rectangle
  for(int i=0;i<nframes;i++)
  {
    inverse_transform(&(H[i*nparams]), H_1, nparams);

    //project the image corners and store the limits
    vec_t s[4]; 
    project(0, 0, H_1,s[0].x,s[0].y,nparams);
    project(nx-1,0, H_1,s[1].x,s[1].y,nparams);
    project(nx-1, ny-1,H_1,s[2].x,s[2].y,nparams);
    project(0,ny-1,H_1,s[3].x,s[3].y,nparams);

    poly_t subject = {4, 0, s};

    //compute the intersection of both polygons using sutherland-hodgman method
    poly res=poly_clip(&subject, clipper);

    poly_free(clipper);
    clipper=res;
  }

  delete []H_1;

  return clipper;
}


/**
  *
  * Create a binary image with 1 inside the common
  * region and 0 in the empty regions 
  *
**/
void create_binary_image(
  poly polygon, //intersecting polygon
  int  *I,      //output binary image with polygon mask
  int  nx,      //number of columns
  int  ny       //number of rows
)
{
  for(int i=0; i<nx*ny; i++) I[i]=0;

  float ymin=ny, ymax=0;

  //find the y-limits of the polygon
  for(int i=0; i<polygon->len; i++)
  {
    if(ymin>polygon->v[i].y) ymin=polygon->v[i].y;
    if(ymax<polygon->v[i].y) ymax=polygon->v[i].y;
  }

  //intersect with the polygon all the rows in between and fill
  for(int y=(int)(ymin+0.5); y<(int)(ymax+0.5); y++)
  {
    float xmin=-1, xmax=-1;

    int i=0;
    int j=polygon->len-1;
    int intersections=0;

    //find the intersection of the current row with the polygon
    //since the polygon is convex we need to find two contacts
    while(i<polygon->len && intersections<2)
    {
      //intersection test
      if((polygon->v[j].y>=y && polygon->v[i].y<=y)||
         (polygon->v[j].y<=y && polygon->v[i].y>=y))
      {
        float x1=polygon->v[j].x;
        float y1=polygon->v[j].y;
        float x2=polygon->v[i].x;
        float y2=polygon->v[i].y;

        //if line is coincident, fill it completely
        if((int)(y1+0.5)==(int)(y2+0.5))
        {
          xmin=x1; xmax=x2;
          intersections=2;
        }
        else if(intersections==0)
          xmin=(y-y1)*(x2-x1)/(y2-y1)+x1;
        else
          xmax=(y-y1)*(x2-x1)/(y2-y1)+x1;

        intersections++;
      }

      j=i;
      i++;
    }

    //swap x limits if not ordered
    if(xmin>xmax)
    {
      float tmp=xmin;
      xmin=xmax;
      xmax=tmp;
    }

    //fill the row 
    if(xmin>=0 && xmin<nx && xmax>=0 && xmax<nx)
      for(int x=(int)xmin; x<=(int)xmax; x++) I[y*nx+x]=255;
  }
}


/**
  *
  * Find the rectangle from the binary image
  * Method of maximum rectangle beneath a histogram
  *
**/
void find_rectangle(
  int *I,    //binary image with polygon mask
  int nx,    //number of columns of the image
  int ny,    //number of rows of the image
  float &x1, //computed x coordinate of top-left corner of the rectangle
  float &y1, //computed y coordinate of top-left corner of the rectangle
  float &x2, //computed x coordinate of bottom-right corner of the rectangle
  float &y2  //computed y coordinate of bottom-right corner of the rectangle
)
{
  //initialize histogram array to 0
  int *H=new int[nx*ny]();

  //create the histogram array
  for(int i=1; i<ny; i++)
    for(int j=0; j<nx; j++)
      if(I[i*nx+j])
        H[i*nx+j]=1+H[(i-1)*nx+j];

  long int max_area = 0;

  //find the maximum area axis-parallel rectangle 
  for(int i=1; i<ny; i++)
  { 
    int *L=new int[nx]();
    int *R=new int[nx]();
    std::stack<int> stack; 

    //find the left limit for this position of the histogram
    for(int j=0; j<nx; j++)
    {
      while (!stack.empty() && H[i*nx+j]<=H[i*nx+stack.top()])
        stack.pop();

      int t;
      if(stack.empty())
        t=0;
      else
        t=stack.top();

      L[j]=t;
      stack.push(j);
    }

    //empty stack 
    std::stack<int> empty;
    std::swap(stack, empty);

    //find the right limit for this position of the histogram
    for (int j=nx-1; j>=0; j--)
    {
      while(!stack.empty() && H[i*nx+j]<=H[i*nx+stack.top()])
        stack.pop();

      int t;
      if(stack.empty())
        t=nx-1;
      else
        t=stack.top();

      R[j]=t;
      stack.push(j);
    }

    //calculate final area in the row and find the maximum
    for (int j=0; j<nx; j++)
    {
      long int A=H[i*nx+j]*(R[j]-L[j]);
      if (A>max_area)
      {
        max_area=A;
        x1=L[j]+1;
        y1=i-H[i*nx+j];
        x2=R[j];
        y2=i;
      }
    }
    delete []L;
    delete []R;
  }

  delete []H;
}


/**
  *
  * Find the optimal inscribed rectangle 
  *
**/
void accurate_rectangle(
  float *H,    //set of transformations to be cropped
  int nx,      //number of columns
  int ny,      //number of rows
  int nframes, //number of frames
  int nparams, //type of transformation
  float &x1,   //x coordinate of top-left corner
  float &y1,   //y coordinate of top-left corner
  float &x2,   //x coordinate of bottom-right corner
  float &y2,   //y coordinate of bottom-right corner
  int verbose  //verbose mode
)
{
  if(verbose)
    printf("   Calculating intersection polygon\n");
  poly polygon=find_intersecting_polygon(H,nx,ny,nframes,nparams);

  int *I=new int[nx*ny];

  if(verbose)
    printf("   Creating binary image with polygon mask\n");
  create_binary_image(polygon, I, nx, ny);

  if(verbose)
    printf("   Computing the inscribed rectangle\n");
  find_rectangle(I, nx, ny, x1, y1, x2, y2);  

  poly_free(polygon);
  delete []I;
}

/**
  *
  * Computes the transformation neccessary to remove the 
  * empty regions (similarity: translate + scale)
  *
**/
void apply_crop_zoom_transform(
  float *H,
  int nx,
  int ny,
  int nframes,
  int nparams,
  float xx1,
  float yy1,
  float xx2, 
  float yy2,
  int verbose
)
{
  //round off
  int x1=(int) (xx1+.5);
  int y1=(int) (yy1+.5);
  int x2=(int) (xx2-.1);
  int y2=(int) (yy2-.1);

  //compute necessary zoom and translation to apply to the transformations
  float s, s1=(float) (y2-y1)/ny;
  float s2=(float) (x2-x1)/nx;

  if(s1<s2)
    s=s1;
  else
    s=s2;

  //translate the center of the rectangle to the origin
  float dx=(float) (x1+x2)/2.;
  float dy=(float) (y1+y2)/2.;
  float dxx=(float) nx/2.;
  float dyy=(float) ny/2.;

  //create a similarity transform for composition
  float matrix[9]={s,0.,dx-s*dxx,0.,s,dy-s*dyy,0.,0.,1.};
  float *T=new float[nparams];

  matrix2params(matrix,T,nparams);

  if(verbose)
  {
    float xp1, xp2, yp1, yp2;
    project(0,0,T,xp1,yp1,nparams);
    project(nx,ny,T,xp2,yp2,nparams);

    printf("   -Frame size: (0, 0)-(%d, %d)\n", nx, ny);
    printf("   -Crop rectangle: (%d, %d)-(%d, %d)\n", x1, y1, x2, y2);
    printf("   -Real rectangle: (%f, %f)-(%f, %f)\n", xp1, yp1, xp2, yp2);
    printf("   -Transformation: zoom=%f, translation=(%f, %f)\n", s, dx, dy);
  }

  for(int i=0;i<nframes;i++) 
    //modify the smoothing transformations
    compose_transform(&(H[i*nparams]), T, &(H[i*nparams]),nparams);

  delete []T;
}


/**
  *
  * Function for dealing with empty regions after stabilization
  *
**/
void crop_and_zoom(
  float *H,    //set of transformations to be cropped
  int nx,      //number of columns
  int ny,      //number of rows
  int nframes, //number of frames
  int nparams, //type of transformation
  int type,    //type of filling strategy
  int verbose  //verbose mode
)
{
  //test if the transformations has enough degrees of freedom
  if(nparams<=3)
  {
    if(verbose)
      printf("   Cannot apply Crop&Zoom to translation "
             "or Euclidean transforms\n");
    return;
  }

  //coordinates of the cropping rectangle
  float x1, y1, x2, y2;

  //calculating the cropping rectangle following one strategy
  switch(type)
  {
    case NO_POSTPROCESS:
      if(verbose)
        printf("  No post-processing:\n");
      break;
    case FAST_CROP_ZOOM:
      if(verbose)
        printf("  Fast Crop and Zoom:\n");
      fast_rectangle(H, nx, ny, nframes, nparams, x1, y1, x2, y2);
      //compute the new transformation for cropping and zooming
      apply_crop_zoom_transform(H, nx, ny, nframes, nparams, x1, y1, x2, y2, verbose);
      break;
    case CROP_ZOOM:
      if(verbose)
        printf("  Crop and Zoom:\n");
      accurate_rectangle(H, nx, ny, nframes, nparams, x1, y1, x2, y2, verbose);
      //compute the new transformation for cropping and zooming
      apply_crop_zoom_transform(H, nx, ny, nframes, nparams, x1, y1, x2, y2, verbose);
      break;
    case NEIGHBOR_INPAINT:
      if(verbose)
        printf("  Neighbor inpainting (not implemented yet):\n");
      break;
  }
}