/* seamless.c - Seam Carving to Resize Images
   Source: http://inamidst.com/code/seamless.c
   
   Written by Arthur B., http://slashdot.org/~Arthur+B./
   Cf. http://slashdot.org/comments.pl?sid=279819&cid=20359241
   Modified by Sean B. Palmer, http://inamidst.com/sbp/

   Contributions of help from:
    * Adam Wendt, http://blog.thelsdj.com/
    * Björn Höhrmann, http://bjoern.hoehrmann.de/
    * Kevin Reid, http://homepage.mac.com/kpreid/

   Requirements: http://www.libgd.org/
   Usage: ./seamless <input> <output> <seams> [ "v" ]
     e.g. ./seamless in.png out.jpg 150 v
          -> Take 150 vertical seams from in.png

   gcc -O3 -o seamless -lgd -L/opt/local/lib -I/opt/local/include seamless.c
   Caveat: segfaults on very large images; dunno why
*/

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

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

int greyscale(int *grey, gdImagePtr img, int sx, int sy) { 
  int i, j; 
  for (i=1; i<sx-1; ++i) {
    for (j=1; j<sy-1; ++j) { 
      int c = gdImageGetPixel(img, i, j); 
      int r = gdImageRed(img, c);
      int g = gdImageGreen(img, c);
      int b = gdImageBlue(img, c);
      grey[i * sy + j] = (int)((r+g+b)/3.+0.5); 
    }
  }
}

void gradient(int *imageGradient, int *grey, int sx, int sy) { 
  int up, down, right, left, center; 
  int i, j; 
  
  // 0,0
  right = grey[1 * sy + 0];
  down = grey[0 * sy + 1];
  center = grey[0 * sy + 0];
  
  imageGradient[0 * sy + 0] = (int)(1./sqrt(3.)*sqrt((right-center)*
     (right-center)+(center-down)*(center-down))+0.5);
  
  // sx-1,0
  left = grey[(sx-2) * sy + 0];
  down = grey[(sx-1) * sy + 1];
  center = grey[(sx-1) * sy + 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 = grey[1 * sy + (sy-1)];
  up = grey[0 * sy + (sy-2)];
  center = grey[0 * sy + (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 = grey[(sx-2) * sy + (sy-1)];
  up = grey[(sx-1) * sy + (sy-2)];
  center = grey[(sx-1) * sy + (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 = grey[i * sy + 0];
    right = grey[(i+1) * sy + 0];
    left = grey[(i-1) * sy + 0];
    down = grey[i * sy + 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 = grey[i * sy + (sy-1)];
    right = grey[(i+1) * sy + (sy-1)];
    left = grey[(i-1) * sy + (sy-1)];
    up = grey[i * sy + (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 = grey[0 * sy + j];
    right = grey[1 * sy + j];
    up = grey[0 * sy + (j-1)];
    down = grey[0 * sy + (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 = grey[(sx-1) * sy + j];
    left = grey[(sx-2) * sy + j];
    up = grey[(sx-1) * sy + (j-1)];
    down = grey[(sx-1) * sy + (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 = grey[i * sy + j];
      up = grey[i * sy + (j-1)];
      down = grey[i * sy + (j+1)];
      right = grey[(i+1) * sy + j];
      left = grey[(i-1) * sy + j];

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

void doc() { 
  printf("Usage: ./seamless <input> <output> <seams> [ v ]\n"); 
  printf("  e.g. ./seamless in.png out.jpg 150 v\n"); 
  printf("    -> Take 150 vertical seams from in.png\n"); 
  exit(1); 
}

int main(int argc, char *argv[]) { 
  int vertical; 
  if (argc == 4) vertical = 0; 
  else if (argc == 5) vertical = 1; 
  else doc(); 

  char *fn = argv[1]; 
  gdImagePtr eimg; 

  FILE *f = fopen(fn, "rb"); 
  if (strcasecmp(strrchr(fn, '.'), ".jpg") == 0) 
    eimg = gdImageCreateFromJpeg(f); 
  else eimg = gdImageCreateFromPng(f);
  fclose(f);

  int sx = gdImageSX(eimg); 
  int sy = gdImageSY(eimg); 

  gdImagePtr img; 
  if (vertical == 0) {
    img = gdImageCreateTrueColor(sy, sx); 
    gdImageCopyRotated(img, eimg, ((float)sy)/2,((float)sx)/2,0,0,sx,sy,90); 
    gdImageDestroy(eimg);
    sx = gdImageSX(img); 
    sy = gdImageSY(img); 
  } else img = eimg; 

  int* imageGradient = malloc(sx * sy * sizeof(*imageGradient)); 
  int* grey = malloc(sx * sy * sizeof(*grey)); 
  greyscale(grey, img, sx, sy); 
  gradient(imageGradient, grey, sx, sy); 

  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;
    }

    gdImagePtr imgout = gdImageCreateTrueColor(sx, sy-1);

    int jj;
    for(i=0; i<sx; ++i) {
      jj = 0; 
      for(j=0; j<sy; ++j) {
        if(j != bestj[i]) {
          int p = gdImageGetTrueColorPixel(img, i, j); 
          gdImageSetPixel(imgout, i, jj, p); 
          jj++;
        }
      }
    }

    sy--; 
    img = imgout; 
    greyscale(grey, img, sx, sy); 
    gradient(imageGradient, grey, sx, sy); 
  }
  
  free(imageGradient); 

  gdImagePtr oimg; 
  if (vertical == 0) { 
    oimg = gdImageCreateTrueColor(sy, sx); 
    gdImageCopyRotated(oimg, img, ((float)sy)/2,((float)sx)/2,0,0,sx,sy,270); 
    gdImageDestroy(img);
  } else oimg = img; 

  char *ofn = argv[2]; 
  FILE *out = fopen(ofn, "wb");
  if (strcasecmp(strrchr(ofn, '.'), ".jpg") == 0) 
    gdImageJpeg(oimg, out, 80); 
  else gdImagePng(oimg, out);
  gdImageDestroy(oimg);
  fclose(out);

  return 0;
}