第十三章：数字图像处理

最新推荐文章于 2022-12-07 09:03:40 发布

victor_cs_

最新推荐文章于 2022-12-07 09:03:40 发布

阅读量1k

点赞数

分类专栏：数学建模算法与应用文章标签：数字图像数学建模

本文链接：https://blog.csdn.net/qq_41000485/article/details/99705979

版权

数学建模算法与应用专栏收录该内容

16 篇文章 25 订阅

订阅专栏

13.1 数字图像概述

13.1.1 图像的基本概念

连续图像：在二维坐标系中连续变化的图像，即图像的像点是无限稠密的，同时具有灰度值
离散图像：用一个数字序列表示的图像

13.1.2 图像的数字化采样

按照图像空间的坐标测量该位置上像素的灰度值。
一幅数字图像在Matlab中可以很自然的表示成矩阵，矩阵中的元素称为像素。每一个像素都有x和y两个坐标，表示其在图像中的位置。另外还有一个值，称灰度值，对应于原始模拟图像在该点处的亮度。

13.1.3 数据类

数值数据类：doubel, uint8, uint16, uint32, int8, int16, int32, single
字符类：char
逻辑数据类：logical

13.1.4 图像类型

（1）二值图像
通常用于文字、线条图的扫描识别和掩膜图像的存储。在Matlab中是一个取值只有0和1的逻辑数组。使用logical函数可以把数值数组转换为二值数组。
（2）灰度图像
矩阵元素的整数取值范围通常为[0, 255]。
若是uint16类，则范围为[0, 65535]。
若是double类，则范围为[0, 1]。
（3）RGB彩色图像
是彩色像素的一个mn3的数组，其中每一个彩色像素点都是在特定空间位置的彩色图像相应的红、绿、蓝三个分量。
（4）索引图像
有两个分量，即数据矩阵X和彩色映射矩阵map。矩阵map是一个大小为m*3且有范围在[0, 1]之间的浮点值构成的double类数组。

13.1.5 数据类和图像类型的转换

用Matlab工具箱，转换函数有ind2gray, gray2ind, rgb2ind, ind2rgb, ntsc2rgb, rgb2ntsc等，可以用imtool命令。

13.2 亮度变换和空间滤波

Matlab中函数imadjust是对图像进行亮度变换的工具。
图像翻转：

f=imread('tu1.bmp');  %读原图像
g=imadjust(f,[0; 1],[1; 0]); %进行图像翻转
subplot(1,2,1), imshow(f) %显示原图像
subplot(1,2,2), imshow(g) %显示翻转图像

13.2.1 线性空间滤波器

使用拉普拉斯滤波器增强图像，使用Matlab函数fspecial

13.2.2 图像恢复实例

利用Matlab工具箱

f=imread('tu2.bmp'); 
h1=fspecial('laplacian',0); 
g1=f-imfilter(f,h1); %中心为-4，c=-1，即从原图像中减去拉普拉斯算子处理的结果
h2=[1 1 1; 1 -8 1; 1 1 1]; 
g2=f-imfilter(f,h2); %中心为-8
subplot(1,3,1),imshow(f) 
subplot(1,3,2),imshow(g1) 
subplot(1,3,3),imshow(g2)

13.2.3 非线性空间滤波器

利用函数ordfilt2，生成统计排序滤波器。

f=imread('Lena.bmp');
f1=imnoise(f,'salt & pepper',0.02); %加椒盐噪声
g=medfilt2(f1); %进行中值滤波
subplot(1,3,1),imshow(f),title('原图像')
subplot(1,3,2),imshow(f1),title('被椒盐噪声污染的图像')
subplot(1,3,3),imshow(g),title('中值滤波图像')

13.3 频域变换

13.3.1 傅里叶变换

（1）二维连续傅里叶变换
（2）二维离散傅里叶变换（DFT）
（3）基于离散傅里叶变换的频域滤波

计算 $1)^{x+y}$ 乘以输入图像进行中心变换
计算图像 $f_c$ (x, y)的离散傅里叶变换
用滤波器H(u, v)作用F(u, v)
计算G(u, v)的离散傅里叶逆变换，并取实部
用 $1)^{x+y}$ 乘以 $g_p$ (x, y)得到中心还原滤波图像

clc, clear
cm=imread('cameraman.tif'); 
[n,m]=size(cm); %计算图像的维数
cf=fft2(cm); %进行傅里叶变换
cf=fftshift(cf); %进行中心变换
u=[-floor(m/2):floor((m-1)/2)] %水平频率
v=[-floor(n/2):floor((n-1)/2)] %垂直频率
[uu,vv]=meshgrid(u,v); %频域平面上的网格节点
bl=1./(1+(sqrt(uu.^2+vv.^2)/15).^2); %构造一阶巴特沃兹低通滤波器
cfl=bl.*cf; %逐点相乘，进行低通滤波
cml=real(ifft2(cfl)); %进行傅里叶逆变换，并取实部
%cml=ifftshift(cml);
cml=uint8(cml); %必须进行数据格式转换
subplot(1,2,1), imshow(cm)  %显示原图像
subplot(1,2,2), imshow(cml) %显示滤波后的图像

