技術標籤：學習記錄matlab

contourlet工具箱的使用之影象分解和重構的簡單介紹

DECDEMO.m檔案：是一個演示檔案，簡單的介紹瞭如何使用contourlet工具箱進行分解、重建一個影象，它提供了一個使用基本功能，如pdfbdec、pdfbrec、showpdfb。
它可以在影象分析等應用程式中進行修改，影象檢索和影象處理。顯示影象時，程式將暫停並等待您的響應。準備好後，只需按Enter鍵繼續。
%decdemo（[im，選項]）
%輸入：
%影象：輸入影象的雙精度或整數矩陣。
%預設為區域板影象。
%選項：演示的選項。預設值為“自動”
%‘auto’----自動演示，無輸入

function coeffs = decdemo( im, option )
disp('Welcome to the contourlet decomposition demo! :)');
disp('Type help decdemo for help' ) ;
disp('You can also view decdemo.m for details.'
 
) ;
disp(' ');

% Input image
if ~exist('im', 'var')
    % Zoneplate image: good for illustrating multiscale and directional
    % decomposition
    im = imread ('zoneplate.png') ;
end

% Show the input image
disp( 'Displaying the input image...');
clf;
imagesc(im, [0, 255]);
title( 'Input image' ) ;
 

axis image off;
colormap(gray);
input( 'Press Enter key to continue...' ) ;
disp( ' ' );


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Image decomposition by contourlets using the 
% pyramidal directional filter bank (PDFB).
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Parameteters:
nlevels = [0, 2, 3, 4] ;        % Decomposition level%分解層次
pfilter = 'pkva' ;              % Pyramidal filter%金字塔過濾器
dfilter = 'pkva' ;              % Directional filter%方向濾波器

% Contourlet decomposition：Contraclet分解
coeffs = pdfbdec( double(im), pfilter, dfilter, nlevels );

% Display the coefficients：顯示係數
disp('Displaying the contourlet coefficients...') ;
imcoeff = showpdfb( coeffs ) ;
title('Contourlet coefficients');
input('Press Enter key to continue...' ) ;
disp(' ');


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Pyramidal directional filter bank (PDFB) reconstruction.
% This is the inverse of pdfbdec, i.e.
% imrec = pdfbrec(coeffs, pfilter, dfilter);
% would reconstruct imrec = im
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Reconstruct image
imrec = pdfbrec( coeffs, pfilter, dfilter ) ;

disp('Displaying the reconstructed image...') ;
disp('It should be a perfect reconstruction' ) ;
disp(' ') ;

% Show the reconstruction image and the original image
subplot(1,2,1), imagesc( im, [0, 255] ); 
title('Original image' ) ;
axis image off;
subplot(1,2,2), imagesc( imrec, [0, 255] );
title('Reconstructed image' ) ;
axis image off;

mse = sum( sum( (imrec - double(im)).^2 ) );
mse = mse / prod(size(im));

disp( sprintf('The mean square error is: %f', mse ) );
disp(' ');

金字塔分解：金字塔方向濾波器組（或輪廓線）分解
%y=pdfbdec（x，pfilt，dfilt，nlevs）
%
%輸入：
%x：輸入影象
%pfilt：金字塔分解步驟的過濾器名稱
%dfilt：定向分解步驟的過濾器名稱
%nlevs：方向濾波器組分解層數的向量
%在每個金字塔層面（從粗到細）。
%當階數為0時，一個臨界取樣的二維小波
%執行分解步驟。
%
%輸出：
%y：長度為（nlevs）+1的單元向量，其中除了y{1}是
%低通子帶，每個單元對應一個金字塔
%level和是一個包含帶通方向的單元向量
%來自該級別的DFB的子帶。
%
%索引約定：
%假設························。
%對於1，則k=j
%y{J+2-J}{k}（n_1，n_2）
%是標度2^j、方向k和位置處的輪廓線係數
%（n_12^{（j+l_j-2），n_2*2}j）對於k<=2^（l_j-1），則，
%（n_12^j，n_2*2（j+l_j-2））對於k>2^（l_j-1）。
%當k從1增加到2^l
%角135度，在cotan中均勻遞增，從-1到1
%k<=2^（l_j-1），然後在tan中均勻遞減，從1到-1
%k>2^（l_j-1）。
%
%另請參見：PFILTERS、DFILTERS、PDFBRE
下面展示一些 內聯程式碼片。

