% inverse transform from the zigzag format to the matrix form
row=round(sqrt(length(x)));
if row*row~=length(x)
disp('invZigzag() fails!! Must be a square matrix!!');
return;
end;
y=zeros(row,row);
count=1;
for s=1:row
if mod(s,2)==0
for m=s:-1:1
y(m,s+1-m)=x(count);
count=count+1;
end;
else
for m=1:s
y(m,s+1-m)=x(count);
count=count+1;
end;
end;
end;
if mod(row,2)==0
flip=1;
else
flip=0;
end;
for s=row+1:2*row-1
if mod(flip,2)==0
for m=row:-1:s+1-row
y(m,s+1-m)=x(count);
count=count+1;
end;
else
for m=row:-1:s+1-row
y(s+1-m,m)=x(count);
count=count+1;
end;
end;
flip=flip+1;
end;
b = imread('D:\matlab\toolbox\images\imdemos\Lena.bmp');
b = im2uint8(b);
figure(1);imshow(b);% 显示原始图像
b = double(b);
b_temp = b-128;
m = size(b_temp,1);n = size(b_temp,2);
%% 离散余弦变换DCT
b_DCT = zeros(m,n);
% 初始化DCT矩阵,便于计算
N = 8;
DCT_matrix = zeros(N,N);
for k = 1:N
for l = 1:N
DCT_matrix(k,l) = cos((2*(k-1)+1)*(l-1)*pi/(2*N));
end
end
% 将原始图像分块(NXN)并进行DCT变换
for k = 1:m/N
for l = 1:n/N
x1 = 1+N*(k-1);x2 = N*k;y1 = 1+N*(l-1);y2 = N*l;
b_DCT(x1:x2,y1:y2) = my_dct2(b_temp(x1:x2,y1:y2),DCT_matrix,N);
end
end
%% 量化
b_Q = zeros(m,n);
% 初始化量化矩阵(低频部分量化系数小,高频部分量化系数大),滤去更多高频信息
Q_matrix = zeros(N,N);Q11 = 4;step = 3;
for k = 1:N
for l = 1:N
Q_matrix(k,l) = Q11+step*(k-1)+step*(l-1);
end
end
% 将DCT后的矩阵进行量化并去整
for k = 1:m/N
for l = 1:n/N
x1 = 1+N*(k-1);x2 = N*k;y1 = 1+N*(l-1);y2 = N*l;
b_Q(x1:x2,y1:y2) = round(b_DCT(x1:x2,y1:y2)./Q_matrix);% 点除
end
end
%% 编码
% % 构建zigzag变换表(可以快速实现变换),注意仅适用于8X8矩阵
% zigzag = [ 0, 1, 8, 16, 9, 2, 3, 10, ...
% 17, 24, 32, 25, 18, 11, 4, 5, ...
% 12, 19, 26, 33, 40, 48, 41, 34, ...
% 27, 20, 13, 6, 7, 14, 21, 28, ...
% 35, 42, 49, 56, 57, 50, 43, 36, ...
% 29, 22, 15, 23, 30, 37, 44, 51, ...
% 58, 59, 52, 45, 38, 31, 39, 46, ...
% 53, 60, 61, 54, 47, 55, 62, 63];
% zigzag = zigzag+1;
b_encode = zeros(m,n);number = 0;j = 0;
for k = 1:m/N
for l = 1:n/N
x1 = 1+N*(k-1);x2 = N*k;y1 = 1+N*(l-1);y2 = N*l;
% 将8X8矩阵zig-zag变换生成1X64数组,使连续的零个数最多
% B = reshape(b_Q(x1:x2,y1:y2),1,N*N);
% B_encode = B(zigzag);
B_encode = toZigzag(b_Q(x1:x2,y1:y2));
% 游程编码
for i = 1:N*N
if B_encode(i)~=0
% 第一个字节存储下一个非零值与上一个非零值之间0的个数
byte1 = uint8(j);% 使用无符号8位整型
% byte1 = bitshift(uint8(j),4);
% byte1 = byte1+uint8(1);
% 第二个字节存储下一个非零值
byte2 = int8(B_encode(i));% 使用有符号8位整型
b_encode(number+1) = byte1;
b_encode(number+2) = byte2;
j = 0;number = number+2;
else
j = j+1;
end
end
end
end
b_encode = b_encode(1:number);% 压缩后的数据
%% 计算压缩比
code_bits = number*8;% 压缩数据大小(比特)
origin_bits = m*n*8;% 源数据大小(比特)
compression_ratio = origin_bits/code_bits;% 压缩比
disp('压缩比:');disp(compression_ratio);
function [B] = my_dct2(A,DCT_matrix,N)
%dct2主函数,使用于NXN矩阵
if (size(A,1)~=N)||(size(A,2)~=N)
error('矩阵A阶数错误');
end
B = zeros(N,N);
for k = 1:N
for l = 1:N
b_temp = 0;
for m = 1:N
for n = 1:N
b_temp = b_temp+A(m,n)*DCT_matrix(m,k)*DCT_matrix(n,l);
end
end
if k ~= 1
c_u = 1;
else
c_u = 1/sqrt(2);
end
if l ~= 1;
c_v = 1;
else
c_v = 1/sqrt(2);
end
B(k,l) = 1/4*c_u*c_v*b_temp;
end
end
%% 解码
% % 反zigzag变换表
% [temp,izigzag] = sort(zigzag);
% 反游程编码
b_decode = zeros(1,m*n);number = 0;
for i = 1:2:size(b_encode,2)
b_decode(number+1:number+b_encode(i)) = 0;%zeros(1,b_encode(i));
b_decode(number+1+b_encode(i)) = b_encode(i+1);
number = number+b_encode(i)+1;
end
% 反zig-zag变换
b_decode2 = zeros(m,n);
for k = 1:m/N
for l = 1:n/N
x1 = 1+N*(k-1);x2 = N*k;y1 = 1+N*(l-1);y2 = N*l;
start = N*N*((k-1)*n/N+(l-1))+1;
% temp = b_decode(start:start+63);
% b_decode2(x1:x2,y1:y2) = reshape(temp(izigzag),N,N);
b_decode2(x1:x2,y1:y2) = invZigzag(b_decode(start:start+N*N-1));
end
end
%% 量化逆
b_deQ = zeros(m,n);
for k = 1:m/N
for l = 1:n/N
x1 = 1+N*(k-1);x2 = N*k;y1 = 1+N*(l-1);y2 = N*l;
b_deQ(x1:x2,y1:y2) = b_decode2(x1:x2,y1:y2).*Q_matrix;
end
end
%% 逆DCT
b_IDCT = zeros(m,n);
for k = 1:m/N
for l = 1:n/N
x1 = 1+N*(k-1);x2 = N*k;y1 = 1+N*(l-1);y2 = N*l;
b_IDCT(x1:x2,y1:y2) = my_idct2(b_deQ(x1:x2,y1:y2),DCT_matrix,N);
end
end
%% 显示
b_show = uint8(b_IDCT+128);
%b_show = mat2gray(b_show);
figure(2);imshow(b_show);
imwrite(b_show,'Lena2.bmp','bmp');
%% 计算PSNR
psnrvalue = PSNR('Lena.bmp','Lena2.bmp');
disp('峰值信噪比:');disp(psnrvalue);
function [B] = my_dct2(A,DCT_matrix,N)
%dct2主函数,使用于NXN矩阵
if (size(A,1)~=N)||(size(A,2)~=N)
error('矩阵A阶数错误');
end
B = zeros(N,N);
for k = 1:N
for l = 1:N
b_temp = 0;
for m = 1:N
for n = 1:N
b_temp = b_temp+A(m,n)*DCT_matrix(m,k)*DCT_matrix(n,l);
end
end
if k ~= 1
c_u = 1;
else
c_u = 1/sqrt(2);
end
if l ~= 1;
c_v = 1;
else
c_v = 1/sqrt(2);
end
B(k,l) = 1/4*c_u*c_v*b_temp;
end
end
function [B] = my_idct2(A,DCT_matrix,N)
%idct2主函数,使用于NXN矩阵
if (size(A,1)~=N)||(size(A,2)~=N)
error('矩阵A阶数错误');
end
B = zeros(N,N);
for k = 1:N
for l = 1:N
b_temp = 0;
for m = 1:N
for n = 1:N
if m ~= 1
c_u = 1;
else
c_u = 1/sqrt(2);
end
if n ~= 1;
c_v = 1;
else
c_v = 1/sqrt(2);
end
b_temp = b_temp+c_u*c_v*A(m,n)*DCT_matrix(k,m)*DCT_matrix(l,n);
end
end
B(k,l) = 1/4*b_temp;
end
end
function psnrvalue = PSNR(original,unzip)
%% 计算原始图像的信号功率
A = imread(original);
A = double(A);
B = imread(unzip);
B = double(B);
%% 计算MSE
[m,n] = size(A);
[m1,n1] = size(B);
if m~=m1||n~=n1
error('图像大小不一致');
end
msevalue = 0;
for i = 1:m
for j = 1:n
msevalue = msevalue+(A(i,j)-B(i,j))^2;
end
end
msevalue = msevalue/(m*n);
if msevalue == 0
error('图像完全相同');
end
%% 计算峰值信噪比
psnrvalue = 255^2/msevalue;
psnrvalue = 10*log10(psnrvalue);
function [y]=toZigzag(x)
% transform a matrix to the zigzag format
[row col]=size(x);
if row~=col
disp('toZigzag() fails!! Must be a square matrix!!');
return
end
y=zeros(row*col,1);
count=1;
for s=1:row
if mod(s,2)==0
for m=s:-1:1
y(count)=x(m,s+1-m);
count=count+1;
end;
else
for m=1:s
y(count)=x(m,s+1-m);
count=count+1;
end
end
end
if mod(row,2)==0
flip=1;
else
flip=0;
end
for s=row+1:2*row-1
if mod(flip,2)==0
for m=row:-1:s+1-row
y(count)=x(m,s+1-m);
count=count+1;
end
else
for m=row:-1:s+1-row
y(count)=x(s+1-m,m);
count=count+1;
end;
end;
flip=flip+1;
end
2 个解决方案
#1
这么长又不说什么错误
#2
发错位置了，这里是C语言板块啊
#1
这么长又不说什么错误
#2
发错位置了，这里是C语言板块啊