imdc.m

Sparse Signal Representation using Overlapping Frames (matlab toolbox)

function varargout = ImDC(I,Idc)
% ImDC       Split an image into DC component and residual, or opposite. 
%            We use a separable low-pass filter whith length 48 and downsample
% by 16 (in each dimension). (that is the synthesis filter is of length 48
% smooth without ringing, but the analysis filter is 'infinite') 
% Thus it is a biorthogonal filter.
% This is an ad-hoc filter, and is not optimal in any way.
%
% [I,Idc]=ImDC(I);             % split image
% I=ImDC(I,Idc);               % put images together again
% ---------------------------------------------------------------------------
% arguments:
%  I        the image, size MxN wher M and N have 16 as factors
%           image should be type double
%  Idc      the DC component image, size (M/16)x(N/16)
% ---------------------------------------------------------------------------

% ---------------------------------------------------------------------------
% Copyright (c) 1999.  Karl Skretting.  All rights reserved.
% Hogskolen in Stavanger (Stavanger University), Signal Processing Group
% Mail:  [e-mail]   Homepage:  http://www.ux.his.no/~karlsk/
% 
% HISTORY:
% Ver. 1.0  20.09.2000 Karl Skretting, function made 
% Ver. 1.1  27.11.2002  KS: moved to ..\FrameTools
% ---------------------------------------------------------------------------

Mfile='ImDC';

if nargin==1
   if nargout~=2
      disp([Mfile,': function called with one input argument.']); 
      error('then it must have two output arguments, see help.'); 
   end
elseif nargin==2
   if nargout~=1
      disp([Mfile,': function called with two input arguments.']); 
      error('then it must have one output argument, see help.'); 
   end
else   
   error([Mfile,': wrong number of input arguments, see help.']); 
end

[M,N]=size(I);
% find the synthesis vectors for low-pass
D=16;
f=blackman(3*D+2);f=f(2:(3*D+1));      % this makes an acceptable filter
f=f*(D/sum(f));
f1=f(1:D)*((1-f(D+1))/(f(D)+f(1)));
f3=flipud(f1);
f2=ones(D,1)-f1-f3;
f=[f1;f2;f3];
% find some matrices
Fm=zeros(M,M/D);
Fm(1:(2*D),1)=[f2+f3;f3];
for k=2:(M/D-1); Fm(((k-2)*D+1):((k+1)*D),k)=f; end;
Fm((M-2*D+1):M,M/D)=[f1;f1+f2];
Fn=zeros(N,N/D);
Fn(1:(2*D),1)=[f2+f3;f3];
for k=2:(N/D-1); Fn(((k-2)*D+1):((k+1)*D),k)=f; end;
Fn((N-2*D+1):N,N/D)=[f1;f1+f2];

if nargin==1
   Fmpi=inv(Fm'*Fm)*Fm';   % the pseudoinverse matrix (analysis)
   Fnpi=inv(Fn'*Fn)*Fn';   % the pseudoinverse matrix (analysis)
   Idc=Fmpi*I*Fnpi';
   Idc=floor(Idc/2+0.5);
   I=I-2*Fm*Idc*Fn';     % the difference image
   varargout(1)={I};
   varargout(2)={Idc};
elseif nargin==2
   I=I+2*Fm*Idc*Fn';     % add the low-pass image
   varargout(1)={I};
end

return


% example of use
clear all;
% Image=GetIm(10,5);              
[Image,ff]=GetIm(1);  % my own function that loads an image            
[Mi,Ni]=size(Image);
%make float of image 
if isa(Image,'uint8');
   Image=double(Image)-128;    % -128 <= x <= 127   (x is pixel value)
end

[Id,Idc]=ImDC(Image);       % the difference image and DC coefficients
Ir=ImDC(Id,Idc);            % the restored image
disp(['Difference is ',num2str(sum(abs(Image(:)-Ir(:))))]);
% the PSNR of the low bit rate image
PSNRnow=-10*log10(Id(:)'*Id(:)/(255*255*Mi*Ni))