% Entry point
function SeamCarving()
  Image = imread('tiger.bmp');  
  Image = double(Image);
  Result = Process(Image, 1, true);
  imwrite(uint8(Result), 'tiger_smaller_1.bmp'); 
end

% Preprocesses the image (rotates back and forth) before repeatedly calling
% the seam carving function
function Result = Process(InputImage, seamCount, isVertical)
  if isVertical == true % vertical resize
    for i = 1 : seamCount
      [InputImage, Debug] = DoSeamCarving(InputImage);
      imwrite(uint8(Debug), strcat(int2str(i), 'debug.bmp'));
    end
    Result = InputImage;
  else % horizontal resize
    % rotate it 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));
    for i = 1 : seamCount      
      [RotatedImage, Debug] = DoSeamCarving(RotatedImage);
      %imwrite(uint8(Debug), strcat(int2str(i), 'debug.bmp'));      
    end        
    % rotate the result back
    if i == seamCount
      Result(:, :, 1) = rot90(RotatedImage(:, :, 1), 3);
      Result(:, :, 2) = rot90(RotatedImage(:, :, 2), 3);
      Result(:, :, 3) = rot90(RotatedImage(:, :, 3), 3);      
    end    
  end
end

% Carves 1 seam
function [Result, Debug] = DoSeamCarving(InputImage)
  % compute the gradient map
  GrayImage = rgb2gray(uint8(InputImage));
  GrayImage = double(GrayImage);
  GradientImage = ComputeGradient(GrayImage);  
  % 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   
  % builds the result by tracking back
  Result(:, :, :) = zeros(height, width - 1, 3);
  Debug(:, :, :) = zeros(height, width, 3);  
  [minValue, minIndex] = min(Memo(height, :));
  for row = height : -1 : 1    
    if minIndex > 1      
      Result(row, 1 : minIndex - 1, :) = InputImage(row, 1 : minIndex - 1, :);          
    end
    if minIndex == 1
      Result(row, 1, :) = InputImage(row, 2, :);      
    end
    if minIndex < width      
      Result(row, minIndex : width - 1, :) = InputImage(row, minIndex + 1 : width, :);              
    end    
    % debug info
    
    if minIndex > 1      
      Debug(row, 1 : minIndex - 1, :) = InputImage(row, 1 : minIndex - 1, :);        
    end
    Debug(row, minIndex, :) = [255, 0, 0];    
    if minIndex < width      
      Debug(row, minIndex + 1 : width, :) = InputImage(row, minIndex + 1 : width, :);              
    end
    
    % advance to the previous row
    if row > 1
      leftCol = max(minIndex - 1, 1);
      rightCol = min(minIndex + 1, width);
      if Memo(row - 1, leftCol) < Memo(row - 1, minIndex) && ...
         Memo(row - 1, leftCol) < Memo(row - 1, rightCol)
         minIndex = leftCol;
      end
      if Memo(row - 1, rightCol) < Memo(row - 1, minIndex) && ...
         Memo(row - 1, rightCol) < Memo(row - 1, leftCol)
         minIndex = rightCol;
      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