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

#include <time.h>
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <gd.h>

int greyscale(gdImagePtr img, int c) {  
  int r = gdImageRed(img, c);
  int g = gdImageGreen(img, c);
  int b = gdImageBlue(img, c);

  return (int)((r+g+b)/3.+0.5);
}

void gradient(int *imageGradient, gdImagePtr img, int sx, int sy) { 
  // static int imageGradient[sx][sy]; 
  int up, down, right, left, center; 
  int i, j; 
  
  // 0,0

  right = greyscale(img, gdImageGetPixel(img, 1, 0));
  down = greyscale(img, gdImageGetPixel(img, 0, 1));
  center = greyscale(img, gdImageGetPixel(img, 0, 0));
  
  imageGradient[0 * sy + 0] = (int)(1./sqrt(3.)*sqrt((right-center)*(right-center)+(center-down)*(center-down))+0.5);
  
  // sx-1,0

  left  = greyscale(img, gdImageGetPixel(img, sx-2, 0));
  down = greyscale(img, gdImageGetPixel(img, sx-1, 1));
  center = greyscale(img, gdImageGetPixel(img, sx-1, 0));
  imageGradient[(sx-1) * sy + 0] = (int)(1./sqrt(3.)*sqrt((center-left)*(center-left)+(center-down)*(center-down))+0.5);

  // 0,sy-1

  right = greyscale(img, gdImageGetPixel(img, 1, sy-1));
  up = greyscale(img, gdImageGetPixel(img, 0, sy-2));
  center = greyscale(img, gdImageGetPixel(img, 0, sy-1));
  imageGradient[0 * sy + (sy-1)] = (int)(1./sqrt(3.)*sqrt((right-center)*(right-center)+(up-center)*(up-center))+0.5);

  // sx-1,sy-1

  left = greyscale(img, gdImageGetPixel(img, sx-2, sy-1));
  up = greyscale(img, gdImageGetPixel(img, sx-1, sy-2));
  center = greyscale(img, gdImageGetPixel(img, sx-1, sy-1));
  imageGradient[(sx-1) * sy + (sy-1)] = (int)(1./sqrt(3.)*sqrt((center-left)*(center-left)+(up-center)*(up-center))+0.5);


  // up line

  for (i=1;i<sx-1;++i) {
    center = greyscale(img, gdImageGetPixel(img, i, 0));
    right = greyscale(img, gdImageGetPixel(img, i+1, 0));
    left = greyscale(img, gdImageGetPixel(img, i-1, 0));
    down = greyscale(img, gdImageGetPixel(img, i, 1));

    imageGradient[i * sy + 0] = (int)(1./sqrt(3.)*sqrt((right-left)*(right-left)/
    4.+(center-down)*(center-down))+0.5);
  }

  // down line

  for (i=1;i<sx-1;++i) {
    center = greyscale(img, gdImageGetPixel(img, i, sy-1));
    right = greyscale(img, gdImageGetPixel(img, i+1, sy-1));
    left = greyscale(img, gdImageGetPixel(img, i-1, sy-1));
    up = greyscale(img, gdImageGetPixel(img, i, sy-2));

    imageGradient[i * sy + (sy-1)] = (int)(1./sqrt(3.)*sqrt((right-left)*(right-left)/
    4.+(up-center)*(up-center))+0.5);
  }

  // left line

  for (j=1;j<sy-1;++j) {
    center  = greyscale(img, gdImageGetPixel(img, 0, j));
    right = greyscale(img, gdImageGetPixel(img, 1, j));
    up = greyscale(img, gdImageGetPixel(img, 0, j-1));
    down = greyscale(img, gdImageGetPixel(img, 0, j+1));

    imageGradient[0 * sy + j] = (int)(1./sqrt(3.)*sqrt((right-center)*(right-center)+
    (up-down)*(up-down))+0.5);
  }

  // right line

  for (j=1;j<sy-1;++j) {
    center  = greyscale(img, gdImageGetPixel(img, sx-1, j));
    left = greyscale(img, gdImageGetPixel(img, sx-2, j));
    up = greyscale(img, gdImageGetPixel(img, sx-1, j-1));
    down = greyscale(img, gdImageGetPixel(img, sx-1, j+1));

    imageGradient[(sx-1) * sy + j] = (int)(1./sqrt(3.)*sqrt((center-left)*
    (center-left)+(up-down)*(up-down))+0.5);
  }

  for (i=1;i<sx-1;++i)
    for (j=1;j<sy-1;++j) {
      center = greyscale(img, gdImageGetPixel(img, i, j));
      up = greyscale(img, gdImageGetPixel(img, i, j-1));
      down = greyscale(img, gdImageGetPixel(img, i, j+1));
      right = greyscale(img, gdImageGetPixel(img, i+1, j));
      left = greyscale(img, gdImageGetPixel(img, i-1, j));

      imageGradient[i * sy + j] = 
        (int)(1./sqrt(3.)*sqrt(
          (right-left)*(right-left)/4. +
          (up-down)*(up-down)/4. )+0.5);  
    }

  // return *imageGradient; 
}

