spiht_wpackets.m

zerotree的寫法,給大家參考

            bitstream(1,count1)=any(polje>=quantstep);
            bitstream(2,count1)=pass;
            bitstream(3,count1)=1; %position bit
            %fwrite(fileout,sn,'ubit1');%f
            if bitstream(1,count1) %tree is significant 
                polje=D(id+(jd-1)*Drows);
                checkbit=0; %indicates if significant coefficients is found between children
                for k=1:length(polje); 
                    count1=count1+1;
                    bitstream(1,count1)=(abs(polje(k))>=quantstep);
                    bitstream(2,count1)=pass;
                    bitstream(3,count1)=1; %position bit
                    if bitstream(1,count1)
                        %fwrite(fileout,1,'ubit1');%f
                        %fwrite(fileout,sign(polje(k)),'ubit1');%f
                        if ~checkbit checkbit=1;end;
                        count1=count1+1;
                        bitstream(1,count1)=(sign(polje(k))==1);
                        bitstream(2,count1)=pass;
                        bitstream(3,count1)=2; %sign bit
                        Drec(id(k),jd(k))=sign(double(bitstream(1,count1))-0.5)*init;
                        %RecEnergy=RecEnergy+RecEadd;
                        endLSP=endLSP+1;
                        LSP(endLSP,:)=[id(k) jd(k)];
                        %fprintf(fid,'%5.2f\n',D(id(k),jd(k)));  %**STATS**
                    else
                        %fwrite(fileout,0,'ubit1');%f
                        endLIP=endLIP+1;
                        LIP(endLIP,:)=[id(k) jd(k)];
                    end;
                    if count1+1>=numbits(cntbitout) %this limit is equivalent in terms of distortion (1 addidional bit wouldn't help)
                        [Arec,PSNR(cntbitout),MSE(cntbitout)]=reconstruct(A,Drec,cntbitout,count1,quantstep,pass,param,p_stream);
                        cntbitout=cntbitout+1;
                    end;
                    if cntbitout>length(numbits) %| abs(RecEnergy-Energy)/Energy<10^-7
                        [stats,fid]=show_stats(fid,'LIS',pass,count1,stats,endLSP); %**STATS**
                        [stats,fid]=show_stats(fid,'LSP',pass,stats(pass,1).LSP-1,stats,endLSP);
                        [stats,fid]=show_stats(fid,'total',pass,count1,stats,endLSP);
                        fprintf('\n...LIS (checking children) loop out');
                        %fclose(fileout);
                        return;
                    end;
                end; 
                if (i<=Dqrows) & (j<=Dqcols)  %if L(i,j) (grandchildren) exist, (i,j) goes at the end of LIS as type D
                    if checkbit
                        endLIS=endLIS+1;
                        LIS(endLIS,:)=[i j Dtype];
                    else %if no child is found significant then some of descendents must be significant
                        tn=length(polje);
                        LIS(endLIS+1:endLIS+tn,1:3)=[id',jd',zeros(tn,1)];
                        endLIS=endLIS+tn;
                    end;  
                end;    
            else  %back to the LIS goes insignificant element  
                liscnt=liscnt+1;
                LISnew(liscnt,:)=LISel; 
            end;
        else %2.2.2. for D type elements
            polje=MagnD(id+(jd-1)*Drows); %children
            count1=count1+1; 
            bitstream(1,count1)=any(polje>=quantstep);
            bitstream(2,count1)=pass;
            bitstream(3,count1)=1; %position bit
            %fwrite(fileout,bitstream(1,count1),'ubit1');%f
            if bitstream(1,count1) %if D(id,jd) is significant then elements O(id,jd) go to LIS as type L    
                endLIStemp=length(polje);
                LIStemp(1:endLIStemp,1:3)=[id',jd',zeros(endLIStemp,1)];
            else %if it's not it stays in LIS  
                liscnt=liscnt+1;
                LISnew(liscnt,:)=LISel; 
            end;  
        end; 
        if count1+1>=numbits(cntbitout) %this limit is equivalent in terms of distortion (1 addidional bit wouldn't help)
            [Arec,PSNR(cntbitout),MSE(cntbitout)]=reconstruct(A,Drec,cntbitout,count1,quantstep,pass,param,p_stream);
            cntbitout=cntbitout+1;
        end;
        if cntbitout>length(numbits) %| abs(RecEnergy-Energy)/Energy<10^-7
            [stats,fid]=show_stats(fid,'LIS',pass,count1,stats,endLSP); %**STATS**
            [stats,fid]=show_stats(fid,'LSP',pass,stats(pass,1).LSP-1,stats,endLSP);
            [stats,fid]=show_stats(fid,'total',pass,count1,stats,endLSP);
            fprintf('\n...LIS loop out');
            %fclose(fileout);
            return;
        end;
    end; %end of LIS loop
    endLIS=liscnt; 
    LIS(1:endLIS,:)=LISnew(1:endLIS,:);
    [stats,fid]=show_stats(fid,'LIS',pass,count1,stats,endLSP); %**STATS**
    %3. REFINEMENT PASS
    if pass>1
        refine=quantstep/2;
        for mls=1:oldendLSP
            i=LSP(mls,1);
            j=LSP(mls,2);
            %RecEnergy=RecEnergy-Drec(i,j).^2;
            count1=count1+1; 
            bitstream(1,count1)=(D(i,j)>Drec(i,j));
            bitstream(2,count1)=pass;
            bitstream(3,count1)=3; %refine bit
            if bitstream(1,count1)  
                Drec(i,j)=Drec(i,j)+refine; 
            else
                Drec(i,j)=Drec(i,j)-refine; 
            end; 
            %RecEnergy=RecEnergy+Drec(i,j).^2;
            if count1>=numbits(cntbitout) 
                [Arec,PSNR(cntbitout),MSE(cntbitout)]=reconstruct(A,Drec,cntbitout,count1,quantstep,pass,param,p_stream);
                cntbitout=cntbitout+1;
            end;
            if cntbitout>length(numbits) %| abs(RecEnergy-Energy)/Energy<10^-7
                [stats,fid]=show_stats(fid,'LSP',pass,count1,stats,endLSP); %**STATS**
                [stats,fid]=show_stats(fid,'total',pass,count1,stats,endLSP);
                fprintf('Refinement pass out...\n');
                return;
            end;
        end;
    end;
    
    [stats,fid]=show_stats(fid,'LSP',pass,count1,stats,endLSP); %**STATS**
    %4. UPDATE
    oldendLSP=endLSP;
    quantstep=quantstep/2;
    pass=pass+1;
end; %end of the main loop
fprintf('Smallest quantization step reached!!\n');

function [LIP,endLIP,LIS,endLIS,LSP,endLSP]=initialize_lists(D,N)
[Drows,Dcols]=size(D);
sizrow=Drows/2^N;
sizcol=Dcols/2^N;
duzina=numel(D);
LIP=zeros(duzina,2);
LIS=zeros(duzina,3);
LSP=zeros(duzina,2); 
Dp=D(1:sizrow,1:sizcol).'; %lowest subband
endLIP=sizrow*sizcol; %LIP contains all lowest subband coefficients
%for EZW type trees, LIS will contain all LIP coefficients, raster scan
ii=repmat(1:sizrow,sizcol,1);
jj=repmat(1:sizcol,[1 sizrow],1);
LIP(1:endLIP,1)=ii(:);
LIP(1:endLIP,2)=jj.';
endLIS=endLIP;
LIS(1:endLIS,1:2)=LIP(1:endLIP,1:2);
endLSP=0; %LSP is initialized empty

function [stats,fid]=show_stats(fid,looptype,pass,count1,stats,endLSP)
stats_index=pass+1;
addstr=' ';
switch looptype
    case 'LIP'    
        stats(stats_index,1).LIP=endLSP-stats(pass,1).LSP;
        stats(stats_index,2).LIP=count1-stats(pass,2).LSP;
        if stats(stats_index,1).LIP
            bpc=stats(stats_index,2).LIP/stats(stats_index,1).LIP;
            if bpc>9 addstr='';end;
        else 
            bpc=Inf;
            addstr='  ';
        end; 
        fprintf(fid,'%6d/%-6d=%s%2.2f bpc',stats(stats_index,2).LIP,stats(stats_index,1).LIP,addstr,bpc);
    case 'LIS'
        stats(stats_index,1).LIS=endLSP-stats(stats_index,1).LIP-stats(pass,1).LSP;
        stats(stats_index,2).LIS=count1-stats(stats_index,2).LIP-stats(pass,2).LSP;
        if stats(stats_index,1).LIS
            bpc=stats(stats_index,2).LIS/stats(stats_index,1).LIS;
            if bpc>9 addstr='';end;
        else
            bpc=Inf;

            addstr='  ';
        end;
        fprintf(fid,'%7d/%-6d=%s%2.2f bpc',stats(stats_index,2).LIS,stats(stats_index,1).LIS,addstr,bpc);
    case 'LSP'
        stats(stats_index,1).LSP=endLSP;
        stats(stats_index,2).LSP=count1;
%       pass_bits=stats(stats_index,2).LSP-stats(pass,2).LSP;
        fprintf(fid,'%12d(%8d)\n',endLSP,count1- stats(pass,1).LSP+1);
    case 'init'
        stats(1,1).LIP=0;      %elements added to LSP list
        stats(1,2).LIP=0;      %bits used
        stats(1,1).LIS=0;
        stats(1,2).LIS=0;
        stats(1,1).LSP=0;      %elements added to LSP list 
        stats(1,2).LSP=count1; %bits used
        %fprintf(fid,'*bpc=bits per significant coefficient\n');
        fprintf(fid,'                   LIP loop                LIS loop              Total significant\n');
        fprintf(fid,'pass(start bit)   bits/coeffs.            bits/coeffs.              coeffs.(LSP start bit)\n');
    case 'total'
        fprintf(fid,'Total:\n');
        totLIPbits=0;
        totLIPcoeffs=0;
        totLISbits=0;
        totLIScoeffs=0;
        for i=1:stats_index
            totLIPbits=totLIPbits+stats(i,2).LIP;
            totLIPcoeffs=totLIPcoeffs+stats(i,1).LIP;
            totLISbits=totLISbits+stats(i,2).LIS;
            totLIScoeffs=totLIScoeffs+stats(i,1).LIS;
        end; 
        if totLIPcoeffs
            bpc=totLIPbits/totLIPcoeffs;
            if bpc>9 addstr='';end;
        else
            bpc=Inf;
            addstr='  ';
        end;
        fprintf(fid,'                %6d/%-6d=%s%2.2f bpc',totLISbits,totLIScoeffs,addstr,bpc);    
        if totLIScoeffs
            bpc=totLISbits/totLIScoeffs;
            if bpc>9 addstr='';end;
        else
            bpc=Inf;
            addstr='  ';
        end;
        fprintf(fid,'%7d/%-6d=%s%2.2f bpc',totLIPbits,totLIScoeffs,addstr,bpc);
        fprintf(fid,'%12d\n',endLSP);
        fclose(fid);
end;


function [Arec,PSNR,MSE]=reconstruct(A,Drec,i,counts,quantsteps,passes,param,p_stream)
Asiz=prod(size(Drec));
fprintf('\n%d. bitrate - %f bpp\n',i,counts/Asiz);
% optimization in regard to probability density function (pdf) of coefficients 
Drec=pdf_opt(Drec,quantsteps); %
%mean value shifting as in original SPIHT /begin
%DLLr=size(A,1)/2^param.N;DLLc=size(A,2)/2^param.N;
%Drec(1:DLLr,1:DLLc)=Drec(1:DLLr,1:DLLc)+param.Dmean*4^param.shift;
%mean value shifting as in original SPIHT /end
%perform reconstruction
Arec=recon_packets2D(Drec,param,p_stream);
Arec=Arec+param.meanvalue;
Arec(Arec>255)=255;
Arec(Arec<0)=0;
Arec=uint8(round(Arec));
%imwrite(Arec,['RecImage_' num2str(i) '.bmp']);
%display some statistics
filebits=8*ceil(counts/8);
fprintf('No. output bits - %d(%d bytes; %d excess bits)-> %f bpp\n',...
    filebits,filebits/8,filebits-counts,filebits/Asiz);
%compute PSNR
[MSE,PSNR]=iq_measures(A,Arec);

function [M,s]=maximum_magnitudes(D,s)
%Input:
% D - wavelet coefficients of the original image
% s - structure containing info on parent-children relationship between
% subbands (see in decomp_packets2D.m)
%
%Output:
% M - contains maximal values within the set composed of a coefficient 
%     and all of its descendants. M(:,:,2) contains subband index.
% s - modified input structure s. The following fields are added: addrows, 
%     addcols and scalediff.
%
%Note:
% This function sorts coefficients trees by amplitudes. Modifies the
% structure by adding additional info used by the SPIHT algorithm.

[Drows,Dcols]=size(D);
M=zeros(Drows,Dcols,3);
M(:,:,1:2)=repmat(abs(D),[1 1 2]);
M(1:s(1).band_abs(4),1:s(1).band_abs(3),3)=1; %the lowest subband coefficients are in subband no.1 (in 's') 
s(1).addrows=[]; %add field addrows
s(1).addcols=[]; %add field addcols
s(1).scalediff=[]; %add field scalediff
[M,s]=sort_subbands(1,s,M);

function [M,s]=sort_subbands(index,s,M)
node=s(index);
col=node.band_abs(1);
row=node.band_abs(2);
cw=node.band_abs(3);
rw=node.band_abs(4);
magn_band=zeros(rw,cw);%
children=node.children;
for i=1:size(children,2)
    child=s(children(i));
    col_c=child.band_abs(1);
    row_c=child.band_abs(2);
    cw_c=child.band_abs(3);
    rw_c=child.band_abs(4);
    M(row_c+1:row_c+rw_c,col_c+1:col_c+cw_c,3)=children(i); %index of the subband
    if child.children(1)>0
        [M,s]=sort_subbands(children(i),s,M);   
    end;    
    scalediff=child.scale-node.scale;
    chdim=2^scalediff;
    mr=repmat(1:chdim,[chdim 1]);
    mc=mr';
    diffrow=row_c-row;
    diffcol=col_c-col;
    s(index).addrows=[s(index).addrows (mr(:)+diffrow-1)'];
    s(index).addcols=[s(index).addcols (mc(:)+diffcol-1)'];
    s(index).scalediff=[s(index).scalediff scalediff*ones(1,numel(mr))];
    subbandM=reshape(M(row_c+1:row_c+rw_c,col_c+1:col_c+cw_c,1),[chdim rw_c/chdim chdim cw_c/chdim]);
    subbandM=permute(subbandM,[2 4 3 1]);
    subbandM=reshape(subbandM,[rw_c/chdim cw_c/chdim chdim^2]);
    subbandM=max(subbandM,[],3);
    magn_band=max(magn_band,subbandM);
end;   
M(row+1:row+rw,col+1:col+cw,2)=magn_band;
M(row+1:row+rw,col+1:col+cw,1)=max(M(row+1:row+rw,col+1:col+cw,1:2),[],3); %sort magnitudes for L