spiht_main.m

对视频序列进行三维小波变换

% reads the luminance (Y) component from the filename specified. The file must
% be of the format CIF, QCIF, SIF. It reads nFrame frames.

nFrame=8;                                               %输入要处理的帧数
filename='stefan.cif';      %输入要处理的视频序列

Y = myfunc_read_cif(nFrame,filename);                   %读取图像序列的亮度成分

%______ 3D - SPIHT encoder______________________        参数输入

 type = 'bior4.4';      % type of filters for the decomposition
 spatial_levels = 3;    % levels of decomposition

%type =input('the type of the wavelet is: ','s');         %输入小波变换类型
%spatial_levels=input('please input the level: ');        %输入空间变换级数
%fprintf('\n the size of the video sequence is %d %d %d \n',size(Y));         %输出视频序列大小
%bitrate=input('please input the bitrate,the unit is bpp: ');                 %输入比特率

%---------Wavelet Decomposition-----------------

nFrame1=nFrame;
cA1=zeros(size(Y,1),size(Y,2),nFrame1);
cD1=zeros(size(Y,1),size(Y,2),nFrame1);
Y1=zeros(size(Y,1),size(Y,2),nFrame1);
j=spatial_levels;

for lel_t=1:j                                                  %按照级数循环
 for m=1:size(Y,1)                                                      % N次时间轴小波变换
    for n=1:size(Y,2)
        [cA1(m,n,1:nFrame1/2),cD1(m,n,1:nFrame1/2)]=dwt(Y(m,n,1:nFrame1),'haar');
       for k=1:nFrame1/2
        Y1(m,n,k)=cA1(m,n,k);
       end
        i=1;
       for k=(nFrame1/2)+1:nFrame1 
        Y1(m,n,k)=cD1(m,n,i);
        i=i+1;
       end
    end
 end
    nFrame1=nFrame1/2;                                                  %时域分解后存Y1中
end                                                      

%———————————————————————————— 
S1=zeros(spatial_levels+2,3,spatial_levels);                      % n级空域分解
nFrame2=nFrame;

for lel_f=1:spatial_levels
    for i=1:nFrame2
       [video_seq(:,:,i) , S] = func_DWT(Y1(:,:,i),lel_f, type);
    end
    S1(1:lel_f+2,1:3,lel_f)=S;
    nFrame2=nFrame2/2;
    fprintf('%d \n',lel_f);
end          
   


figure;  

                                                      %时间轴小波变换与空间小波变换交替进
for i=1:nFrame                                       %显示变换后的各帧图像                空域编码后存video_seq中
  subplot(2,nFrame/2,i);   
  imshow(uint8(video_seq(:,:,i)));
end

%__________________________________________________________________  %3D_SPIHT阶段
% maximum bits for the 3D-SPIHT encoding

max_bits = nFrame*size(Y,1)*size(Y,2)*bitrate;      % eg: 8 frames 388*352 0.5

% max_bits = 675840; % for exaple if we have 16 frames each 240x352 and use 0.5 bits per pixel max_bits = 675840

out = func_SPIHT_Enc(video_seq, max_bits,spatial_levels);

% 3D - SPIHT decoder

in = func_SPIHT_Dec(out,size(Y));
video_seq2 = in;

%--------------Wavelet Reconstruction------------

nFrame3=nFrame/(2^spatial_levels);                                                %空域解码
nFrame4=nFrame;
                                                                                  %空域解码后存在Y_new中
for j=1:spatial_levels
  for i=1:nFrame4
     Y_new(:,:,i) = func_InvDWT(video_seq(:,:,i), S1(1:j+2,:,j), type, j);
  end
  nFrame4=nFrame4/2;
end


A1=zeros(size(Y,1),size(Y,2),nFrame);
D1=zeros(size(Y,1),size(Y,2),nFrame);
Y2=zeros(size(Y,1),size(Y,2),nFrame);


for re_lel=1:spatial_levels                                                    %按照级数循环
  for m=1:size(Y,1)                                                            %时域解码   重建之后存Y2
    for n=1:size(Y,2)
        A1(m,n,1:nFrame3*2)=upcoef('a',cA1(m,n,1:nFrame3),'haar',1,nFrame3*2);
        D1(m,n,1:nFrame3*2)=upcoef('d',cD1(m,n,1:nFrame3),'haar',1,nFrame3*2);
       for k=1:nFrame3*2
        Y2(m,n,k)=A1(m,n,k)+D1(m,n,k);
        cA1(m,n,k)=Y2(m,n,k);
       end
       if (nFrame3*2) < nFrame
         di=1;
       for k1=nFrame3*2+1:nFrame3*2*2
        cD1(m,n,di)=Y_new(m,n,k1);
        di=di+1;
       end
       end
    end
  end
 nFrame3=nFrame3*2;
end 


for i=1:nFrame
    subplot(2,nFrame/2,i);
    imshow(uint8(Y2(:,:,i)));
end


for i=1:nFrame
 Q = 255;
MSE = sum(sum((double(Y(:,:,i))-double(Y2(:,:,i))).^2))/288 / 352;
fprintf('The %d frame psnr performance is %.2f dB\n',i,10*log10(Q*Q/MSE));
end