function y = pdfbdec(x, pfilt, dfilt, nlevs)
if length(nlevs) == 0
    y = {x};
    
else
    % Get the pyramidal filters from the filter name
    [h, g] = pfilters(pfilt);
    
    if nlevs(end) ~= 0
        % Laplacian decomposition
        [xlo, xhi] = lpdec(x, h, g);
    
        % DFB on the bandpass image
        switch dfilt        % Decide the method based on the filter name
            case {'pkva6', 'pkva8', 'pkva12', 'pkva'}   
                % Use the ladder structure (whihc is much more efficient)
                xhi_dir = dfbdec_l(xhi, dfilt, nlevs(end));
            
            otherwise       
                % General case
                xhi_dir = dfbdec(xhi, dfilt, nlevs(end));                
        end
        
    else        
        % Special case: nlevs(end) == 0
        % Perform one-level 2-D critically sampled wavelet filter bank
        [xlo, xLH, xHL, xHH] = wfb2dec(x, h, g);
        xhi_dir = {xLH, xHL, xHH};
    end
    
    % Recursive call on the low band
    ylo = pdfbdec(xlo, pfilt, dfilt, nlevs(1:end-1));

    % Add bandpass directional subbands to the final output
    y = {ylo{:}, xhi_dir};
end
%%%%%%%%%%%%%%%%%%%%%%%%%pfilter子函式%%%%%%%%%%%%%%%%%%%%%
function [h, g] = pfilters(fname)
% PFILTERS    Generate filters for the Laplacian pyramid
%pfilter為拉普拉斯金字塔生成濾波器
%	[h, g] = pfilters(fname)
%
% Input:
%	fname:	Name of the filters, including the famous '9-7' filters
%		    and all other available from WFILTERS in Wavelet toolbox
%
% Output:
%	h, g:	1D filters (lowpass for analysis and synthesis, respectively)
%		    for seperable pyramid

switch fname
    case {'9-7', '9/7'}
	    h = [.037828455506995 -.023849465019380 -.11062440441842 ...
	         .37740285561265];	
	    h = [h, .85269867900940, fliplr(h)];
    
	    g = [-.064538882628938 -.040689417609558 .41809227322221];
	    g = [g, .78848561640566, fliplr(g)];
        
    case {'5-3', '5/3'}
	    h = [-1, 2, 6, 2, -1] / (4 * sqrt(2));
        g = [1, 2, 1] / (2 * sqrt(2));
        
    case {'Burt'}
	    h = [0.6, 0.25, -0.05];
	    h = sqrt(2) * [h(end:-1:2), h];
	
	    g = [17/28, 73/280, -3/56, -3/280];
	    g = sqrt(2) * [g(end:-1:2), g];
 	    
    case {'pkva'}	% filters from the ladder structure	
	    % Allpass filter for the ladder structure network
	    beta = ldfilter(fname);
	
	    lf = length(beta);
	    n = lf / 2;
	
	    if n ~= floor(n)
	        error('The input allpass filter must be even length');
        end
	
	    % beta(z^2)
	    beta2 = zeros(1, 2*lf-1);
	    beta2(1:2:end) = beta;
	
	    % H(z)
	    h = beta2;
	    h(2*n) = h(2*n) + 1;
	    h = h / 2;
	
	    % G(z)
	    g = -conv(beta2, h);
	    g(4*n-1) = g(4*n-1) + 1;
	    g(2:2:end) = -g(2:2:end);
	
	    % Normalize
	    h = h * sqrt(2);
	    g = g * sqrt(2);

    otherwise
	    [h, g] = wfilters(fname, 'l');
