function slide=optdeblur
% This demo shows how the Optimization Toolbox can be used to 
% recover a blurred image.  It does this by solving a 
% large-scale constrained linear least-squares optimization
% problem.

% Javier Pena, 8-18-95, Matthew J. Simoneau, 09-04-98
% Copyright 1990-2000 The MathWorks, Inc.
% $Revision: 1.8 $ $Date: 2000/06/16 22:25:15 $

if nargout<1,
  playshow optdeblur
else
  %========== Slide 1 ==========

  slide(1).code={
   'if ~isappdata(gcf,''Initialized'')',
   '  set(gcf,''defaultaxesposition'',get(gca,''position''));',
   '  setappdata(gcf,''Initialized'',1);',
   'end;',
   '',
   'load optdeblur;',
   'subplot(1,1,1);',
   'image(colorcar);',
   'axis image;',
   'axis off;',
   '' };
  slide(1).text={
     ['This demo shows how the Optimization Toolbox can be used to solve a '...
           'large-scale constrained linear least-squares optimization '...
           'problem to recover a blurred image.  We will add motion '...
           'blur to this photo of Mary Ann and Matthew sitting in Joe''s '...
           'car, then try to restore the original.'],
   '>> load optdeblur;',
   '>> image(colorcar);',
   '>> axis image;',
   '>> axis off;',
   ''};

  %========== Slide 2 ==========

  slide(2).code={
   'imagesc(P);',
   '[m,n]=size(P);',
   'mn=m*n;',
   'colormap(gray);',
   'axis image;',
   'axis off;',
   'title([num2str(m) '' x '' num2str(n) '' ('' num2str(mn) '' pixels)'']);' };
slide(2).text={
   ['Our starting image is the black and white image shown above, '...
         'contained the m x n matrix P.  Each element in the matrix '...
         'represents a pixel''s gray intensity between black and white '...
         '(0 and 1).'],
   '>> [m,n]=size(P);',
   '>> mn=m*n;',
   '>> imagesc(P);',
   '>> colormap(gray);',
   '>> axis image;  axis off;',
   '>> title([num2str(m) '' x '' num2str(n) '' ('' num2str(mn) '' pixels)'']);',
   ''};

  %========== Slide 3 ==========

slide(3).code={
   'cla;',
   'axis ij',
   'axis off;',
   'xs=25;',
   'ys=15;',
   'xlim([0,xs+1]);',
   'ylim([0,ys+1]);',
   '[ix,iy]=meshgrid(1:(xs-1),1:(ys-1));',
   'l=abs(ix-iy)<=5;',
   'text(ix(l),iy(l),''x'')',
   'text(ix(~l),iy(~l),''0'')',
   'text(xs*ones(ys,1),1:ys,''...'');',
   'text(1:xs,ys*ones(xs,1),''...'');',
   'title(''Blurring Operator D ( x = 1/11)'')',
   '' };
slide(3).text={
   ['We can simulate the effect of vertical motion blurring by averaging '...
         'each pixel with the 5 above and below.  We construct a sparse '...
         'matrix D, shown in part above, that will do this with a single '...
         'matrix multiply.'],
   '>> blur=5;  mindex=1:mn;  nindex=1:mn;',
   '>> for i=1:blur,',
   '>>  mindex=[mindex i+1:mn 1:mn-i]; nindex=[nindex 1:mn-i i+1:mn];',
   '>> end',
   '>> D=sparse(mindex,nindex,1/(2*blur+1));'};

  %========== Slide 4 ==========

  slide(4).code={
   'subplot(1,1,1);',
   'imagesc(blurred);',
   'axis image;',
   'axis off;',
   '' };
  slide(4).text={
     ['We multiply the image by this operator to create the blurred image.  '...
        'P is the original image, D is the operator, and G is the blurred '...
        'image shown above.'],
   '>> G=D*P(:);',
   '>> imagesc(reshape(G,m,n));',
   '>> axis image;',
   '>> axis off;'};

  %========== Slide 5 ==========

  slide(5).code={
   'subplot(1,2,1);',
   'imagesc(blurred); axis image; axis off;',
   'title(''G, the blurred image'');',
   'subplot(1,2,2);',
   'imagesc(reshape(P(:),m,n)); axis image; axis off;',
   'title(''P, the original image'');' };
  slide(5).text={
     ['Now, let''s pretend Joe took this blurred picture G from a moving '...
           'elevator.  Assume we know how fast the elevator is moving, so we '...
           'know the blurring operator D.  How well can we remove the blur '...
           'and recover the original image P?'],
     ''
     'The simplest approach is to solve the least squares problem:',
     '   min( || D*P(:) - G||^2 )'};

  %========== Slide 6 ==========

