%COLORSAMP   CCD color subsampling and linear interpolation reconstruction.
%   Given an RGB image, performs color subsampling by selecting a single
%   color channel for each column in alternating order (RGBRGB...).
%   Constructs and plots a grayscale image showing these sampled values.
%   Constructs and plots an approximation to the original RGB image using
%   linear interpolation.
%
%   For problem 2, run this file twice as follows:
%     colorsamp(imread('brettan.tif'));
%     colorsamp(imread('brettanBlurred.tif'));
%   You should see that blurring noticeably reduces the color fringes
%   obviously present in the unblurred reconstruction, at the price of
%   a loss in high-resolution detail.
%
%     [sampImage, interpImage] = colorsamp(origImage)
%    
% Parameters:
%   origImage = RGB color image
% Outputs:
%   sampImage = grayscale image showing intensity of color channel 
%               extracted for each pixel
%   interpImage = linearly interpolated fullcolor reconstruction of origImage,
%                 based on the samples in sampImage

% 6.801 Problem Set 0: Problem 2
% Erik Sudderth

function [sampImage,interpImage] = colorsamp(origImage)

  % Determine input image size, and convert to double format
  [M,N,d] = size(origImage);
  if (d ~= 3)
    error('Input image must be RGB.');
  end
  origImage = im2double(origImage);
  
  % Plot original image
  figure; imshow(origImage);
  title('Original');

  % Extract color subsamples
  rSamp = origImage(:,1:3:end,1);
  gSamp = origImage(:,2:3:end,2);
  bSamp = origImage(:,3:3:end,3);
  
  % Create single grayscale image with color subsampled by columns
  sampImage = zeros(M,N);
  sampImage(:,1:3:end) = rSamp;
  sampImage(:,2:3:end) = gSamp;
  sampImage(:,3:3:end) = bSamp;
 
  figure; imshow(sampImage);
  title('Color Subsampling');
  
  % Transform back to full color by linear interpolation
  interpImage = zeros(M,N,3);
  
  lastRed = 3*floor((N-1)/3) + 1;
  interpImage(:,1:3:lastRed,1) = rSamp;
  interpImage(:,2:3:lastRed,1) = (2/3)*rSamp(:,1:end-1) + (1/3)*rSamp(:,2:end);
  interpImage(:,3:3:lastRed,1) = (1/3)*rSamp(:,1:end-1) + (2/3)*rSamp(:,2:end);
  for (i = lastRed+1:N)
    interpImage(:,i,1) = rSamp(:,end);
  end
  
  lastGreen = 3*floor((N-2)/3) + 2;
  interpImage(:,1,2) = gSamp(:,1);
  interpImage(:,2:3:lastGreen,2) = gSamp;
  interpImage(:,3:3:lastGreen,2) = (2/3)*gSamp(:,1:end-1) + (1/3)*gSamp(:,2:end);
  interpImage(:,4:3:lastGreen,2) = (1/3)*gSamp(:,1:end-1) + (2/3)*gSamp(:,2:end);
  for (i = lastGreen+1:N)
    interpImage(:,i,2) = gSamp(:,end);
  end
  
  lastBlue = 3*floor((N-3)/3) + 3;
  interpImage(:,1:2,3) = [bSamp(:,1) bSamp(:,1)];
  interpImage(:,3:3:lastBlue,3) = bSamp;
  interpImage(:,4:3:lastBlue,3) = (2/3)*bSamp(:,1:end-1) + (1/3)*bSamp(:,2:end);
  interpImage(:,5:3:lastBlue,3) = (1/3)*bSamp(:,1:end-1) + (2/3)*bSamp(:,2:end);
  for (i = lastBlue+1:N)
    interpImage(:,i,3) = bSamp(:,end);
  end
  
  figure; imshow(interpImage);
  title('Linear Interpolation');