matlab学习笔记

一、序言

1. 图像的输入输出和显示

f=imread("test.png");
f=rgb2gray(f);%rgb图像转化为灰度图像
imshow(f);
imwrite(f,"result.jpg","quality",50);%50代表jpg形式压缩质量0-100

2. matlab支持的四种图像类别

• 灰度级图像(Gray-scale images)
  当灰度级图像是uint8或者uint16类型时，分别具有[0,255]和[0,65535]的整数值
• 二值图像(Binary images)
  矩阵元素为0或者1，使用B=logical(A)将数值数组转化为二值图像
• 索引图像(Indexed images)
• RGB图像(RGB images)

3. 算数运算符

运算符	含义
$+$	数组和矩阵加
$-$	数组矩阵减
$.\ast$	数组元素乘
$\ast$	矩阵乘积
$./$	数组元素右除
$.\setminus$	数组元素左除
$/$	矩阵右除；A/B相当于A*inv(B)
$\setminus$	矩阵左除；A\B相当于inv(A)*B
.^	数组元素幂运算
^	矩阵幂运算
.’	转置
’	共轭转置

4. 语法特性

4.1 数组索引

• matlab与其他语言索引不同，是从1开始！
• 数组元素之间之间不用逗号间隔
• 矩阵每行之间用分号间隔

v = [1 3 5 7 9];%v(1)=1
a = [1 2 3 4;5 6 7 8];%a(1,:)=1 2 3 4

4.2 函数

函数形式：function [output] = name(input)

function result = mysum(a,b)
result = a+b;
end

4.3 函数句柄

• 简单函数句柄

f = @sin; %f(pi)就相当于sin(pi),可以用f()来间接调用sin()

• 匿名函数句柄
  @(input-arg-list) expression

g = @(x) x.^2; %g(2)=2.^2=4
f = @() sin(pi); %f()=sin(pi)=0
timeit(f); %可以用来测试运行时间

tic
fun();% 也可以用来计算函数运行时间
toc

4.4 可变数量输入和输出

• 为检测输入到M函数的参数数目，可以使用nargin，它返回输入参数的个数
• 同样的，可以使用nargout返回输出参数个数
• msg=nargchk(low,high,number)能够检测传递的参数数量是否正确
  例如：

