实验六 数学形态及图像压缩
实验要求
对test 目录下的图像进行图像压缩测试，调节参数查看效果。
实验内容
形态学的基本操作(膨胀，腐蚀，开闭等运算)
如何选择形态学结构元素（保持基本形状，去除不相关特征）
JPEG的编码过程
Huffman编码的过程
形态学用到matlab自带的函数：
imdilate ：膨胀
imerode：腐蚀
strel ：生成结构元素
imopen: 开操作
imclose: 闭操作
bwhitmiss: 击中击不中
bwmorph: 实现多种形态学操作
图像压缩编码用到的matlab函数：
huffman **编解码：**mat2huff，huff2mat
JPEG 编解码： im2jpeg，jpeg2im
JPEG2000编解码： im2jpeg2k，jpeg2k2im
参看例子代码： code目录下的CompressionTest.m
Code
im2jpeg2k.m 使用JPEG2000近似压缩一幅图像
im2jpeg2k使用JPEG2000近似值压缩图像
im2jpeg2k Compresses an image using a JPEG 2000 approximation.
Y = im2jpeg2k(X, N, Q) compresses image X using an N-scale JPEG
2K wavelet transform, implicit or explicit coefficient
quantization, and Huffman symbol coding augmented by zero
run-length coding. If quantization vector Q contains two
elements, they are assumed to be implicit quantization
parameters; else, it is assumed to contain explicit subband step
sizes. Y is an encoding structure containing Huffman-encoded
data and additional parameters needed by JPEG2K2IM for decoding.
%im2jpeg2k源程序
function y=im2jpeg2k(x,n,q)
%im2jpeg2k使用JPEG2000近似值压缩图像
%y=im2jpeg2k(x,n,q)使用N尺度jpeg2k小波变换，不明确和明确系数量化和零步长解码去压缩图像
%如果量化向量仅包含两个元素，它们被认为不明确量化参数。
%否则被认为包含明确子带步长。
%y是包含霍夫曼编码数据和加上为编码通过JPEG2K2IM得到的参数的一个编码结构。
global RUNS
error(nargchk(3,3,nargin));%检查输入参数
if nidms(x) ~=2 | ~isreal(x) | ~isnumeric(x) | ~isa(x,'uint8')
error('The input must be a UINT8 image.');
end
if length(q) ~=2 & length(q) ~=3*n+1
error('The quantization step size vector is bad.');
end
%水平平移输入和计算它的小波变量
x=double(x)-128;
[c,s]=wavefast(x,n,'jpeg9.7');
%量化小波系数
q=stepsize(n,q);
sgn=sign(c);sgn(find(sgn==0))=1;c=abs(c);
for k=1:n
qi=3*k-2;
c=wavepaste('h',c,s,k,wavecopy('h',c,s,k)/q(qi));
c=wavepaste('v',c,s,k,wavecopy('v',c,s,k)/q(qi+1));
c=wavepaste('d',c,s,k,wavecopy('d',c,s,k)/q(qi+2));
end
c=wavepaste('a',c,s,k,wavecopy('a',c,s,k)/q(qi+3));
c=floor(c);c=c.*sgn;
%通过创造一个特别的编码为0执行和结束编码来开始和做一个步长表
zrc=min(c(:))-1;eoc=zrc-1;RUNS=[65535];
%找到步长变换指针：'plus'包含零运行开始的指针与相适应的minus是它的结束+1
z=c==0;z=z-[0 z(1:end-1)];
plus=find(z==1);minus=find(z==-1);
%从c中删除零运行
if length(plus) ~=length(minus)
c(plus(end):end)=[];c=[c eoc];
end
%删除从c中所有其他零运行(建立在'plus'和'minus')
for i=length(minus):-1:1
run=minus(i)-plus(i);
if run>10
ovrflo=floor(run/65535);run=run-ovrflo*65535;
c=[c(1:plus(i)-1) repmat([zrc 1],1,ovrflo) zrc ...
runcode(run) c(minus(i):end)];
end
end
%霍夫曼编码和加混合。解码信息
y.runs=uint16(RUNS);
y.s=uint16(s(:));
y.zrc=uint16(-zrc);
y.q=uint16(100*q');
y.n=uint16(n);
y.huffman=mat2huff(c);
所用函数：
wavefast
function [c, s] = wavefast(x, n, varargin)
%wavefast源函数
function [c,s]=wavefast(x,n,varargin)
%wavefast执行一个二维快速小波变换
%[c,l]=wavefast(x,n,lp,hp)用带有最大值和最小值的重构滤波器执行一个FWT图像
%[c,l]=wavefast(x,n,wname)执行同样的程序，但是使用wavefilter为小波wname取得一个带有最大值和最小值的滤波器
%其中尺度参数n必须小于或等于图像维数最大值的log2。滤波器的最大值和最小值必须是偶数。
%为了减少边界失真，x必须是对称扩展的。例如x=[c1 c2 c3 ... cn](一维的)，而它的对称扩展是[... c3 c2 c1 c1 c2 c3 ... cn cn cn-1 cn-2 ...]
%输出：
%矩阵c是系数分解向量：
%c=[a(n) h(n) v(n) d(n) h(n-1) ... v(1) d(1)]
%其中a,h,v和d分别是包含近似值，水平的，垂直的，和二维系数矩阵的列向量。c有3n+1个部分，其中n是小波重构数。
%矩阵s是一个(n+2)*2的簿记矩阵
%s=[sa(n,:);sd(n,:);sd(n-1,:);... sd(1,:);sx],其中sa和sd是近似值和细节大小条目
%检查输入参数
error(nargchk(3,4,nargin));
%nargin用来判断输入变量个数的函数，
%当nargin的值大于4时，nargchk返回字符串'Too many input arguments';当nargin的值小于3时，nargchk返回字符串'Not enough input arguments';
%当nargin的值在3到4之间，nargchk返回空字符串。error以错误的方式显示警告字符串(以红色字体显示)。
if nargin==3
if ischar(varargin{1})
[lp,hp]=wavefilter(varargin{1},'d');
else
error('Missing wavelet name.');
end
else
lp=varargin{1};hp=varargin{2};
end
f1=length(lp);sx=size(x);
if (ndims(lp)~=2) | (min(sx)<2) | ~isreal(x) | ~isnumeric(x)
error('X must be a real,numeric matrix.')
end
if(ndims(lp)~=2) | ~isreal(lp) | ~isnumeric(lp)...
| (ndims(hp)~=2) | ~isreal(hp) | ~isnumeric(hp)...
| (f1~=length(hp)) | rom(f1,2)~=0
error(['LP and HP must be even and equal length real,''numeric filter vectors.']);
end
if ~isreal(n) | ~isnumeric(n) | (n<1) | (n>log2(max(sx)))
error(['N must be a real scalar between 1 and' 'log2(max(size(X))).']);
end
%设置最初的输出数据结构和最初近似值
c=[];s=sx;app=double(x);
%对于每个重构
for i=1:n
%对称地扩展近似值
[app,keep]=symextend(app,f1);
%用hp和采样率计算卷积行。然后用hp和lp得到的卷积列去获取对角线和垂直的系数。
rows=symconv(app,hp,'row',f1,keep);
coefs=symconv(rows,hp,'col',f1,keep);
c=[coefs(:)' c];s=[size(coefs);s];
coefs=symconv(rows,lp,'col',f1,keep);
c=[coefs(:)' c];
% 用hp和采样率计算卷积行。然后用hp和lp得到的卷积列去获取垂直的和下一个近似的系数。
rows=symconv(app,lp,'row',f1,keep);
coefs=symconv(rows,hp,'col',f1,keep);
c=[coefs(:)'c];
app=symconv(rows,lp,'col',f1,keep);
end
%附加最后的近似值结构
c=[app(:)'c];s=[size(app);s];
stepsize
function q=stepsize(n,p)
%stepsize源程序
function q=stepsize(n,p)
%通过分解创造步骤大小顺序的一个子带量化和子带(水平的、垂直的、对角线的和近似值子带的最后分解)
if length(p)==2%不明确量化
q=[];
qn=2^(8-p(2)+n)*(1+p(1)/2^11);
for k=1:n
qk=2^-k*qn;
q=[q (2*qk) (2*qk)(4*qk)];
end
q=[q qk];
else%明确量化
q=p;
end
q=round(q*100)/100;%大约1/100分配
if any(100*q>65535)
error('The quantizing steps are not UNIT16 representable.')
end
if any(q==0)
error('A quantizing step of 0 is not allowed.');
end
wavepaste
function nc=wavepaste(type,c,s,n,x)
%wavepaste源程序
function nc=wavepaste(type,c,s,n,x)
%wavepaste是在小波重构结构是投入的系数
%粘贴x之后返回一个新的重构结构
%输入：
%类型 系数类型
%'a' 近似值系数
%'h' 水平细节
%'v' 垂直细节
%'d' 对角线细节
%[c,s]是一个小波数据结构。
%n是指定的重构水平(type='a'则忽略)。
%x是二维近似值或者细节系数矩阵
error(nargchk(5,5,nargin))
nc=wavework('paste',type,c,s,n,x);
wavecopy
function y=wavecopy(type,c,s,n)
%wavecopy源程序
function y=wavecopy(type,c,s,n)
%wavecopy是在小波重构结构是取得的系数
%决定在type和n，返回一个系数数组
%输入：
%类型 系数类型
%'a' 近似值系数
%'h' 水平细节
%'v' 垂直细节
%'d' 对角线细节
%[c,s]是一个小波数据结构。n是指定的重构水平(type='a'则忽略)。
error(nargchk(3,4,nargin));
if nargin==4
y=wavework('copy',type,c,s,n);
else
runcode
function y=runcode(x)
%runcode源程序
function y=runcode(x)
%在步长表中找到一个零运行，如果没找到，在表中创造一个新的入口。
%返回运行的指针
global RUNS
y=find(RUNS==x);
if length(y) ~=1
RUNS=[RUNS;x];
y=length(RUNS);
end
mat2huff
function =mat2huff(x)
%mat2huff源函数
function =mat2huff(x)
%mat2huff编码一个矩阵
%使用符号概率在最大值和最小值之间建立单位宽度直方图。
%编码数据返回一个结构y
%y:
%y.code x的Huffman编码值，储存在16比特向量。y的其他领域包括额外的解码信息，包括：
%y.min x的最小值加32768
%y.size x的大小
%y.hist x的直方图
%如果x是logical,uint8,uint16,uint32,int8,int16,或者double的整数值，可以直接输入到mat2muff中。
%x的最小值必须是int16.
%如果x是非整值的double。例如：一个图像的值在0和1之间，首先x的大小应该接近整数的范围。
%例如y=mat2huff(255*x)是256灰度级编码
%note:
%Huffman的编码数是round(max(x(:)))-round(min(x(:)))+1.输出x的产生合理的长度.x的行数和列数的最大维数和最小维数是65535.
if ndims(x)~=2 | ~isreal(x) | (~isnumeric) & ~islogical(x))
error('X must be a 2-D real numeric or logical matrix.');
end
%输入x的类型
y.size=uint32(size(x));
%找到x的范围并储存最小值+32768作为一个uint16.
x=round(double(x));
xmin=min(x(:));
xmax=max(x(:));
pmin=double(int16(xmin)); y.min=pmin;
%使用x的最小值和最大值之间之间的单位宽度bin计算x的直方图h，并缩放该直方图，以使其为uint16向量。
x=x(:)';
h=histc(x,xmin,xmax);
if max(h)>65535
h=65535*h/max(h);
end
h=uint16(h);y.hist=h;
%编码输入矩阵和储存结果
map=huffman(double(h));%做Huffman编码图
hx=map(x(:)-xmin+1);%产生单元数组
hx=hx(:)';%转成字符数组
hx(hx=='')=[];%删除空格符
ysize=ceil(length(hx)/16);%计算编码大小
hx16(1:length(hx))=hx;%分配模16向量
hx16=hx16'-'0';%在长度上做hx模16
twos=pow2(15:-1:0);%重建
y.code=uint16(sum(hx16.*twos(ones(ysize,1),:),2))';%转变二维字符串为十进制
jpeg2k2im.m 解码IM2JPEG2K压缩的图像
jpeg2k2im Decodes an IM2JPEG2K compressed image.
X = jpeg2k2im(Y) decodes compressed image Y, reconstructing an
approximation of the original image X. Y is an encoding
structure returned by IM2JPEG2K.
%jpeg2k2im源程序
function x=jpeg2k2im(y)
%解码一个im2jpeg压缩图像
%x=jpeg2k2im(y)解码压缩图像y，重构初始图像x的一个近似值。
%y是通过im2jpeg2k返回编码结构
%也看im2jpeg2k
error(nargchk(1,1,nargin));%检查输入参数
%得到解码参数：尺度、量化向量、步长表大小、零运行编码、数据结束编码、小波簿记数组和步长表
n=double(y.n);
q=double(y.q)/100;
runs=double(y.runs);
rlen=length(runs);
zrc=-double(y.zrc);
eoc=zrc-1;
s=double(y.s);
s=reshape(s,n+2,2);
%计算小波变换大小
c1=prod(s(1,:));
for i=2:n+1
c1=c1+3*prod(s(i,:));
end
%通过零运行编码执行霍夫曼编码
r=huff2mat(y.huffman);
c=[];zi=find(r==zrc);i=1;
for j=1:length(zi)
c=[c r(i:zi(j)-1) zeros(1,runs(r(zi(j)+1)))];
i=zi(j)+2;
end
zi=find(r==eoc);%撤消终止零运行或最后非零运行
if length(zi)==1
c=[c r(i:zi-1)];
c=[c zeros(1,c1-length(c))];
else
c=[c r(i:end)];
end
%非标准化系数
c=c+(c>0)-(c<0);
for k=1:n
qi=3*k-2;
c=wavepaste('h',c,s,k,wavecopy('h',c,s,k)*q(qi));
c=wavepaste('v',c,s,k,wavecopy('v',c,s,k)*q(qi+1));
c=wavepaste('d',c,s,k,wavecopy('d',c,s,k)*q(qi+2));
end
c=wavepaste('a',c,s,k,wavecopy('a',c,s,k)*q(qi+3));
%计算逆小波变换和水平平移
x=waveback(c,s,'jpeg9.7',n);
x=uint8(x+128);
huff2mat
function x=huff2mat(y)
%huff2man源程序
function x=huff2mat(y)
%huff2man解码一个Huffman编码的矩阵
%x=huff2man(y)解码一个Huffman编码的16比特的结构y
%field:
%y.min x的最小值+32768
%y.size x的大小
%y.hist x的直方图
%y.code Huffman编码
%输出x是双精度
if ~isstruct(y) | ~isfield(y,'min') | ~isfield(y,'size') | ...
~isfield(y,'hist') | ~field(y,'code')
error('The input must be a structure as returned by MAT2HUFF');
end
sz=double(y.size);m=sz(1);n=sz(2);
xmin=double(y.min)-32768;%得到x的最小值
map=huffman(double(y.hist));%得到Huffman编码
%为Huffman解码程序创造一个二维调查表
%code包含源符号字符串与连接编码符合，
%当'link'包含映射(+)是给编码符号字符串加上0和1，映射(-)是给解码Huffman编码词放在map中。
%数组'left'是一列编码，目前给link入口加工的。
code=cellstr(char('','0','1'));%设置开始条件
link=[2;0;0];left=[2 3];
fiund=0;tofind=length(map);%处理变量
while length(left) & (found<tofind)
look=find(strcmp(map,code{left(1)}));%map是否是字符串？
if look%是的
link(left(1))=-look;%指出 Huffman图
left=left(2:end);%删除当前编码
found=found+1;%found编码的增量
else
len=length(code);%不是的话，加上2编码与指标
link(left(1))=len+1;%在编码中加指标
link=[link;0;0];%加未加工的编码
code{end+1}=strcat(code{left(1)},'0');
code{end+1}=strcat(code{left(1)},'1');
left=left(2:end);%删除加工编码
left=[left len+1 len+2];%加2个未加工编码
end
end
x=unravel(y.code',link,m*n);%使用C拆散解码
x=x+xmin-1;%调整x最小值的补偿
x=reshape(x,m,n);%是向量成为一个数组
waveback
function [varargout]=waveback(c,s,varargin)
%%wavebac源程序
function [varargout]=waveback(c,s,varargin)
%waveback函数执行一个二维n阶水平部分或完整的小波分解结构[c,s]的重构
%syntax:
%y=waveback(c,s,'wname');输出FWT的逆矩阵y
%y=waveback(c,s,lr,hr);使用低通和高通重构滤波器或名为'wname'的波形滤波器
%[nc,ns]=waveback(c,s,'wname',n);输出新的小波
%[nc,ns]=waveback(c,s,lr,hr,n);n步重构后的分解结构[nc,ns]
%检查输入和输出参数
error(nargchk(3,5,narargin));
error(nargchk(1,2,nargout));
if(ndims(c)~=2) | (size(c,1)~=1)
error('C must be a row vector.');
end
elements=prod(s,2);
if(length(c)<elements(end)) | ~(elements(1)+3*sum(elements(2:end-1))>=elements(end))
error(['[c s] must be a standard wavelet''decomposition structure.']);%标准小波分解结构
end
%在[c,s]中的最大水平
nmax=size(s,1)-2;
%得到第三个输入参数和开始检查标记
wname=varargin{1};filterchk=0;nchk=0;
switch nargin
case 3
if ischar(wname)
[lp,hp]=wavefilter(wname,'r');n=nmax;
else
error('Undefined filter.');
end
if nargout~=1
error('Wrong number of output arguments.');
end
case 4
if ischar(wname)
[lp,hp]=wavefilter(wname,'r');
n=varargin{2};nchk=1;
else
lp=varargin{1};hp=varargin{2};
filterchk=1;n=nmax;
if nargout~=1
error('Wrong number of output arguments.');
end
end
case 5
lp=varargin{1};hp=varargin{2};filterchk=1;
n=varargin{3};nchk=1;
otherwise
error('Improper number of input arguments.');
end
f1=length(lp);
if filterchk
if (ndims(lp)~=2) | ~isreal(lp) | ~isnumeric(lp)...
| (ndims(hp)~=2) | ~isreal(hp) | ~isnumeric(hp)...
|(f1~=length(hp)) | rem(f1,2)~=0
error(['LP and HP must be even and equal length real,''numeric filter vectors.']);
end
end
if nchk & (~isnumeric(n) | isreal(n))
error('N must be a real numeric.');
end
if (n>nmax) & (nargout~=2)
error('Not enough output arguments.');
end
nc=c;ns=s;nnmax=nmax;
for i=1:n
%计算新的近似值
a=symconvup(wavecopy('a',nc,ns),lp,lp,f1,ns(3,:))+...
symconvup(wavecopy('h',nc,ns,nnmax),...
hp,lp,f1,ns(3,:))+...
symconvup(wavecopy('v',nc,ns,nnmax),...
lp,hp,f1,ns(3,:))+...
symconvup(wavecopy('d',nc,ns,nnmax),...
hp,hp,f1,ns(3,:));
%更新分解
nc=nc(4*prod(ns(1,:))+1:end);nc=[a(:)' nc];
ns=ns(3:end,:);
nnmax=size(ns,1)-2;
end
%为完成重构，重定输出格式为二维
if nargout==1
a=nc;nc=repmat(0,ns(1,:));nc(:)=a;
end
varargout{1}=nc;
if nargout==2
varargout{2}=ns;
end
imratio.m 计算压缩比
该函数用于表示两幅图像文件或者变量的比特数的比率。
imratio Computes the ratio of the bytes in two images/variables.
CR = imratio(F1, F2) returns the ratio of the number of bytes in
variables/files F1 and F2. If F1 and F2 are an original and
compressed image, respectively, CR is the compression ratio.
%imratio源函数
function cr=imratio(f1,f2)
%imratio计算两幅图像或两个变量的比特数的比率
%返回f1和f2这两幅图像或两个变量的比特数的比率
%f1是原始图像，f2是压缩后图像，cr是压缩比率
error(nargchk(2,2,nargin));%检查输入参数
cr=bytes(f1)/bytes(f2);%计算比率
%--------------------------------------------------%
function b=bytes(f)
%返回输入f的比特数
%如果f是字符串，则它是一个图像文件名；如果不是，它是一个图像变量
if ischar(f)
info=dir(f);b=info.bytes;
elseif isstruct(f)
%matlab的函数报告每个结构领域的124个字节的内存，因为matlab在内存中储存结构。
%不要计算这额外的内存；相反，将这些内存与每个领域联系
b=0;
fields=fieldnames(f);
for k=1:length(fields)
b=b+bytes(f.(fields{k}));
end
else
info=whos('f');b=info.bytes;
end
compare.m 计算均方误差
%compare源程序
function rmse=compare(f1,f2,scale)
%compare计算和展示两个矩阵的误差
%返回输入f1和f2平方值得平方根误差和展示误差的直方图和程度误差图像。
%当scale省略时，默认为一度程度。
%检查输入参数和设置默认值
error(nargchk(2,3,nargin));
if nargin<3
scale=1;
end
%计算平均值的平方根误差
e=double(f1)-double(f2);
[m,n]=size(e);
rmse=sqrt(sum(e(:).^2)/(m*n));
%如果rmse不等于0，输出误差图像和直方图
if rmse
%形成错误直方图
emax=max(abs(e(:)));
[h,x]=hist(e(:),emax);
if length(h)>=1
figure;bar(x,h,'k');
%对称地显示误差图像
emax=emax/scale;
e=mat2gray(e,[-emax,emax]);
figure,imshow(e);
end
end
基础知识 ——总有一些东西是需知的
数字图像压缩中，三种基本的数据冗余：
1、编码冗余
压缩：平均最优编码长度；霍夫曼编码
2、心理视觉冗余
量化压缩，会导致数据有损压缩。心理视觉冗余是指那些不十分重要的信息。这些冗余在不会削弱图像感知质量的情况下可以消除。
3、像素间冗余
基于预测编码可以消除像素间冗余。
JPEG压缩——JPEG
该算法基于离散余弦变换DCT
JPEG压缩——JPEG2000
基于小波变换
形态学图像处理
膨胀与腐蚀
开操作与闭操作
开操作：先腐蚀再膨胀
闭操作：先膨胀再腐蚀
CompressionTest.m 压缩测试
clc
clear all
f=imread('E:\University\Digital image\实验6\实验六 数学形态学及图像压缩\实验六 数学形态学及图像压缩\test\euro2016.jpg');
figure(1),imshow(f);
f=rgb2gray(f); %转换为灰度图像
figure(2),imshow(f);
c1=im2jpeg2k(f,5,[8 8.5]); % 使用JPEG2000近似压缩一幅图像Imratio
f1=jpeg2k2im(c1); %解码IM2JPEG2K压缩的图像
figure(3),imshow(f1);
cr = imratio(f,c1); %计算两幅图像或变量中的比特率 计算压缩比 原图像/压缩后图像
rms = compare(f,f1); %计算均方误差
f:
f1:
压缩比：15.4382
rms:5.1758
c1=im2jpeg2k(f,8,[8 8.5]); % 使用JPEG2000近似压缩一幅图像Imratio
cr:114.1585
rms:18.4194
1855
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
