Copyright 2010 Mitsubishi Electric Research Laboratories All Rights Reserved.

Permission to use, copy and modify this software and its documentation without fee for educational, research and non-profit purposes, is hereby granted, provided that the above copyright notice and the following three paragraphs appear in all copies.

To request permission to incorporate this software into commercial products contact: Vice President of Marketing and Business Development; Mitsubishi Electric Research Laboratories (MERL), 201 Broadway, Cambridge, MA 02139 or .

IN NO EVENT SHALL MERL BE LIABLE TO ANY PARTY FOR DIRECT, INDIRECT, SPECIAL, INCIDENTAL, OR CONSEQUENTIAL DAMAGES, INCLUDING LOST PROFITS, ARISING OUT OF THE USE OF THIS SOFTWARE AND ITS DOCUMENTATION, EVEN IF MERL HAS BEEN ADVISED OF THE POSSIBILITY OF SUCH DAMAGES.

MERL SPECIFICALLY DISCLAIMS ANY WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. THE SOFTWARE PROVIDED HEREUNDER IS ON AN "AS IS" BASIS, AND MERL HAS NO OBLIGATIONS TO PROVIDE MAINTENANCE, SUPPORT, UPDATES, ENHANCEMENTS OR MODIFICATIONS.



%==================================================================
% Same code is available as part of Complete Matlab Code and Datasets
% Motion Deblurring Source Code for a Simple Case Shown in Figure 3
% Solves for 1D motion blur assuming object is moving from top to bottom
% Using the 52 element binary sequence for fluttering
 
%==============================================================
m = 52; % Coded Sequence length
CodeSeq=double('1010000111000001010000110011110111010111001001100111')-'0';
% Read input image
InputImage = double(readpfm_color('InputTaxi.pfm')); 
[H,W,CH] = size(InputImage);
k = [235]; % Assume Blur Size in pixels = 235
 
% Resize image height by m/k so that the effective blur is m
InputImage1 = imresize(InputImage,[ceil(H*m/k) W],'bicubic');
 
% Get object size, n, knowing blurredSize = n + m - 1
n = size(InputImage1,1) - m + 1;
 
% Foreground: Get A matrix for foreground which encodes the blurring
Af = ComposeMotionBlurMatrix(CodeSeq, n);
 
% Background: bk is contribution vector of bkgrnd in blurred img for all pix
bk = abs(1 - Af*ones(n,1));
 
% Assume constant background in first m and last m pixels (effective blur is m)
bkLeft = zeros(size(bk)); bkLeft(1:m)=bk(1:m);
bkRite = zeros(size(bk)); bkRite(end-m+1:end)=bk(end-m+1:end);


% Ready to Solve AX=B for each color channel
A = [Af  bkLeft  bkRite];

for colorchannel = 1:CH
    B = InputImage1(:,:,colorchannel); % coded img for channel

        X = A\B; %Least square solution
    
    OutColorImage(:,:,colorchannel) = X(1:end-2,:);
end

% Expand/shrink the deblurred image to match blur size of k 
OutColorImage = imresize(OutColorImage,[H-k+1 W],'bicubic');



%==================================================================

function [A] = ComposeMotionBlurMatrix(CodeSeq,n)
CodeSeq = CodeSeq(:)/sum(CodeSeq);
k = size(CodeSeq,1);
ZeroPaddedCodeSeq = sparse([CodeSeq; zeros(n-1,1)]);
A = gallery('circul',ZeroPaddedCodeSeq)';% Create Circulant Matrix
A = A(:,1:n); % Keep first n columns out of (n+k-1) columns
 
%=================================================================