end

lpdec子函式內聯程式碼片`。

function [c, d] = lpdec(x, h, g)
% LPDEC   Laplacian Pyramid Decomposition
%拉普拉斯金字塔分解
%	[c, d] = lpdec(x, h, g)
%
% Input:
%   x:      input image
%   h, g:   two lowpass filters for the Laplacian pyramid
%
% Output:
%   c:      coarse image at half size
%   d:      detail image at full size
%
% See also:	LPREC, PDFBDEC

% Lowpass filter and downsample
xlo = sefilter2(x, h, h, 'per');
c = xlo(1:2:end, 1:2:end);    
    
% Compute the residual (bandpass) image by upsample, filter, and subtract
% Even size filter needs to be adjusted to obtain perfect reconstruction
adjust = mod(length(g) + 1, 2);

xlo = zeros(size(x));
xlo(1:2:end, 1:2:end) = c;
d = x - sefilter2(xlo, g, g, 'per', adjust * [1, 1]);

PDFB：%SHOWPDFB顯示PDFB係數。
%顯示PDFB系統
%showpdfb（y，[縮放模式，顯示模式。。。
%低比率、高比率、係數模式、子帶隙]）
%輸入：
%y：長度為n+1的細胞向量，每層一個
%來自DFB的子帶影象，y{1}是低通影象
%縮放模式：
%縮放模式（字串或數字）：
%如果是一個數，則表示最顯著的
%要顯示的係數。其預設值為“auto2”。
%“auto1”–所有圖層使用一個比例。它反映了真實的
%係數的值。
%然而，能見度會很差。
%“auto2”–低通使用第一個刻度。所有高通量的使用
%第二個音階。
%“auto3”–低通使用第一個刻度。
%所有的小波高通量都使用第二個尺度。
%所有的輪廓線高通量使用第三個刻度。
%顯示模式：
%顯示模式（字串）：
%“matlab”–在matlab環境下顯示。
%它使用背景色作為邊緣
%影象。
%“其他”–在其他環境中顯示或列印。
%它使用了白色的邊緣影象。
%它是預設值。
%低比率：
%低通濾波器的顯示比率（預設值為2）。
%範圍從1.2到4.0。
%高比率：
%高通濾波器的顯示比率（預設值為6）。
%範圍從1.5到10。
%共模：
%係數模式（字串）：
%‘real’——高通濾波器使用實係數。
%‘abs’——高通濾波器使用絕對係數。
%它是預設值
%子帶間距：
%子帶之間的間距（以畫素為單位）。範圍從1到4。
%輸出：
%displayIm：顯示影象的矩陣。
%另請參見：PDFBDEC、DFBIMAGE、COMPUTESCALE
下面展示一些 內聯程式碼片。

function displayIm = showpdfb(y, scaleMode, displayMode, ...
                              lowratio, highratio, coefMode, subbandgap)
if ~exist('scaleMode', 'var')
    scaleMode = 'auto2' ;
elseif isnumeric( scaleMode ) % Denote the number of significant coefficients to be displayed
    if scaleMode < 2
        display( 'Warning! The numbe of significant coefficients must be positive!' ) ;
        scaleMode = 50 ;
    end
elseif ~strcmp(scaleMode,'auto1') & ~strcmp(scaleMode, 'auto2') & ~strcmp(scaleMode, 'auto3')
    display ('Warning! There are only two scaleMode mode: auto1, auto2, auto3! Its defualt value is "auto2"!');
    scaleMode = 'auto2' ;
end

% Display ratio for the lowpass band
if ~exist('lowratio', 'var')
    lowratio = 2 ;
elseif highratio < 1
    display ('Warning! lowratio must be larger than 1!Its defualt value is 2!');
end

% Display ratio for the hiphpass band
if ~exist('highratio', 'var')
    highratio = 6 ;
elseif highratio < 1
    display ('Warning! highratio must be larger than 1! Its defualt value is 6!');
end

% Gap between subbands
if ~exist('subbandgap', 'var')
    subbandgap = 1;  
elseif subbandgap < 1
    display ('Warning! subbandgap must be no less than 1! Its defualt value is 1!');
    subbandgap = 1;
end

% Display mode
if ~exist('displayMode', 'var')
    displayMode = 'others' ;
elseif ~strcmp(displayMode,'others') & ~strcmp(displayMode, 'matlab')
    display ('Warning! There are only two display mode: matlab, others! Its defualt value is "others"!');
    displayMode = 'others' ;
end

% Coefficient mode
if ~exist('coefMode', 'var')
    coefMode = 'abs' ;
elseif ~strcmp(coefMode,'real') & ~strcmp(coefMode, 'abs')
    display ('Warning! There are only two coefficients mode: real, abs! Its defualt value is "abs"!');
    coefMode = 'abs' ;
end

% Parameters for display
layergap = 1 ; % Gap between layers

% Input structure analysis. 
nLayers = length(y); %number of PDFB layers
% Compute the number of wavelets layers.
% We assume that the wavelets layers are first several consecutive layers.
% The number of the subbands of each layer is 3.
fWaveletsLayer = 1;
nWaveletsLayers = 0 ; %Number of wavelets layers.
nInxContourletLayer = 0 ; %The index of the first contourlet layer.
i = 2 ;
while fWaveletsLayer > 0 & i <= nLayers
    if length( y{i} ) == 3
        nWaveletsLayers = nWaveletsLayers + 1 ;
    else
        fWaveletsLayer = 0 ;
    end;
    i = i + 1 ;
end;
nInxContourletLayer = 2 + nWaveletsLayers ;
    
% Initialization 
% Since we will merge the wavelets layers together, 
% we shall decrease the number of display layers.
nDisplayLayers = nLayers - nWaveletsLayers ;
cellLayers = cell (1, nDisplayLayers); % Cell for multiple display layers
vScalesTemp = zeros (1, 2) ;  % Temporary scale vector.
vScales = zeros (nLayers, 2); % Scale vectors for each layer 
nAdjustHighpass = 2 ; % Adjustment ratio for the highpass layers.


if ~isnumeric( scaleMode )
    switch( scaleMode )% Compute the scales for each layer
    case 'auto1'
        vScalesTemp = computescale( y, lowratio, 1, nLayers, coefMode ) ;
        for i = 1 : nLayers
            vScales( i, :) = vScalesTemp ;
        end
    case 'auto2'
        vScales( 1, :) = computescale( y, lowratio, 1, 1, coefMode ) ;
        vScalesTemp = computescale( y, highratio, 2, nLayers, coefMode ) ;
        % Make a slight adjustment. Compared to the lowpass, the highpass shall be insignificant.
        % To make the display more realistic, use a little trick to make the upper bound a little bigger. 
        vScalesTemp (2) = nAdjustHighpass * ( vScalesTemp (2) - vScalesTemp (1) ) + vScalesTemp (1) ;
        for i = 2 : nLayers
            vScales( i, :) = vScalesTemp ;
        end
    case 'auto3'
        vScales( 1, :) = computescale( y, lowratio, 1, 1, coefMode ) ;
        vScalesTemp = computescale( y, highratio, 2, 1+nWaveletsLayers, coefMode ) ;
        % Make a slight adjustment. Compared to the lowpass, the highpass shall be insignificant.
        % To make the display more realistic, use a little trick to make the upper bound a little bigger. 
        vScalesTemp (2) = nAdjustHighpass * ( vScalesTemp (2) - vScalesTemp (1) ) + vScalesTemp (1) ;
        for i = 2 : nWaveletsLayers + 1
            vScales( i, :) = vScalesTemp ;
        end
        vScalesTemp = computescale( y, highratio, nInxContourletLayer, nLayers, coefMode ) ;
        % Make a slight adjustment. Compared to the lowpass, the highpass shall be insignificant.
        % To make the display more realistic, use a little trick to make the upper bound a little bigger. 
        vScalesTemp (2) = nAdjustHighpass * ( vScalesTemp (2) - vScalesTemp (1) ) + vScalesTemp (1) ;
        for i = nInxContourletLayer : nLayers
            vScales( i, :) = vScalesTemp ;
        end
    otherwise % Default value: 'auto2'.
        vScales( 1, :) = computescale( y, lowratio, 1, 1, coefMode ) ;
        vScalesTemp = computescale( y, highratio, 2, nLayers, coefMode ) ;
        for i = 2 : nLayers
            vScales( i, :) = vScalesTemp ;
        end
    end
    % Verify that they are reasonable
    for i = 1: nLayers
        if vScales (i, 2) < vScales (i, 1) + 1.0e-9
            display ('Error! The scale vectors are wrong! Exit!' ) ;
            exit ;
        end;
    end;
    %display ( vScales ) ;
else % Compute the threshold for the display of coefficients
    % Convert the output into the vector format
    [vCoeff, s] = pdfb2vec(y);
    
    % Sort the coefficient in the order of energy.
    vSort = sort( abs( vCoeff ));
    clear vCoeff;
    vSort = fliplr(vSort);
    
    % Find the threshold value based on number of keeping coeffs
    dThresh = vSort( scaleMode );
    clear vSort;
end

% Prepare for the display
colormap(gray);
cmap = get(gcf,'Colormap');
cColorInx = size(cmap,1);

%  Find background color index:
if strcmp( displayMode, 'matlab' )
    % Get the background color (gray value)
    dBgColor = get( gcf, 'Color' ) ;
    
    % Search the color index by 2-fold searching method.
    % This method is only useful for the gray color!
    nSmall = 1 ;
    nBig = cColorInx ;
    while nBig > nSmall + 1
        nBgColor = floor ((nSmall + nBig) / 2) ;
        if dBgColor(1) < cmap (nBgColor, 1)
            nBig = nBgColor ;
        else
            nSmall = nBgColor ;
        end
    end;
    if abs( dBgColor(1) - cmap (nBig, 1) ) > abs ( dBgColor(1) - cmap( nSmall, 1) )
        nBgColor = nSmall ;
    else
        nBgColor = nBig ;
    end
end

% Merge all layers to corresponding display layers.
% Prepare the cellLayers, including the boundary.
% Need to polish with real boudary later!
% Now we add the boundary, but erase the images!!
% First handle the lowpass filter
% White line around subbands

% 1. One wavelets layers. 
gridI = cColorInx - 1 ;
cell4Wavelets = cell(1, 4) ; %Store 4 wavelets subbands.
if isnumeric ( scaleMode ) %Keep the significant efficients
    waveletsIm = cColorInx * double(abs(y{1}) >= dThresh) ;
else
    dRatio = (cColorInx-1) / (vScales(1,2)-vScales(1,1));
    if strcmp( coefMode, 'real' )
        waveletsIm = double( 1 + (y{1}-vScales(1,1))*dRatio ) ;
    else
        waveletsIm = double( 1 + (abs(y{1})-vScales(1,1))*dRatio ) ;
    end
end
% Merge other wavelets layers
if nWaveletsLayers > 0
    for i=2 : nWaveletsLayers + 1
        cell4Wavelets{1} = waveletsIm ; 
        % Compute with the scale ratio.
        if ~isnumeric( scaleMode )
            dRatio = (cColorInx-1) / (vScales(i,2)-vScales(i,1));
        end
        m = length(y{i});
        if m ~= 3
            display('Error! Incorect number of wavelets subbands! Exit!');
            exit;  
        end
        for k = 1:m
            if isnumeric ( scaleMode ) %Keep the significant efficients
                cell4Wavelets{k+1} = cColorInx * double(abs(y{i}{k}) >= dThresh) ;
            else
                if strcmp( coefMode, 'real' )
                    cell4Wavelets{k+1} = double( 1 + (y{i}{k}-vScales(i,1))*dRatio );
                else
                    cell4Wavelets{k+1} = double( 1 + (abs(y{i}{k})-vScales(i,1) )* dRatio );
                end
            end
        end
        waveletsIm = dfbimage(cell4Wavelets, subbandgap, gridI);
    end
end
cellLayers{1} = waveletsIm ;
nHeight = size( cellLayers{1}, 1 );


% 2. All the contourlet layers
for i = nInxContourletLayer : nLayers
    % Compute with the scale ratio.
    if ~isnumeric( scaleMode )
        dRatio = (cColorInx-1) / (vScales( i, 2)-vScales( i, 1));
    end
    m = length(y{i});
    z = cell(1, m);
    for k = 1:m
        if isnumeric ( scaleMode ) %Keep the significant efficients
            z{k} = cColorInx * double(abs(y{i}{k}) >= dThresh) ;
        else
            if strcmp( coefMode, 'real' )
                z{k} = double( 1 + (y{i}{k}-vScales(i,1)) * dRatio );
            else
                z{k} = double( 1 + (abs(y{i}{k})-vScales(i,1) )* dRatio );
            end
        end
    end
    cellLayers{i-nWaveletsLayers} = dfbimage(z, subbandgap, gridI);
    nHeight = nHeight + size(cellLayers{i-nWaveletsLayers}, 1);
end
% Compute the width of the dispaly image.
nWidth = size(cellLayers{nDisplayLayers}, 2);

% Merge all layers and add gaps between layers
nHeight = nHeight + layergap * (nDisplayLayers - 1) ;
% Set the background for the output image
if strcmp( displayMode, 'matlab' )
    displayIm = nBgColor * ones( nHeight, nWidth);
else
    displayIm = (cColorInx-1) * ones( nHeight, nWidth);
end
nPos = 0; %output image pointer
for i = 1 : nDisplayLayers
    [h, w] = size( cellLayers{i} );
    displayIm( nPos+1: nPos+h, 1:w) = cellLayers{i};
    if i < nDisplayLayers
        % Move the position pointer and add gaps between layers
        nPos = nPos + h + layergap ;
    end        
end

hh = image( displayIm );
% title('decompostion image');
axis image off

%PDFBREC金字塔定向濾波器組重構
%x=pdfbrec（y，pfilt，dfilt）
%輸入：
%y：長度為n+1的細胞向量，每層一個
%來自DFB的子帶影象，y{1}是低頻帶影象
%pfilt：金字塔的過濾器名稱
%dfilt：定向過濾器組的過濾器名稱
%輸出：
%x：重建影象
%另請參見：PFILTERS、DFILTERS、PDFBDEC

下面展示一些 內聯程式碼片。

function x = pdfbrec(y, pfilt, dfilt)
n = length(y) - 1;
if n <= 0
    x = y{1};
    
else
    % Recursive call to reconstruct the low band
    xlo = pdfbrec(y(1:end-1), pfilt, dfilt);
    
    % Get the pyramidal filters from the filter name
    [h, g] = pfilters(pfilt);
    
    % Process the detail subbands
    if length(y{end}) ~= 3
        % Reconstruct the bandpass image from DFB
        
        % Decide the method based on the filter name
        switch dfilt        
            case {'pkva6', 'pkva8', 'pkva12', 'pkva'}	
                % Use the ladder structure (much more efficient)
                xhi = dfbrec_l(y{end}, dfilt);
                
            otherwise	
                % General case
                xhi = dfbrec(y{end}, dfilt); 
        end
        
        x = lprec(xlo, xhi, h, g);
   
    else    
        % Special case: length(y{end}) == 3
        % Perform one-level 2-D critically sampled wavelet filter bank
        x = wfb2rec(xlo, y{end}{1}, y{end}{2}, y{end}{3}, h, g);
    end
end

時間：2020年12月17日15:56
註釋：部分子函式篇幅問題不能一一展現，網路資源很多，可自行下載。http://www.matlabsky.com/thread-879-1-3.html