function G=testhv2(x,y,z)
error(nargchk(2,3,nargin); %当只有一个参数:testv2(6)会报错！

function g=intrans(f,method,varargin)
···
if length(varargin) == 2 %从第三个参数开始数，如果有2个参数
	...
end

4.5 一些常用函数

g=gscale(f,method,low,high);%先将f归一化（0,1）,再转化为uint8（默认）类型
method:
'full8'（默认）,'full16'，这两种参数，low和high会被忽略；
'minmax'，灰度会被映射到(low,high)

fplot(fhandle,limits,'LineSpec');%绘图函数，fhandle是要绘制函数的句柄

4.6 快捷键

ctrl+r 快速注释
ctrl+t 取消注释

二、灰度变换与空间滤波

1. 灰度变换

1.1 利用imadjust()函数变换

g = imadjust(f,[low_in high_in],[low_out high_out],gamma)

low_high=strechlim(f); %返回图像灰度值的最小值和最大值[low high]
low_high=strechlim(f，[low high]); %指定以low和high充满的部分

将[low_in high_in]的灰度值映射为[low_out high_out]，gamma控制曲线形状，默认为1，即线性映射。通常和strechlim(f)函数一起用。

f=imread("test.png");
f=mat2gray(f);
imshow(f);
g=imadjust(f,stretchlim(f,[0.5 1]),[0 1]); %将f中[0.5 1]范围的像素映射为[0 1]上
figure,imshow(g);

1.2 利用对数变换

$g=c\ast log(1+f)$

代码实现：

g=im2uint8(mat2gray(log(1+double(f) )));

1.3 指定任意灰度变换

T=[0 1];
Z=linspace(1,0,numel(T)); %创建[1,0]之间长度为numel(T)的数组
g=interp1(Z,T,g);

2. 直方图

2.1 直方图表示

$$h(r_k)=n_k$$
其中$r_k$是第k级灰度，$n_k$ 为出现$r_k$这种灰度级的像素数。matlab画直方图的函数为：

imhist(f,b); %b为像素级数目，默认为256

2.2 图像归一化

$$p(r_k)=h(r_k)/n$$
$n$为像素总数，即numel(f)。

f=imread("test.png");
f=rgb2gray(f);
[h,x]=imhist(f,25);%每25个灰度值分为一组
h=h/numel(f);
bar(horz,h);%画出条形直方图

2.3 利用bar画条形直方图

f=imread("test.png");
f=rgb2gray(f);
h=imhist(f,25);%每25个灰度值分为一组
horz=linspace(0,255,25);%在0 255之间插25个值
bar(horz,h);%画出条形直方图
axis([0 255 0 60000]);%设置横纵坐标最大最小值
%在坐标添加标记
set(gca,"xtick",0:50:255); %gca表示获取当前轴
set(gca,"ytick",0:20000:60000);

2.4 利用杆状图画直方图

f=imread("test.png");
f=rgb2gray(f);
h=imhist(f,25);%每25个灰度值分为一组
horz=linspace(0,255,25);%在0 255之间插25个值
stem(horz,h,'fill');%画出条形直方图
axis([0 255 0 60000]);%设置横纵坐标最大最小值
%在坐标添加标记
set(gca,"xtick",0:50:255); %gca表示获取当前轴
set(gca,"ytick",0:20000:60000);

2.5 利用曲线画直方图

f=imread("test.png");
f=rgb2gray(f);
h=imhist(f);
plot(h,"r--o");%画出折线，r表示红色，--表示虚线，o表示点
axis([0 255 0 4000]);%设置横纵坐标最大最小值

2.6 直方图均衡化

$$S_k=T(r_k)=\sum _{j=0}^k{p}_r(r_j)=\sum _{j=0}^k{n}_j/n$$
$r_k$表示第$k$级灰度，$p_r$表示r灰度的概率，$T(r_k)$即归一化直方图的累加求和

均衡化过程：
假设有如下图像：

得图像的统计信息如下图所示，并根据统计信息完成灰度值映射：

映射后的图像如下所示：

以8个灰度的图像为例说明：

f=imread("result.jpg");
figure,imshow(f),title("原图");
figure,imhist(f),title("原图直方图");

g=histeq(f); %均衡化函数

figure,imshow(g),title("均衡化图");
figure,imhist(g),title("均衡化直方图");

结果：

2.7 直方图规定化

含义：给定特定的直方图，将图像按照给定直方图均衡化输出。

工具函数：g = histeq(f,hspec);
hspec为直方图（一个规定值的行向量），g为均衡化输出图像

代码示例：

function p = twomodegauss(m1, sig1, m2, sig2, A1, A2, k)
%  产生一个双峰的高斯函数
%  P = TWOMODEGAUSS(M1, SIG1, M2, SIG2, A1, A2, K)在区间[0, 1]上产生一个
%  双峰类高斯的函数。P是标准化的包含256个元素的向量，从而SUM(P)=1.双峰各自的
%  均值和标准差是(M1, SIG1) 和 (M2, SIG2).A1和A2是双峰的峰值。当输出为标准化时
%  A1和A2的相对值比较重要。K是基线，一个比较好的取值集是M1 = 0.15,
%  SIG1 = 0.05, M2 = 0.75, SIG2 = 0.05, A1 = 1, A2 =0.07和 K = 0.002.
 
c1 = A1 * (1 / ((2 * pi) ^ 0.5) * sig1);
k1 = 2 * (sig1 ^ 2);
c2 = A2 * (1 / ((2 * pi) ^ 0.5) * sig2);
k2 = 2 * (sig2 ^ 2);
z = linspace(0, 1, 256);

p = k+c1*exp(-((z-m1).^2)./k1)+c2*exp(-((z-m2).^2)./k2);
p = p./sum(p(:));

p=twomodegauss(0.15,0.05, 0.75, 0.05,1, 0.07, 0.002);
plot(p);

%使用定义的p直方图均衡化图像f
g=histeq(f,p)

2.8 自适应直方图均衡

用直方图匹配方法逐个处理图像中的较小区域（小片），然后使用双线性内插法将相邻的小片组合起来。

工具函数：g = adapthisteq(f,param1,val1,param2,val2,...);

代码示例：

% 参数说明此处略
g1 = adapthisteq(f);
g2 = adapthisteq(f,'NumTiles',[25 25]); %使用双线性内插将相邻25行25列的小片组合起来
g3 = adapthisteq(f,'NumTiles',[25 25],'ClipLimit',0.05);

3. 空间滤波

3.1 线性滤波

滤波函数：g = imfilter(f,w,filter_mode,boundary_options,size_options)

• filter_mode有两种：corr和conv两种，rot90(w,2)则两种模式效果相同
• boundary有：replicate、symmetric、circular
• size_options有：full、same

产生模板w函数：fspecial(type,params)

• H =fspecial('average',HSIZE)
• H=fspecial('gaussian',[r c],sig)
• H=fspecial('laplacian',aplha)

3.2 非线性滤波

工具函数：nlfilter不常用，速度慢和g = colfilt(f,[m n],'sliding',fun);

• sliding表示处理过程是m*n区域输入图像f中逐像素滑动
• fun是一个句柄，函数fun必须对A的每一列进行操作，并返回一个行向量v，行向量的第k个值是对A中的第k列进行fun操作后的结果。

在使用colfilt()时，滤波前必须显示地填充输入图像，为此需要使用padarray:

fp=padarray(f,[r c],method,direction)

• [r c]表示要充填f的行数和列数
• method：replicate、symmetric、circular
• direction：pre、post、both

代码示例：

% prod(A,1)函数表示列元素相乘形成一个新的行向量
gmean = @(A) prod(A,1)^1/size(A,1);
r = padarray(f,[m n],'replicate');
g = colfilt(f,[m n],'sliding',@gmean);

%最后删除充填的部分
[M,N]=size(f);
g=g((1:M)+m,(1:N)+n);

非线性滤波器：g=ordfilt2(f,order,domin)

• order表示用第几个数去代替该元素g=ordfilt2(f,(m*n+1)/2,ones(m,n)表示中值滤波

中值滤波函数：g=medfilt2(f,[m n],padopt)

• padopt可选：zeros(默认)、symmetric、indexed（f是double类用1填充，否则用0）

三、频率域滤波

四、图像配准函数

图像配准的过程：
见官方文档

fitgeotrans函数，用来计算浮动点与固定点之间的映射矩阵，无法设置参数，感觉不如estimateGeometricTransform函数准确，函数具体形式如下：

• tform = fitgeotrans(movingPoints,fixedPoints,transformationType)
• tform = fitgeotrans(movingPoints,fixedPoints,'polynomial',degree)
• tform = fitgeotrans(movingPoints,fixedPoints,'pwl')
• tform = fitgeotrans(movingPoints,fixedPoints,'lwm',n)

estimateGeometricTransform函数，同样用来评估浮动点与固定点之间的映射矩阵，可以通过参数设定最大实验次数，映射最大差值，输出满足条件的对应映射点和映射矩阵，函数具体形式如下：

• tform = estimateGeometricTransform(matchedPoints1,matchedPoints2,transformType)
• [tform,inlierpoints1,inlierpoints2] = estimateGeometricTransform(matchedPoints1,matchedPoints2,transformType)
• [___,status] = estimateGeometricTransform(matchedPoints1,matchedPoints2,transformType)

imwarp函数，对图像应用变换矩阵，函数具体形式如下：

• B = imwarp(A,tform)

transformPointsForward函数，点应用变换矩阵，即矩阵左乘运算，函数具体形式如下：

• [x,y] = transformPointsForward(tform,u,v)

imregister函数，将浮动图像与固定图像进行自适应性配准，可以设置优化器、迭代次数、检测半径，达到更好的效果，是一种非刚性图像配准方法，函数具体形式如下：

• moving_reg = imregister(moving,fixed,transformType,optimizer,metric)

imregtform函数，计算浮动图像与固定图像之间的变换矩阵，该函数底层与imregister相同，（由于隐藏了检测点和映射点的过程，不敢确定特别准确，精确度要求不高的情况用该方法），函数具体形式如下：

• tform = imregtform(moving,fixed,transformType,optimizer,metric)

其他相关函数(详情查看matlab文档)：
maketform
tformfwd
imtransform
tformfwd
tforminv
注：monomodal单模、multimodal多模图像的区别： 单模图像指具有相似的亮度和对比度的两幅图像，来源于同一类型的扫描仪或传感器；多模图像具有不同的亮度和对比度。这些图像可以来自两种不同类型的设备，也可以来自单个设备，使用不同曝光或者不同的波长拍摄。