13.3.2 离散余弦变换(DCT)

实现：
（1）快速傅里叶变换（FFT），利用dct2函数实现
（2）DCT变换矩阵，利用dctmtx函数实现
画图的程序：

clc, clear, T=dctmtx(8); %8*8的DCT变换矩阵
colormap('gray'); %设置颜色映射矩阵
for m = 1:8
for n = 1:8
subplot(8,8,(m-1)*8+n);
Y=zeros(8); Y(m,n)=1; %8*8矩阵中只有一个元素为1，其余元素都为0
X = T'*Y*T; %作逆DCT变换
imagesc(X); 
axis square %画图区域是方形
axis off %不显示轴线和标号
end
end

例一：

I = imread('cameraman.tif'); 
I = im2double(I); %数据转换成double类型
T = dctmtx(8); %T为8*8的DCT变换矩阵
dct = @(block_struct) T * block_struct.data * T'; %定义正DCT变换的匿名函数，这里的block_struct是Matlab内置的结构变量
B = blockproc(I,[8 8],dct); %作正DCT变换
mask = [1 1 1 1 0 0 0 0
1 1 1 0 0 0 0 0
1 1 0 0 0 0 0 0
1 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0]; %给出掩膜矩阵
B2 = blockproc(B,[8 8],@(block_struct) mask .* block_struct.data); %提取低频系数
invdct = @(block_struct) T' * block_struct.data * T; 
I2 = blockproc(B2,[8 8],invdct); %作逆变换
subplot(1,2,1), imshow(I) 
subplot(1,2,2), imshow(I2)

例二：对RGB彩色图像做压缩

clc, clear
f0=imread('tu3.bmp'); 
f1=double(f0); 
for k=1:3
g(:,:,k)=dct2(f1(:,:,k));  %对RGB各个分量分别作离散余弦变换
end
g(abs(g)<0.1)=0;  %°把DCT系数小于0.1的变成0
for k=1:3
f2(:,:,k)=idct2(g(:,:,k)); %作DCT逆变换
end
f2=uint8(f2); %把数据转换成uint8格式
imwrite(f2,'tu4.bmp'); 
subplot(1,2,1),imshow(f0)
subplot(1,2,2),imshow(f2)

13.3.3 图像保真度和质量

保真度（逼真度）准则：对信息损失的测度以描述解码图像相对于原始图像的偏离程度
（1）客观保真度准则
（2）主观保真度准则

13.4 数字图像的水印防伪

13.4.1 基于矩阵奇异值分解

（1）矩阵的奇异值分解（SVD）与图像矩阵的能量

f=imread('Lena.bmp');
f=double(f);   %uint8类型数据无法作奇异值分解，必须转换成double类型
[u,s,v]=svd(f); %进行奇异值分解，s为对角矩阵
s=diag(s); %提出对角矩阵的对角线元素，得到一个向量
smax=max(s), smin=min(s) %求最大奇异值和最小奇异值
s1=s; s1(21:end)=0; %只保留前20个大的奇异值，其他奇异值置0
s1=diag(s1);  %吧向量变成对角矩阵
g=u*s1*v';    %计算压缩以后的图像矩阵
g=uint8(g);   %必须转换成原数据类
imwrite(g,'Lena2.bmp') 
subplot(1,2,1), imshow('Lena.bmp') 
subplot(1,2,2), imshow(g) 
figure, plot(s,'.','Color','k') %画出奇异值对应的点

（2）水印嵌入
需要选择合适的嵌入强度因子，小的有利于水印的透明性，大的有利于增强算法的鲁棒性
（3）水印提取

clc, clear
A=imread('tu5.bmp'); 
W=imread('tu6.bmp'); 
[m1,m2,m3]=size(W);  
A0=A([1:m1],[1:m2],:); 
A0=double(A0); W=double(W); 
a=0.05; %嵌入强度因子为0.05
for i=1:3
    [U1{i},S1{i},V1{i}]=svd(A0(:,:,i)); 
    A1(:,:,i)=S1{i}+a*W(:,:,i); 
    [U2{i},S2{i},V2{i}]=svd(A1(:,:,i)); 
    A2(:,:,i)=U1{i}*S2{i}*V1{i}'; 
end
AW=A;  %整体水印合成图片初始化
AW([1:m1],[1:m2],:)=A2; 
AW=uint8(AW); W=uint8(W); 
subplot(1,3,1), imshow(A) 
subplot(1,3,2), imshow(W) 
subplot(1,3,3), imshow(AW) 
%以下是水印的提取
AWstar=imnoise(AW,'gaussian',0,0.01); %加入高斯噪声
A2star=AWstar([1:m1],[1:m2],:);
A2star=double(A2star); 
for i=1:3
    [U3{i},S2star{i},V3{i}]=svd(A2star(:,:,i)); 
    A1star(:,:,i)=U2{i}*S2star{i}*V2{i}';  
    Wstar(:,:,i)=(A1star(:,:,i)-S1{i})/a;  