int main(int argc, char *argv[]) {
  FILE *f, *out;
  gdImagePtr img, imgout; 
  // int c, sx, sy; 

  f = fopen(argv[1], "rb"); 
  img = gdImageCreateFromPng(f);
  fclose(f);

  int sx, sy; 
  sx = gdImageSX(img);
  sy = gdImageSY(img);

  // int c = gdImageGetPixel(img, 0, 0);
  // int imageGradient[sx][sy]; 
  int* imageGradient; 
  imageGradient = malloc(sx * sy * sizeof(*imageGradient)); 
  gradient(imageGradient, img, sx, sy); 
  
  // int* newImageGradient; 
  // newImageGradient = malloc(sx * sy * sizeof(*newImageGradient)); 
  
  // @@ this is where we need to repeat from

  int seams, i, j; 
  int maxseams = atoi(argv[3]); 
  for (seams=0; seams<maxseams; seams++) { 
    // find a path from left to right with minimal energy
    int weight[sx][sy];
    i = sx-1;
  
    // if we finish on that pixel we have to pay
    for (j=0;j<sy;++j)
      weight[sx-1][j] = imageGradient[i * sy + j];

    for (i=sx-2;i>=0;--i) {
      weight[i][0] = min(weight[i+1][0],weight[i+1][1]) + imageGradient[i * sy + 0];
      weight[i][sy-1] = min(weight[i+1][sy-1], weight[i+1][sy-2]) +                         
                          imageGradient[i * sy + (sy-1)];

      for (j=1;j<sy-1;++j)
        weight[i][j] = min(min(weight[i+1][j], weight[i+1][j-1]), 
          weight[i+1][j+1]) + imageGradient[i * sy + j];
    }
  
    int bestj[sx];
    bestj[0] = 0;

    // find best starting point
    for(j=0;j<sy;++j) { 
      if (weight[0][j] < weight[0][bestj[0]])
        bestj[0] = j;
    }

    srand((unsigned long)time(NULL));

    for (i=1;i<sx;++i) { 
      if (bestj[i-1] == 0)
        bestj[i] =  (weight[i][0] <= weight[i][1])?0:1;
      else if (bestj[i-1] == sy-1)
        bestj[i] =  (weight[i][sy-1] <= weight[i][sy-2])?(sy-1):(sy-2);
      else if (weight[i][bestj[i-1]] <= weight[i][bestj[i-1]-1] &&     
               weight[i][bestj[i-1]] <= weight[i][bestj[i-1]+1])
        bestj[i] = bestj[i-1];
      else if (weight[i][bestj[i-1]-1]==weight[i][bestj[i-1]+1])
        bestj[i] = (rand()/(RAND_MAX+1.)>=0.5)?(bestj[i-1]+1):(bestj[i-1]-1);
      else if (weight[i][bestj[i-1]-1] <= weight[i][bestj[i-1]] && 
               weight[i][bestj[i-1]-1] <= weight[i][bestj[i-1]+1])
        bestj[i] = bestj[i-1]-1;
      else bestj[i] = bestj[i-1]+1;
    }

    imgout = gdImageCreateTrueColor(sx, sy-1);

    int jj;
    // int newImageGradient[sx-1][sy-1]; 

    for(i=0;i<sx;++i) {
      jj = 0;
      for(j=0;j<sy;++j) {
        if(j != bestj[i]) {
          int c = gdImageGetPixel(img, i, j);
          gdImageSetPixel(imgout, i, jj, 
               gdTrueColor(gdImageRed(img, c),
                           gdImageGreen(img, c),
                           gdImageBlue(img, c)));
          // newImageGradient[i * sy + j] = imageGradient[i * sy + j]; 
          jj++;
        }
      }
    }
    
    // c = gdImageGetPixel(img, 0, 0);
    img = imgout; 
    // int imageGradient[sx-1][sy-1]; 
    // int *const tempImageGradient = imageGradient; 
    // imageGradient = newImageGradient; 
    // newImageGradient = tempImageGradient; 
    /* sx--; */ sy--; 
    gradient(imageGradient, img, sx, sy); 
    // (int[sx-1][sy-1])*newImageGradient; 
    // gradient(*imageGradient, img, sx, sy); 
  }
  
  free(imageGradient); 
  // free(newImageGradient); 

  out = fopen(argv[2], "wb");
  gdImagePng(imgout, out);
  gdImageDestroy(imgout);
  fclose(out);

  return 0;
}