www.pudn.com > SPIHT(Matlab).rar > SPIHT_Enc.m
function out=SPIHT_Enc(m,level,rate)
% Matlab implementation of SPIHT (without Arithmatic coding stage)
%
% Encoder
%
% input: m : input image in wavelet domain
% level : wavelet decomposition level
% rate : bits per pixel
% output: out : bit stream
%
% Jing Tian
% Contact me : scuteejtian@hotmail.com
% This program is part of my undergraduate project in GuangZhou, P. R. China.
% April - July 1999
%----------- Initialization -----------------
m=double(m);
[nRow nColumn]=size(m);
OrigSize = nRow*nColumn;%size(Orig_I, 1);%经过小波变换后小波系数矩阵中元素的个数
max_bits = floor(rate * OrigSize);%SPIHT算法得到的码个数
bitctr = 0;
out = 2*ones(1,max_bits - 14);%initialization 2
n_max = floor(log2(double(max(max(abs(m))))));
Bits_Header = 0;
Bits_LSP = 0;
Bits_LIP = 0;
Bits_LIS = 0;
%----------- output bit stream header ----------------
% image size, number of bit plane, wavelet decomposition level should be
% written as bit stream header.
out(1,[1 2 3 4]) = [nRow n_max level rate];
%个人认为应该将24改为4
index = 5;%%%out向量前四个已经被占用，编码从第五个空间开始
Bits_Header = Bits_Header + 24;%%个人认为应该将24改为4
%----------- Initialize LIP, LSP, LIS ----------------
temp = [];
bandsize = 2.^(log2(size(m, 1)) - level + 1);%第level层小波系数的维数
temp1 = 1 : bandsize;
for i = 1 : bandsize
temp = [temp; temp1];
end
LIP(:, 1) = temp(:);
temp = temp';
LIP(:, 2) = temp(:);
LIS(:, 1) = LIP(:, 1);
LIS(:, 2) = LIP(:, 2);
LIS(:, 3) = zeros(length(LIP(:, 1)), 1);
pstart = 1;
pend = bandsize / 2;
for i = 1 : bandsize / 2
LIS(pstart : pend, :) = [];
pdel = pend - pstart + 1;%删去
pstart = pstart + bandsize - pdel;%pstart = pstart + bandsize/2
pend = pend + bandsize - pdel;%pend = pend + bandsize/2
end
LSP = [];
n = n_max;
%----------- coding ----------------
while(bitctr < max_bits)
% Sorting Pass
LIPtemp = LIP; temp = 0;
for i = 1:size(LIPtemp,1)%%%%表示LL层中总的元素个数%%%%%此循环用来扫描LIP中的数据 返回是行数
temp = temp+1;
if (bitctr + 1) >= max_bits
if (bitctr < max_bits)
out(length(out))=[];
end
return
end
if abs(m(LIPtemp(i,1),LIPtemp(i,2))) >= 2^n % 1: positive; 0: negative
out(index) = 1; bitctr = bitctr + 1;%%%%bitstr用来记录编码位数
index = index +1; Bits_LIP = Bits_LIP + 1;
sgn = m(LIPtemp(i,1),LIPtemp(i,2))>=0;
out(index) = sgn; bitctr = bitctr + 1;
index = index +1; Bits_LIP = Bits_LIP + 1;
LSP = [LSP; LIPtemp(i,:)];%%%%%将重要系数的坐标加到LSP的链尾
LIP(temp,:) = []; temp = temp - 1;%%%%在LIP中删去此项
else
out(index) = 0; bitctr = bitctr + 1;
index = index +1;Bits_LIP = Bits_LIP + 1;
end
end
LIStemp = LIS; temp = 0; i = 1;
while ( i <= size(LIStemp,1))%%%%%%%size(LIStemp）表示直系子孙的元素个数
temp = temp + 1;
if LIStemp(i,3) == 0%%%参照180行 %%在处理（i,j）的直系子孙序列
if bitctr >= max_bits
return
end
max_d = func_MyDescendant(LIStemp(i,1),LIStemp(i,2),LIStemp(i,3),m);
if max_d >= 2^n
out(index) = 1; bitctr = bitctr + 1;
index = index +1; Bits_LIS = Bits_LIS + 1;
x = LIStemp(i,1); y = LIStemp(i,2);
if (bitctr + 1) >= max_bits
if (bitctr < max_bits)
out(length(out))=[];
end
return
end
if abs(m(2*x-1,2*y-1)) >= 2^n%%%对(2*x-1,2*y-1)操作
LSP = [LSP; 2*x-1 2*y-1];%%%将重要系数坐标加入到LSP中
out(index) = 1; bitctr = bitctr + 1;
index = index +1; Bits_LIS = Bits_LIS + 1;
sgn = m(2*x-1,2*y-1)>=0;
out(index) = sgn; bitctr = bitctr + 1;
index = index +1; Bits_LIS = Bits_LIS + 1;
else
out(index) = 0; bitctr = bitctr + 1;
index = index +1; Bits_LIS = Bits_LIS + 1;
LIP = [LIP; 2*x-1 2*y-1];%%%将不重要系数坐标加入到LIP中
end
if (bitctr + 1) >= max_bits
if (bitctr < max_bits)
out(length(out))=[];
end
return
end
if abs(m(2*x-1,2*y)) >= 2^n%%%对(2*x-1,2*y)操作
LSP = [LSP; 2*x-1 2*y];
out(index) = 1; bitctr = bitctr + 1;
index = index +1; Bits_LIS = Bits_LIS + 1;
sgn = m(2*x-1,2*y)>=0;
out(index) = sgn; bitctr = bitctr + 1;
index = index +1; Bits_LIS = Bits_LIS + 1;
else
out(index) = 0; bitctr = bitctr + 1;
index = index +1; Bits_LIS = Bits_LIS + 1;
LIP = [LIP; 2*x-1 2*y];
end
if (bitctr + 1) >= max_bits
if (bitctr < max_bits)
out(length(out))=[];
end
return
end
if abs(m(2*x,2*y-1)) >= 2^n%%%对(2*x,2*y-1)操作
LSP = [LSP; 2*x 2*y-1];
out(index) = 1; bitctr = bitctr + 1;
index = index +1; Bits_LIS = Bits_LIS + 1;
sgn = m(2*x,2*y-1)>=0;
out(index) = sgn; bitctr = bitctr + 1;
index = index +1; Bits_LIS = Bits_LIS + 1;
else
out(index) = 0; bitctr = bitctr + 1;
index = index +1; Bits_LIS = Bits_LIS + 1;
LIP = [LIP; 2*x 2*y-1];
end
if (bitctr + 1) >= max_bits
if (bitctr < max_bits)
out(length(out))=[];
end
return
end
if abs(m(2*x,2*y)) >= 2^n%%%对(2*x,2*y)操作
LSP = [LSP; 2*x 2*y];
out(index) = 1; bitctr = bitctr + 1;
index = index +1; Bits_LIS = Bits_LIS + 1;
sgn = m(2*x,2*y)>=0;
out(index) = sgn; bitctr = bitctr + 1;
index = index +1; Bits_LIS = Bits_LIS + 1;
else
out(index) = 0; bitctr = bitctr + 1;
index = index +1; Bits_LIS = Bits_LIS + 1;
LIP = [LIP; 2*x 2*y];
end
if ((2*(2*x)-1) < size(m) & (2*(2*y)-1) < size(m))
LIS = [LIS; LIStemp(i,1) LIStemp(i,2) 1];%%%%类型1表示
LIStemp = [LIStemp; LIStemp(i,1) LIStemp(i,2) 1];
end
LIS(temp,:) = [];%删去一个节点
temp = temp-1;
else %%%%(lis(i,1),lis(i,2))小于2^n时
out(index) = 0; bitctr = bitctr + 1;
index = index +1; Bits_LIS = Bits_LIS + 1;
end
else %%%%%%%%对应93行：LIStemp(i,3) == 0
if bitctr >= max_bits
return
end
max_d = func_MyDescendant(LIStemp(i,1),LIStemp(i,2),LIStemp(i,3),m);
if max_d >= 2^n
out(index) = 1; bitctr = bitctr + 1;
index = index +1;
x = LIStemp(i,1); y = LIStemp(i,2);
LIS = [LIS; 2*x-1 2*y-1 0; 2*x-1 2*y 0; 2*x 2*y-1 0; 2*x 2*y 0];
LIStemp = [LIStemp; 2*x-1 2*y-1 0; 2*x-1 2*y 0; 2*x 2*y-1 0; 2*x 2*y 0];
LIS(temp,:) = []; temp = temp - 1;
else
out(index) = 0; bitctr = bitctr + 1;
index = index +1; Bits_LIS = Bits_LIS + 1;
end
end
i = i+1;
end
% Refinement Pass
temp = 1;
value = floor(abs(2^(n_max-n+1)*m(LSP(temp,1),LSP(temp,2))));
while (value >= 2^(n_max+2) & (temp <= size(LSP,1)))%%%依次输出LSP中各个元素绝对值的第n位有效数字（n为当前的阈值指数）
if bitctr >= max_bits
return
end
s = bitget(value,n_max+2);
out(index) = s; bitctr = bitctr + 1;
index = index +1; Bits_LSP = Bits_LSP + 1;
temp = temp + 1;
if temp <= size(LSP,1)
value = floor(abs(2^(n_max-n+1)*m(LSP(temp,1),LSP(temp,2))));
end
end % Refinement Pass
temp = 1;
n = n - 1;
end
%%for i=1:max_bits-24
% out(i)
% end
%
% fid=fopen(out_stream,'wb');
% fwrite(fid,out(1),'short');
% fwrite(fid,out(2),'uchar');
% fwrite(fid,out(3),'uchar');
%
% for i=1:(max_bits-21)/8
% a=0;
% for j=0:7
% a=a+out(8*i+j+4)*2^j;
% end
% a
% fwrite(fid,a,'uchar');
% end
%
% fclose(fid)