end
for i=1:3
Wstar(:,:,i)=medfilt2(Wstar(:,:,i));  
end
Wstar=uint8(Wstar); 
figure, subplot(1,2,1), imshow(AWstar) 
subplot(1,2,2), imshow(Wstar)

13.4.2 基于DCT变换

（1）水印嵌入
（2）水印提取

clc, clear 
a=imread('Lena.bmp'); 
[M1,N1]=size(a); 
a=im2double(a);
knum1=M1/8; knum2=N1/8; %把载体图像划分成8*8的子块，高和长方向划分的块数
b0=imread('tu7.bmp'); 
b=im2double(b0); 
subplot(1,2,1), imshow(a)  
subplot(1,2,2), imshow(b)  
mask1=[1 1 1 1 0 0 0 0
       1 1 1 0 0 0 0 0
       1 1 0 0 0 0 0 0
       1 0 0 0 0 0 0 0
       0 0 0 0 0 0 0 0
       0 0 0 0 0 0 0 0
       0 0 0 0 0 0 0 0
       0 0 0 0 0 0 0 0]; %给出低频的掩膜矩阵
ind1=find(mask1==1); %低频系数的位置
mask2=[0 0 0 0 1 1 0 0
       0 0 0 1 1 0 0 0
       0 0 1 1 0 0 0 0
       0 1 1 0 0 0 0 0
       1 1 0 0 0 0 0 0
       1 0 0 0 0 0 0 0
       0 0 0 0 0 0 0 0
       0 0 0 0 0 0 0 0]; %给出中频的掩膜矩阵
ind2=find(mask2==1); %中频系数的位置
[M2,N2]=size(b); 
L=M2*N2; 
knum3=ceil(M2*N2/11); 
b=b(:); b(L+1:11*knum3)=0; 
T=dctmtx(8); 
ab=zeros(M1,N1); 
k=0; 
for i=0:knum1-1  %该两层循环进行水印嵌入
    for j=0:knum2-1
        xa=a([8*i+1:8*i+8],[8*j+1:8*j+8]);  
        ya=T*xa*T';  
        coef1=(mask1+mask2).*ya; 
        if k<knum3
        coef1(ind2)=coef1(ind2)+0.05*b(11*k+1:11*k+11); 
        end
        ab([8*i+1:8*i+8],[8*j+1:8*j+8])=T'*coef1*T;
        k=k+1;
    end
end
acha=ab-a;
k=0; tb=zeros(11*knum3,1); 
for i=0:knum1-1
    for j=0:knum2-1
        xa2=acha([8*i+1:8*i+8],[8*j+1:8*j+8]);  
        ya2=T*xa2*T'; 
        coef2=mask2.*ya2;
        if k<knum3
        tb(11*k+1:11*k+11)=20*coef2(ind2); 
        end
        k=k+1;
    end
end
tb(L+1:end)=[]; 
tb=reshape(tb,[M2,N2]); 
figure, subplot(1,2,1), imshow(ab) 
subplot(1,2,2), imshow(tb)

13.5 图像的加密和隐藏

clc, clear
a=imread('tu8.bmp'); ws1=size(a);    
b=imread('tu9.bmp'); ws2=size(b);   
nb=imresize(b,ws1(1:2));  %把载体图像变换成与保密图像同样大小
key=-0.400001;  %给出密匙，即混沌序列的初始值
L=max(ws1); x(1)=key; y(1)=key; alpha=1.4; beta=0.3;
for i=1:L-1     %生成两个混沌序列
     x(i+1)=1-alpha*x(i)^2+y(i); y(i+1)=beta*x(i);
end
x(ws1(1,1)+1:end)=[]; 
[sx,ind1]=sort(x); [sy,ind2]=sort(y); %对混沌序列从小到大排序
ea(ind1,ind2,:)=a; %打乱保密图像的行序和列序，生成加密图像矩阵
imshow(ea) 
nb2=bitand(nb,240); 
ea2=bitand(ea,240); 
ea2=ea2/16; %加密图像高四位移到低4位
da=bitor(nb2,ea2); %把加密图像嵌入载体图像的低4位，构造合成图像
da2=bitand(da,15)*16; 
da3=da2(ind1,ind2,:); %对加密图像进行解密
figure, subplot(1,2,1), imshow(da3) 
subplot(1,2,2),imshow(da)

victor_cs_

关注

0
点赞
踩
9

收藏

觉得还不错? 一键收藏
0
评论
第十三章：数字图像处理

13.1 数字图像概述13.1.1 图像的基本概念连续图像：在二维坐标系中连续变化的图像，即图像的像点是无限稠密的，同时具有灰度值离散图像：用一个数字序列表示的图像13.1.2 图像的数字化采样按照图像空间的坐标测量该位置上像素的灰度值。一幅数字图像在Matlab中可以很自然的表示成矩阵，矩阵中的元素称为像素。每一个像素都有x和y两个坐标，表示其在图像中的位置。另外还有一个值，称灰度值...
复制链接

扫一扫