slide(6).code={
   'subplot(1,1,1) ;',
   'delete(gca);',
   'axis ij',
   'axis off;',
   'xs=11;',
   'ys=11;',
   'xlim([0,xs+1]);',
   'ylim([0,ys+1]);',
   '[ix,iy]=meshgrid(1:(xs-1),1:(ys-1));'
   'four=(ix==iy);',
   'one=(abs(ix-iy)==1);',
   'text(ix(one),iy(one),''-1'')',
   'text(ix(four),iy(four),''+4'')',
   'text(ix(~four & ~one),iy(~four & ~one),'' 0'')',
   'text(xs*ones(ys,1),1:ys,''...'');',
   'text(1:xs,ys*ones(xs,1),''...'');',
   'title(''L, The Discrete Laplacian'')',
   '' };
slide(6).text={
   ['In practice the results obtained with this simple approach tend '...
         'to be noisy. To compensate for this, a regularization term is added:'],
   '      0.0004*|| L*P ||^2',
   ['L is the discrete Laplacian, shown '...
         'above, which relates each pixel to those surrounding it.  Since '...
         'we know we are looking for a gray intensity, we also impose the '...
         'constraint that the elements of P must fall between 0 and 1.'],
   '>>L = sparse( [1:mn,2:mn,1:mn-1],  [1:mn,1:mn-1,2:mn], ...',
   '       [4*ones(1,mn) -1*ones(1,2*(mn-1))]  );'};

  %========== Slide 7 ==========

  slide(7).code={
   'delete(gca);',
   'axis on;',
   'line([0 mn],[0 mn],''linewidth'',3);',
   'line([mn 2*mn],[0 mn],''linewidth'',3);',
   'xlabel([''nz = 648714'']);',
   ['set(gca,''xlim'',[0 2*(mn+1)],''ylim'',[0 mn+1],''ydir'','...
   	'''reverse'',''grid'',''none'');'],
   'title(''Structure of Matrix A'''''');',
   '' };
  slide(7).text={
   'To obtain the deblurred picture we want to solve for P:',
   '     min( || D*P(:) - G(:) ||^2 + 0.0004*|| L*P(:) ||^2 )',
   'We can simplify this expression by defining A and b:',
   '>> A = [D ; 0.02*L];  b = [ G(:) ; zeros(mn,1) ];',
   'which changes the above equation to:',
   '     min( || A*P(:) - b ||^2 )',
   ['subject to 0<=P<=1.  Both matrices D and L relate each pixel to a few '...
         'of its neighbors.  This makes A structured and sparse.'],
   '>> spy(A'');   (simulated above)'};

  %========== Slide 8 ==========

slide(8).code={
   'subplot(1,1,1);',
   'delete(gca)',  
   'imagesc(deblurred);',
   'axis image;',
   'axis off;',
   '' };
slide(8).text={
   ['Because A is sparse, we can use a large-scale algorithm to solve '...
         'this linear least squares optimization problem.  We call '...
         'LSQLIN with A, b, lower bounds, upper bounds, and options:'],
   '>> options = optimset(''LargeScale'', ''on'',''Display'', ''off'' )',
   '>> x = lsqlin(A, b, [], [], [], [], zeros(mn,1), ones(mn,1), [], options);'
   ['Due to the size of the optimization problem this process takes '...
         'about 10 minutes.  For this demo, the solution has been '...
         'previously calculated and saved.'],
   '>> imagesc(reshape(x,m,n)); axis image; axis off;'};

  %========== Slide 9 ==========

  slide(9).code={
   'subplot(1,2,1);',
   'imagesc(blurred); axis image; axis off;',
   'title(''Blurred'');',
   '',
   'subplot(1,2,2);',
   'imagesc(deblurred); axis image; axis off;',
   'title(''De-Blurred'');' };
  slide(9).text={
     ['Let''s compare the blurred and deblurred pictures.  You can '...
        'enlarge this window for a better view.'],
   '>> subplot(1,2,1);',
   '>> imagesc(reshape(G,m,n)); axis image; axis off;',
   '>> title(''Blurred'');',
   '',
   '>> subplot(1,2,2);',
   '>> imagesc(reshape(x,m,n)); axis image; axis off;',
   '>> title(''De-Blurred'');',
   ''};
end