decom2d.m

Sparse Signal Representation using Overlapping Frames (matlab toolbox)

function Y = Decom2D(TransMet, X, arg3, arg4)
% Decom2D    Decompose a 2D-signal (image) into expansion coefficients, 
%            or the inverse function if TransMet is negative.
% Note: if X is uint8 and Y is uint16, then Y is multiplied by 2 (to keep
% it more accurate, and rounded) and thus energy of Y is approx. 4 times energy
% in X, if exact values are needed use double type.
% Note: Decom2D is "more similar to" Ttimes than to Decom1D
%
% Examples:
% Y = Decom2D(TransMet, X);      % decomposition (analysis) by TransMet
% X = Decom2D(-TransMet, Y);     % inverse decomposition (reconstruction) by TransMet
% ---------------------------------------------------------------------------
% arguments:
%  X        the signal, usually an MxN image, may be uint8 or double
%  Y        the decomposition coefficients, usually also of size MxN
%           may be uint16 or double
%  TransMet An integer for the decomposition method
%           use negative integer for the inverse
%            0    : identity transform, Y=X
%            1    : 2x2 Haar wavelet, LP component is split recursively
%            2    : a dyadic wavelet (given by arg3, name as in wfilters.m) 
%                   number of levels is given by arg4,
%                   arg3 and arg4 are used in GetWave.m
%                   works only for orthogonal wavelets and the db79 filters!
%            3- 4 : may be used for other special transforms
%            5- 8 : NxN DCT transform, N=[2,4,8,16]
%            9-12 : 2N*N LOT,          N=[2,4,8,16]
%           13-16 : 4N*N ELT,          N=[2,4,8,16]
% ---------------------------------------------------------------------------
% some problems exist for biorthogonal wavelets, TransMet==2
% the function works only for orthogonal wavelets and the db79 filters!
% The problems may be in GetWave

% ---------------------------------------------------------------------------
% 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:  dd.mm.yyyy
% Ver. 1.0  03.08.2000  Karl Skretting, function made 
% Ver. 1.1  02.12.2002  KS: moved from ..\Frames2D to ..\FrameTools
% ---------------------------------------------------------------------------

Mfile='Decom2D';
Display=1;

% check input and output arguments, and assign values to arguments
if (nargin < 2); 
   error([Mfile,': function must have two input arguments, see help.']); 
end
if (nargout ~= 1); 
   error([Mfile,': function must have one output arguments, see help.']); 
end
TransMet=floor(TransMet);

[M,N]=size(X);
if isa(X,'uint8');
   Y=double(X)-2^7;
elseif isa(X,'uint16');
   Y=double(X)-2^15;
elseif isa(X,'double');
   Y=X;
else
   error([Mfile,': unknown type of input argument.']); 
end

if abs(TransMet)>4
    Nt=2^(mod(abs(TransMet)-1,4)+1);
    if (rem(M,Nt) | rem(N,Nt))
        % check that M and N have Nt as factor
        disp([Mfile,': size is not blocks of ',int2str(Nt),'x',...
                int2str(Nt),', decomposition using method ',int2str(TransMet),...
                'is not possible.']);
        return   % no decomposition is done
    end
    %
    if (N*M)>(2^18)         % resonable sized block
        Nn=2^ceil(log2(N*M)-18);
        if (mod(N,Nn) | mod(M,Nn)); Nn=Nt; end;
        for n=1:Nn
            Y(:,((n-1)*N/Nn+1):(n*N/Nn))=Ttimes(TransMet,Y(:,((n-1)*N/Nn+1):(n*N/Nn)));
        end
        for n=1:Nn
            Y(((n-1)*M/Nn+1):(n*M/Nn),:)=Ttimes(TransMet,Y(((n-1)*M/Nn+1):(n*M/Nn),:)')';
        end
    else
        Y=Ttimes(TransMet,Y);
        Y=Ttimes(TransMet,Y')';
    end
elseif abs(TransMet)==2
    % dyadic wavelet filter bank, implementad as a "transform"
    if (nargin < 2); 
        error([Mfile,': for TransMet==2, function must have four ',...
                'input arguments, see help.']); 
    end
    wname=[arg3,'       '];
    if (TransMet<0)
        if strcmp(wname(1:4),'db79')
            T=GetWave('db79a',arg4);
        else
            T=GetWave(arg3,arg4);
        end
        [NtP,Nt]=size(T);
        P=NtP/Nt;
        T=T';
        T=reshape(T,Nt,Nt,P);    % the analysis filter bank
        for p=1:P; T(:,:,p)=T(:,:,p)'; end;  % the synthesis filter bank
    else
        if strcmp(wname(1:4),'db79')
            T=GetWave('db79s',arg4);
        else
            T=GetWave(arg3,arg4);
        end
        [NtP,Nt]=size(T);
        P=NtP/Nt;
        T=T';
        T=reshape(T,Nt,Nt,P);    % the analysis filter bank
    end
    if (rem(M,Nt) | rem(N,Nt))
        % check that M and N have Nt as factor
        disp([Mfile,': size is not blocks of ',int2str(Nt),'x',...
                int2str(Nt),', decomposition using method ',int2str(TransMet),...
                'is not possible.']);
        return   % no decomposition is done
    end
    %
    if (N*M)>(2^18)         % resonable sized block
        Nn=2^ceil(log2(N*M)-18);
        if (mod(N,Nn) | mod(M,Nn)); Nn=Nt; end;
        for n=1:Nn
            Y(:,((n-1)*N/Nn+1):(n*N/Nn))=Ttimes(T,Y(:,((n-1)*N/Nn+1):(n*N/Nn)),(TransMet<0));
        end
        for n=1:Nn
            Y(((n-1)*M/Nn+1):(n*M/Nn),:)=Ttimes(T,Y(((n-1)*M/Nn+1):(n*M/Nn),:)',(TransMet<0))';
        end
    else
        Y=Ttimes(T,Y,(TransMet<0));
        Y=Ttimes(T,Y',(TransMet<0))';
    end
elseif abs(TransMet)==1
    Y=Ttimes(TransMet,Y,sqrt(1/2));
    Y=Ttimes(TransMet,Y',sqrt(1/2))';
end

if 0
    disp([min(Y(:)),max(Y(:))]);
end
if isa(X,'uint8');                      % used for image
    Y=uint16(round(Y*2+2^15));
elseif isa(X,'uint16');                 % used for expansion coefficients
    Y=uint8(round(Y/2+2^7));
end

return;