% Entry point
function SeamCarving()
  Image = imread('beach.bmp');
  Image = double(Image);
  Result = Process(Image, 128, false);
  imwrite(uint8(Result), 'beach_larger_vertical_32.bmp');  
end

% Pre and post processes the input before enlarging
function Result = Process(InputImage, seamCount, isVertical)
  if isVertical == true % horizontal resize
    [InputImage, Debug] = DoSeamCarving(InputImage, seamCount);
    %imwrite(uint8(Debug), strcat(int2str(i), 'debug.bmp'));    
    Result = InputImage;
  else % vertical resize
    % rotate the input 90 degree
    RotatedImage = zeros(size(InputImage, 2), size(InputImage, 1), 3);
    RotatedImage(:, :, 1) = rot90(InputImage(:, :, 1));
    RotatedImage(:, :, 2) = rot90(InputImage(:, :, 2));
    RotatedImage(:, :, 3) = rot90(InputImage(:, :, 3));    
    [RotatedImage, Debug] = DoSeamCarving(RotatedImage, seamCount);
    %imwrite(uint8(Debug), strcat(int2str(i), 'debug.bmp'));                
    % rotate the output back
    Result(:, :, 1) = rot90(RotatedImage(:, :, 1), 3);
    Result(:, :, 2) = rot90(RotatedImage(:, :, 2), 3);
    Result(:, :, 3) = rot90(RotatedImage(:, :, 3), 3);    
  end
end

% Enlarges the image by seamCount seams
function [Result, Debug] = DoSeamCarving(InputImage, seamCount)
  % compute the gradient map
  GrayImage = rgb2gray(uint8(InputImage));
  GrayImage = double(GrayImage);
  GradientImage = ComputeGradient(GrayImage);
  %AverageImage = imfilter(InputImage, [0.5 0.0 0.5], 'replicate');    
  % dynamic programming to compute seams
  height = size(GrayImage, 1);
  width = size(GrayImage, 2);
  Memo = zeros(height, width);
  Memo(1, :) = GradientImage(1, :);
  for row = 2 : height
    for col = 1 : width
      leftCol = max(col - 1, 1);
      rightCol = min(col + 1, width);      
      Memo(row, col) = GradientImage(row, col) + ...
        min([Memo(row - 1, leftCol), Memo(row - 1, col), Memo(row - 1, rightCol)]);
    end
  end
  
  % build the output by tracking back
  Result(:, :, :) = zeros(height, width + seamCount, 3);
  Debug(:, :, :) = zeros(height, width, 3);  
  Temp = Memo(height, :);
  [Temp, Ix] = sort(Temp);
  minIndex = Ix(1 : 2 : seamCount * 2);  
  for row = height : -1 : 1     
    i = 1;
    j = 1;
    k = 1;    
    minIndex = sort(minIndex);
    while i <= width + seamCount
      if k <= seamCount && j == minIndex(k)        
        t = 2;
        while k < seamCount && minIndex(k) == minIndex(k + 1)
          t = t + 1;
          k = k + 1;
        end                                       
        Result(row, i, :) = InputImage(row, j, :);        
        for u = 2 : t
          % duplicate to the right                
          Result(row, i + u - 1, :) = InputImage(row, j, :);
          % advance to the previous row
          if row > 1
            leftCol = max(minIndex(k) - 1, 1);
            rightCol = min(minIndex(k) + 1, width);
            col = minIndex(k);
            if Memo(row - 1, leftCol) < Memo(row - 1, minIndex(k)) && ...
               Memo(row - 1, leftCol) < Memo(row - 1, rightCol)
               minIndex(k) = leftCol;
            end
            if Memo(row - 1, rightCol) < Memo(row - 1, minIndex(k)) && ...
               Memo(row - 1, rightCol) < Memo(row - 1, leftCol)
               minIndex(k) = rightCol;
            end
            
            
            if k > 1
              if minIndex(k) == minIndex(k - 1)
                if minIndex(k) == leftCol
                  if Memo(row - 1, col) <= Memo(row - 1, rightCol)
                    minIndex(k) = col;
                  else
                    minIndex(k) = rightCol;
                  end
                elseif minIndex(k) == col
                  minIndex(k) = rightCol;
                end                
                if minIndex(k - 1) > 1
                  minIndex(k - 1) = minIndex(k - 1) - 1;
                end
              end
            end
                   
          end                            
        end                
        i = i + t;
        j = j + 1;
        k = k + 1;
      else
        % just copy
        Result(row, i, :) = InputImage(row, j, :);        
        i = i + 1;
        j = j + 1;
      end          
    end
  end
end

% Computes the gradient map
function GradientImage = ComputeGradient(InputImage)
  SobelFilter = fspecial('sobel');
  HorizontalGradientImage = imfilter(InputImage, SobelFilter);
  VerticalGradientImage = imfilter(InputImage, SobelFilter');
  GradientImage = abs(HorizontalGradientImage) + abs(VerticalGradientImage);    
end