数字图像处理-第十章和第十一章
学习了第10章的后半部分和第11章
plot
x=0:pi/100:2*pi;
y=2exp(-0.5x).sin(2pi*x);
plot(x,y)
程序执行后，打开一个图形窗口，在其中绘制出如下曲线
注意：指数函数和正弦函数之间要用点乘运算，因为二者是向量。
例52 绘制曲线
这是以参数形式给出的曲线方程，只要给定参数向量，再分别求出x,y向量即可输出曲线：
t=-pi:pi/100:pi;
x=t.cos(3t);
y=t.sin(t).sin(t);
plot(x,y)
程序执行后，打开一个图形窗口，在其中绘制出如下曲线
以上提到plot函数的自变量x,y为长度相同的向量，这是最常见、最基本的用法。实际应用中还有一些变化。分别说明：
①
2． 含多个输入参数的plot函数
plot函数可以包含若干组向量对，每一组可以绘制出一条曲线。含多个输入参数的plot函数调用格式为：plot(x1，y1，x2，y2，…，xn，yn)
如下列命令可以在同一坐标中画出3条曲线。
x=linspace(0,2*pi,100);
plot(x,sin(x),x,2sin(x),x,3sin(x))
当输入参数有矩阵形式时，配对的x,y按对应的列元素为横坐标和纵坐标绘制曲线，曲线条数等于矩阵的列数。
x=linspace(0,2*pi,100);
y1=sin(x);
y2=2*sin(x);
y3=3*sin(x);
x=[x;x;x]';
y=[y1;y2;y3]';
plot(x,y,x,cos(x))
x,y都是含有三列的矩阵，它们组成输入参数对，绘制三条曲线；x和cos(x)又组成一对，绘制一条余弦曲线。
利用plot函数可以直接将矩阵的数据绘制在图形窗体中，此时plot函数将矩阵的每一列数据作为一条曲线绘制在窗体中。如
A=pascal(5)
A =
1 1 1 1 1
1 2 3 4 5
1 3 6 10 15
1 4 10 20 35
1 5 15 35 70
plot(A)
3． 含选项的plot函数
Matlab提供了一些绘图选项，用于确定所绘曲线的线型、颜色和数据点标记符号。这些选项如表所示：
线型
颜色
标记符号
- 实线
b蓝色
. 点
- s 方块
- 虚线
g绿色
o 圆圈
d 菱形
-. 点划线
r红色
× 叉号
∨朝下三角符号
-- 双划线
c青色
- 加号
∧朝上三角符号
m品红
- 星号
<朝左三角符号
y黄色
朝右三角符号
k黑色
p 五角星
w白色
h 六角星
例 用不同的线型和颜色在同一坐标内绘制曲线 及其包络线。
x=(0:pi/100:2*pi)';
y1=2exp(-0.5x)*[1,-1];
y2=2exp(-0.5x).sin(2pi*x);
x1=(0:12)/2;
y3=2exp(-0.5x1).sin(2pi*x1);
plot(x,y1,'k:',x,y2,'b--',x1,y3,'rp');
在该plot函数中包含了3组绘图参数，第一组用黑色虚线画出两条包络线，第二组用蓝色双划线画出曲线y，第三组用红色五角星离散标出数据点。
集合运算
f = imread('Fig0903(utk).tif');
g = imread('Fig0903(gt).tif');
subplot(2, 3, 1), imshow(f), title('UTK');
subplot(2, 3, 2), imshow(g), title('GT');
subplot(2, 3, 3), imshow(~f), title('UTK的补集');
subplot(2, 3, 4), imshow(f | g), title('UTK和GT的合集');
subplot(2, 3, 5), imshow(f & g), title('UTK和GT的交集');
subplot(2, 3, 6), imshow(f & ~g), title('UTK和GT的差集');
霍夫变换的线检测
S=sparse(r, c, s, m, n)将矩阵变为稀疏矩阵
full(S) 可以将稀疏矩阵还原
函数 intline(x1, x2, y1, y2) 将二个点连起来，并返回这条线的向量
[x, y] = intline(0, 100, 0, 100);
imshow(bound2im(cat(2, x, y)));
函数 [H,theta,rho] = hough(f)
f = zeros(101, 101);
f(1, 1) = 1;
f(101, 1) = 1;
f(1, 101) = 1;
f(101, 101) = 1;
f(51, 51) = 1;
[H,theta,rho] = hough(f);
subplot(1, 2, 1), imshow(f), title('原图');
subplot(1, 2, 2), imshow(H,[],'XData',theta,'YData',rho,...
'InitialMagnification','fit'), title('计算后');
xlabel('\theta'), ylabel('\rho');
axis on, axis normal;
函数[r,c,hnew] = houghpeaks(h,numpeaks,threshold,nhood)
hough峰值检测
1.找到包含有最大值的hough变换单元并记下它的位置
2.把第一步中找到的最大值点的邻域中的hough变换单元设为零
3.重复该步骤，直到找到需要的峰值数时为止，或者达到一个指定的阈值时为止
函数[r,c] = houghpixels(f,theta,rho,rbin,cbin)
找到影响峰值相关位置
lines = houghlines(f,theta,rho,rr,cc,fillgap,minlength)
相关位置像素合成线段
hough变换的应用
直线检测
一条直线在图像中是一系列离散点的集合，通过一个直线的离散极坐标公式，
可以表达出直线的离散点几何等式如下：X *cos(theta) + y * sin(theta) = r 其中角度theta指r与X轴之间的夹角，r为到直线几何垂直距离。如果我们能绘制每个(r, theta)值根据像素点坐标P(x, y)值的话，那么就从图像笛卡尔坐标系统转换到极坐标霍夫空间系统，这种从点到曲线的变换称为直线的霍夫变换.
close all;
f = imread('Fig1006.tif');
f = edge(f, 'canny', [0.04, 0.10], 1.5 );
[H,theta,rho] = hough(f,'ThetaRes',0.2);
imshow(H,[], 'XData', theta, 'YData', rho, 'InitialMagnification', 'fit');
axis on,axis normal;
xlabel('\theta'),ylabel('\rho');
peaks=houghpeaks(H,5);
hold on;
plot(theta(peaks(:,2)), rho(peaks(:,1)),...
'LineStyle', 'none', 'marker', 's', 'color', 'w');
lines=houghlines(f,theta,rho,peaks);
figure;imshow(f), hold on;
for k = 1:length(lines)
xy = [lines(k).point1 ; lines(k).point2];
plot(xy(:, 1),xy(:,2),'LineWidth', 4, 'Color', [.8 .8 .8]);
圆检测
(x-a)2+(y-b)2=r^2;
x=a+rcos(theta),y=b+rsin(theta)
(x和y为实际的图像空间某个边缘点的坐标，a和b为其对应的参数空间的坐标)
中心点处的坐标值必定最强.
function [hough_space,hough_circle,para] = hough_circle(BW,step_r,step_angle,r_min,r_max,thresh);
% hough_space:参数空间，h(a,b,r)表示圆心在(a,b)半径为r的圆上的点数
% hough_circle:二值图像，检测到的圆
% para:检测到的圆的圆心、半径
[m,n] = size(BW);
size_r = round((r_max-r_min)/step_r)+1;
size_angle = round(2*pi/step_angle);
hough_space = zeros(m,n,size_r); %对应于参数平面，将其所有数据置为0.
[rows,cols] = find(BW);
ecount = size(rows); %非零坐标的个数
% Hough变换
% 将图像空间(x,y)对应到参数空间(a,b,r)
% a = x-r*cos(angle)
% b = y-r*sin(angle)
%hough_space(a,b,r)确立了一个三维的立体坐标系，每个坐标点存放着一个数值，此数值代表着圆心为（a，b）半径为r的圆的可信度
for i=1:ecount
for r=1:size_r %单个点在所有半径空间内检测
for k=1:size_angle
a = round(rows(i)-(r_min+(r-1)*step_r)*cos(k*step_angle));
b = round(cols(i)-(r_min+(r-1)*step_r)*sin(k*step_angle));
if(a>0&&a<=m&b>0&b<=n)
hough_space(a,b,r) = hough_space(a,b,r)+1; %让圆对应的参数平面的点的值都加１。
end
end
end
end
% 搜索超过阈值的聚集点
max_para = max(max(max(hough_space))); %返回值就是这个矩阵的最大值
index = find(hough_space>=max_para*thresh );
%一个矩阵中，想找到其中大于max_para*tresh数的位置,找到可信度较高的圆
% index 是一个一维数组每个值代表着 hough_space(a,b,r) 一维化后的位置
length = size(index); %符合阈值的个数
hough_circle = false(m,n);
%通过位置求半径和圆心.根据数据 一维化后的位置 定位三维的坐标 index=m*n*(r-1)+(a-1)*m+b
for i=1:ecount
for k=1:length
par3 = floor(index(k)/(m*n))+1; %转化为圆的半径
par2 = floor((index(k)-(par3-1)*(m*n))/m)+1;%转换为圆心的a值
par1 = index(k)-(par3-1)*(m*n)-(par2-1)*m; %转换为圆心的b值
if((rows(i)-par1)^2+(cols(i)-par2)^2<(r_min+(par3-1)*step_r)^2+5&...
(rows(i)-par1)^2+(cols(i)-par2)^2>(r_min+(par3-1)*step_r)^2-5)
hough_circle(rows(i),cols(i)) = true;
%为了判断原图中的圆边界点 以便建立 hough_circle ，判断原理是基于欧式距离，允许的误差范围是sqrt（5）。
end
end
end
% 打印结果
for k=1:length
par3 = floor(index(k)/(m*n))+1;
par2 = floor((index(k)-(par3-1)*(m*n))/m)+1;
par1 = index(k)-(par3-1)*(m*n)-(par2-1)*m;
par3 = r_min+(par3-1)*step_r;
fprintf(1,'Center %d %d radius %d\n',par1,par2,par3);
para(:,k) = [par1,par2,par3];
end
阈值处理
f = imread('Fig17.tif');
T = 0.5 * (double(min(f(:))) + double(max(f(:))));
done = false;
while ~done
g = f >= T;
Tnext = 0.5 * (mean(f(g)) + mean(f(~g)));
done = abs(T - Tnext) < 0.5;
T = Tnext;
end
subplot(1, 2, 1), imshow(f), title('原图');
subplot(1, 2, 2), imshow(g), title('利用迭代策略');
表示与描述
单元数组与结构
表示区域设计到两个基本选择
用外部特征（区域的边界）表示区域
用内部特征（组成区域的像素）表示区域
函数floor(x):取最小的整数
floor(2.8)=2,floor(3.8)=3
函数ceil(x):取最大的整数
ceil(3.2)=4,ceil(5.18)=6
函数g = imfill(fI, conn, 'holes')
将填充输入灰度图像fI 的孔洞
gB = imfill(fB, locations, conn)
在输入二值图像fB的背景像素上从参数locations指定的点开始，执行填充操作。
gI = imfill(fI, conn, 'holes')
函数 [L,num]=bwlabel(f,conn)
bwlabel计算二值图像中所有的连通分量（区域）.[L,num]=bwlabel(f,conn),其中f是输入图像，conn指定了期望的连通性（4连接或8连接，默认为8连接）。num是找到的连通分量数，L是标记矩阵，L对每个连通分量分配唯一的1到num的整数
函数find 可以和 bwlabel一起使用，返回构成某个指定对象的像素的坐标向量。
[gB, num] = bwlabel(fB)
产生了多个连续区域（即num > 1）
则使用一下语法可以获得第二个区域的坐标:
[r, c] = find(g == 2)
函数 g=bwperim(f,conn)
g=bwperim(f,conn)这个函数返回二值图像g，其中仅包含f中所有区域的周界（边界）像素。这个函数中conn指定背景的连通性：4连接（默认）或8连接。这样，为了得到4连接的区域边界，可以把conn指定为8,8连接的边界可通过将conn指定为4得到
函数 z=sortrows(S)
对数组进行排序
z = sortrows(S)
S 必须是矩阵或列向量
函数 [z, m, n] = unique(S, 'rows')
既对数组进行排序，又去除重复行
函数 z=circshift(S, [ud lr])
z = circshift(S, [ud lr])
对数组进行向上、向下或侧移指定位置数的移位操作
ud 是 S 向上或向下移位的元素数
若S 是一幅图像，则circshift 就是对图像进行的卷动操作或平移操作。
函数 B = bwboundaries(f, conn, options)
第一版书的语法是 B = boundaries(f, conn, dir)
第二版书已改成 B = bwboundaries(f, conn, options)
conn相对于边界本身，并且值为4或8（默认），参数options的值可以为‘holes’和‘noholes’.使用holes会提取区域和孔洞的边界，也可以提取嵌套在区域内的区域边界。使用noholes只能得到区域或其子区域的边界。首先在B中列出区域，紧跟着是孔洞。
输出的B是px1的单元数组，其中，P是物体数。单元数组中每个单位都包含一个npx2的矩阵，其中的行是边界像素的行列坐标，np是相应区域的边界像素数。每个边界坐标都是以顺时针方向安排的，并且边界的最后一点与第一个点相同，这样就提供了闭合的边界。
[B,L]=bwboundaries(...)L是标记矩阵（与f的尺寸一样），使用不同的整数来标记f的每个元素（无论是区域还是孔洞）。背景像素标记为0,。区域和孔洞数由max(L(😃)给出。
[B,L,NR,A]=bwboundaries(...)返回找到的孔洞数NR,逻辑稀疏矩阵A，它详细描述了父-子-孔洞的依赖关系。由B{k}闭合的更直接的边界由下列语句给出：boundaryEnclosed=find(A(:,k))或boundaryEnclosing=find(A(k,:)).
clc;
clear all;
I=imread('coin.png');
figure;
imshow(I);
title('原始图像');
bw=im2bw(I,graythresh(I)); %图像二值化
figure;
imshow(bw);
title('二值图像');
[b,l]=bwboundaries(bw,'noholes'); %搜索物体的外边界
figure;
imshow(label2rgb(l,@jet,[.5 .5 .5])); %以不同的颜色标记不同的区域，0像素点用0.50.050.5标记
hold on;
for k=1:length(b)
boundary=b{k}; %取出每一个线条的坐标画在在图像上
plot(boundary(:,2),boundary(:,1),'w','LineWidth',3);
end
原始图像：
参数为holes图像 提取区域和孔洞的边界也提取包含嵌套区域的区域边界
参数为noholes图像 仅提取区域和子区域的边界
函数 b8 = bound2eight(b)
b8 = bound2eight(b)
函数 b4=bound2four(b)
去除除了4链接外的像素
函数 g = bound2im(b, M, N, x0, y0)
g = bound2im(b, M, N, x0, y0)
生成一幅二值图像g ，参数 x0 y0 决定图像中b 的最小x 和y 坐标位置
函数[s, su] = bsubsamp(b, gridsep)
[s, su] = bsubsamp(b, gridsep)
输出s 是一个比b有更少的点的边界，点的数目由gridsep的值确定，su 是按比例取得的边界点的集合，这样可使坐标的转换趋于一致。
函数 connectpoly(x, y)
将未连接的点连接起来，并且比原来更加平滑。
x, y 为坐标，可以为数组，但是要一一对应
f = imread('coins.png');
f = im2bw(f);
f = imfill(f, 4, 'holes');
g = boundaries(f, 4, 'cw');
d = cellfun('length', g);
[max_d, k] = max(d);
v = g{k(1)};
[w, h] = size(f);
subplot(1, 3, 1), imshow(bound2im(v)), title('原图');
v = bsubsamp(v, 2);
subplot(1, 3, 2), imshow(bound2im(v)), title('gridsep 2');
subplot(1, 3, 3), imshow(bound2im(connectpoly(v(:, 1), v(:, 2)))), title('connectpoly');
链码
函数 qtdecomp(B, threshold, [mindim maxdim])
将图像分割成等大小的部分方块
函数 qtgetblk(B, Q, middim)
获得四叉树中的分解的块
f = imread('Fig11.tif');
h = fspecial('average', 9);
f = imfilter(f, h, 'replicate');
f = im2bw(f, 0.5);
l = f(1:512, 1:512);
b = bwperim(l, 8);
Q = qtdecomp(b, 0, 2);
[vals, r, c] = qtgetblk(b, Q, 2);
l = zeros(512, 512);
for i=1:length(r)
l(r(i), c(i)) = vals(1, 1, i);
l(r(i) + 1, c(i)) = vals(2, 1, i);
l(r(i), c(i) + 1) = vals(1, 2, i);
l(r(i) + 1, c(i) + 1) = vals(2, 2, i);
end
subplot(1, 2, 1), imshow(f), title('原图');
subplot(1, 2, 2), imshow(l), title('分块后');
函数 c = fchcode (b, conn, dir)
函数 inpolygon(X, Y, xy, yx)
f = imread('Fig12.tif');
subplot(3, 3, 1), imshow(f), title('原图');
B = im2bw(f);
subplot(3, 3, 2), imshow(B), title('二值图');
b = bwboundaries(B);
b = b{1};
[M, N] = size(B);
bim = bound2im(b, M, N, min(b(:, 1)), min(b(:, 2)));
subplot(3, 3, 3), imshow(imdilate(bim, ones(3, 3))), title('boundaries');
[x, y] = minperpoly(B, 2);
b = connectpoly(x, y);
bim= bound2im(b, M, N, xmin, ymin);
subplot(3, 3, 4), imshow(imdilate(bim, ones(3, 3))), title('MMP 2');
[x, y] = minperpoly(B, 3);
b = connectpoly(x, y);
bim= bound2im(b3, M, N, xmin, ymin);
subplot(3, 3, 5), imshow(imdilate(bim, ones(3, 3))), title('MMP 3');
[x, y] = minperpoly(B, 4);
b = connectpoly(x, y);
bim= bound2im(b, M, N, xmin, ymin);
subplot(3, 3, 6), imshow(imdilate(bim, ones(3, 3))), title('MMP 4');
[x, y] = minperpoly(B, 8);
b = connectpoly(x, y);
bim = bound2im(b, M, N, xmin, ymin);
subplot(3, 3, 7), imshow(imdilate(bim, ones(3, 3))), title('MMP 8');
补充之前的
interp1q(x, y, xi) 线性内插
x 为列向量 设置x轴各个点
y 为列向量 设置y轴的各个点
xi 为列向量，代表取xi的位置的值
z = interp1q([0 10]', [0 5]', [0 1 2]') % 返回[0; 0.5; 1.0]
图像的代数运算
imadd,imsubtract,immultiply,imdivide
加法可以提高信号噪音比，乘法可以用1和0组成的掩膜图遮住图像的指定部分，减法可以用于背景去除和运动目标检测等
图像文件I/O函数命令
imfinfo 返回图形图像文件信息
语法：info=imfinfo(filename,fmt)
info=imfinfo(filename)
imread 从图像文件中读取（载入）图像
语法：A=imread(filename,fmt)
[X,map]=imread(filename,fmt)
[...]=imread(...,idx) (CUR,ICO,and TIFF only)
[...]=imread(...,'frames',idx) (GIF only)
[...]=imread(...,ref) (HDF only)
[...]=imread(...,'BackgroundColor',BG) (PNG only)
[A,map,alpha] =imread(...) (ICO,CUR,PNG only)
imwrite 把图像写入（保存）图像文件中
语法：imwrite(A,filename,fmt)
imwrite(X,map,filename,fmt)
imwite(...,Param1,Val1,Param2,Val2...)
imcrop 剪切图像
语法：I2=imcrop(I)
imresize 改变图像大小
语法：B=imresize(A,m,method)
imrotate 旋转图像
语法：B=imrotate(A,angle,method) B=imrotate(A,angle,method,'crop')
rgb和lab空间的相互转换
（1）RGB转XYZ
假设r,g,b为像素三个通道，取值范围均为[0,255]，转换公式如下：
M=
0.4124,0.3576,0.1805
0.2126,0.7152,0.0722
0.0193,0.1192,0.9505
等同于如下公式：
X = var_R * 0.4124 + var_G * 0.3576 + var_B * 0.1805
Y = var_R * 0.2126 + var_G * 0.7152 + var_B * 0.0722
Z = var_R * 0.0193 + var_G * 0.1192 + var_B * 0.9505
上面的gamma函数，是用来对图象进行非线性色调编辑的，目的是提高图像对比度。这个函数不是唯一的，但是我在网上查到的基本都使用上式。
（2）XYZ转LAB
L，a，b*是最终的LAB色彩空间三个通道的值。X，Y，Z是RGB转XYZ后计算出来的值，Xn，Yn，Zn一般默认是95.047，100.0，108.883
rgb2lab.m
function Image=rgb2lab(Image)
% %判断输入的图像是否为3维
if size(Image,3) ~= 3
error('image must have three color channels');
end
%判断输入的类是否为double 类
if ~isa(Image,'double')
Image=double(Image)/255;
end
% Undo gamma correction
R = invgammacorrection(Image(:,:,1));
G = invgammacorrection(Image(:,:,2));
B = invgammacorrection(Image(:,:,3));
% Convert RGB to XYZ
% Y = inv（X）返回矩阵X的倒数
T = inv([3.2406, -1.5372, -0.4986; -0.9689, 1.8758, 0.0415; 0.0557, -0.2040, 1.057]);
Image(:,:,1) = T(1)*R + T(4)*G + T(7)*B; % X
Image(:,:,2) = T(2)*R + T(5)*G + T(8)*B; % Y
Image(:,:,3) = T(3)*R + T(6)*G + T(9)*B; % Z
% Convert XYZ to CIE L*a*b*
WhitePoint = [0.950456,1,1.088754];
X = Image(:,:,1)/WhitePoint(1);
Y = Image(:,:,2)/WhitePoint(2);
Z = Image(:,:,3)/WhitePoint(3);
fX = f(X);
fY = f(Y);
fZ = f(Z);
Image(:,:,1) = 116*fY - 16; % L*
Image(:,:,2) = 500*(fX - fY); % a*
Image(:,:,3) = 200*(fY - fZ); % b*
end
function R = invgammacorrection(Rp)
R = zeros(size(Rp));
i = (Rp <= 0.0404482362771076);
R(i) = Rp(i)/12.92;
R(~i) = real(((Rp(~i) + 0.055)/1.055).^2.4);
end
function fY = f(Y)
fY = real(Y.^(1/3));
i = (Y < 0.008856);
fY(i) = Y(i)*(841/108) + (4/29);
end
function Y = invf(fY)
Y = fY.^3;
i = (Y < 0.008856);
Y(i) = (fY(i) - 4/29)*(108/841);
end
lab2rgb.m
lab转rgb图像相当于以上公式的逆推
function Image = lab2rgb(Image)
if size(Image,3) ~= 3
error('im must have three color channels');
end
if ~isa(Image,'double')
Image = double(Image)/255;
end
WhitePoint = [0.950456,1,1.088754];
% Convert CIE L*ab to XYZ
fY = (Image(:,:,1) + 16)/116;
fX = fY + Image(:,:,2)/500;
fZ = fY - Image(:,:,3)/200;
Image(:,:,1) = WhitePoint(1)*invf(fX); % X
Image(:,:,2) = WhitePoint(2)*invf(fY); % Y
Image(:,:,3) = WhitePoint(3)*invf(fZ); % Z
% Convert XYZ to RGB
T = [3.2406, -1.5372, -0.4986; -0.9689, 1.8758, 0.0415; 0.0557, -0.2040, 1.057];
R = T(1)*Image(:,:,1) + T(4)*Image(:,:,2) + T(7)*Image(:,:,3); % R
G = T(2)*Image(:,:,1) + T(5)*Image(:,:,2) + T(8)*Image(:,:,3); % G
B = T(3)*Image(:,:,1) + T(6)*Image(:,:,2) + T(9)*Image(:,:,3); % B
end
function Y=invf(fY)
Y=fY.^3;
i=(Y<0.008856);
Y(i) = (fY(i) - 4/29)*(108/841);
end
function Rp = gammacorrection(R)
Rp = zeros(size(R));
%临界值当X=0.04045时，gamma(x)= 0.0031306684425005883
i = (R <= 0.0031306684425005883);
Rp(i) = 12.92*R(i);
% x=（gamma(x)^(1/2.4))*1.055-0.055
Rp(~i) = real(1.055*R(~i).^0.416666666666666667 - 0.055);
end
clc,clear,close all
im = imread('bee.jpg');
im1 = rgb2lab(im); % RGB转化为LAB
figure('color',[1,1,1])
%subplot(121),
imshow(im1,[]);title('RGB转成LAB')
figure;
im2 = lab2rgb(im1);
%subplot(122),
imshow(im2,[]);title('LAB转回RGB')
根据公式借鉴网上别人写的代码，发现原rgb图像与图像经lab转回来的rgb图像还是有点差异，有可能是提供的参数精确度的问题。
