4d_spiht_main.m

4-D SPHIT 压缩Matlab 程序

function 4D_SPIHT_main(nRow,nColumn,nV,nFrame)

tic;
nRow = 512;
nColumn = 512;
nV = 20;
nFrame = 16;

fprintf('-----------------read data---------------\n');
init_Y = myfunc_read_dicom(nV,nFrame);  %进行数据读取

fprintf('-----------------read data done---------------\n');

%---------------------4D-SPIHT begin---------------

type = 'bior4.4';    %waveler basis
spatial_levels = 3;  %wavelet decomposition level for z-dimension
t_levels = 3;        %wavelet decomposition level for t-dimension

%---------------------Wavelet Decomposition-----------
fprintf('---------------Wavelet Decomposition--------------\n ');
%在t-dimension上进行一维小波变换
new_Y = [];
IW1_Y = zeros(1:size(init_Y,1),1:size(init_Y,2),1:size(init_Y,3),1:size(init_Y,4));  %
for k = 1:size(init_Y,3)
    for i = 1:size(init_Y,1)
        for j = 1:size(init_Y,2)
            for t = 1:size(init_Y,4)
               new_Y = [new_Y,init_Y(i,j,k,t)];  %
            end
            [IW1_Y,L] = wavedec(new_Y,t_levels,type);  %
            new_Y = [];
        end
    end
end

%在z-dimension上进行一维小波变换?????

%在每个v内每一层进行二维小波变换
for t = 1:size(init_Y,4)
    for z = 1:size(init_Y,3)
        [IW2_Y,S] = func_DWT2(IW1_Y(:,:,z,t),spatial_levels,type);  %
    end
end
fprintf('---------------Wavelet Decomposition done--------------\n ');

%-------------------4D-SPIHT Encoding---------------------------
fprintf('-----------------Encoding------------------\n');
max_bits =41943040;  %for exaple if we have 16 frames each 512*512*20 and use 0.5 bits per pixel max_bits = 41943040;
img_enc = func_SPIHT_Enc(IW2_Y,max_bits,spatial_levels,t_levels);
fprintf('-----------------Encoding done------------------\n');

%--------------------4D-SPIHT Decoding------------------------
fprintf('-----------------Decoding-----------------\n');
img_dec = func_SPIHT_Dec(img_enc,size(init_Y));
fprintf('-----------------Decoding done-----------------\n');

%--------------------Wavelet Reconstruction----------------
fprintf('---------------Wavelet Reconstruction--------------\n ');
%在t-dimension上进行一维小波反变换
for k = 1:size(img_dec,3)
    for i = 1:size(img_dec,1)
        for j = 1:size(img_dec,2)
            for t = 1:size(img_dec,4)
               new_Y = [new_Y,img_dec(i,j,k,t)];  %
            end
            IW1_Y = waverec(new_Y,t_levels,type);  %
            new_Y = [];
        end
    end
end

%在z-dimension上进行一维小波变换?????

%在每个v内每一层进行二维小波变换
for t = 1:size(init_Y,4)
    for z = 1:size(init_Y,3)
        img_reconstruct = func_InvDWT2(IW1_Y(:,:,z,t),S,spatial_levels,type);  %
    end
end
fprintf('---------------Wavelet Reconstruction done--------------\n ');

Q = 255;
MSE = sum(sum(sum(sum(img_reconstruct - init_Y).^2))) / (nRow * nColumn * nV *nFrame);
fprintf('The psnr performance is %.2f db\n',10 * log10(Q*Q/MSE));
toc;