Matlab 图像处理编程实践初步

Matlab 图像处理编程实践初步

课程时间：大一学期短学期，7.3~7.10

matlab 自带帮助，可用于语法、特殊符号等查询

本文章作用：巩固梳理知识点，标注学习重难点

0x00 基础知识&矩阵运算

0x01 Matlab 基础知识

Matlab 本质上是一种方便编写代码的工具（矩阵等），类似于 IDE，不过里面用于编写脚本、函数之类的东西，语法也类似于 python。

使用界面如上图，左边一列是当前文件目录，中间一列包括函数脚本编辑器以及命令行窗口，右边一列可以实时显示各变量的值。

数值

与 C 语言类似，可以带小数点以及正负号。可以采用字符 e（或 E）的形式表示 10 的幂，即科学计数法的表示。支持虚数，末尾加 i 或 j，如 2i，-3i 等。

不同形式下，变量会显示不同格式。（默认状态下：没有明确变量而进行命令输入，默认赋值给变量 ans（matlab 自带），因此不能随便把 ans 当作变量名进行使用）

变量

• 变量名区分大小写
• 变量名最多包含 63 个字符，超出忽略
• 变量名必须以字母开头，之后任意数字 + 字母 + 下划线
• 变量名中不允许出现标点符号（很多标点符号在 matlab 中有特殊意义）
• 变量不需要进行提前声明，且可以随时更改类型
• 系统关键字不能作为变量名，在 matlab 命令行中输入 iskeyword 可以调出所有的系统关键字
• 系统中的特殊变量

注：只有在数值后面加上 i 和 j 才表示虚数，其余情况可以视作一般变量处理

表达式

• 运算符优先级与表达式
  类似于 C 语言的运算符优先级。不过因为在 matlab 上使用比较多的还有矩阵，因此会出现矩阵乘法、矩阵加法、矩阵除法（左除、右除）等特殊的运算，这些运算需要满足一些特殊的条件才能运行。

注意：&与&&以及 I 与 II 的区别，一般来说表达式中利用&和 I，这一点与 C 不太一样

程序流程控制语句

• 顺序结构

a=input('a=');
b=input('b=');
c=input('c=');
d=b*b-4*a*c;
x=[(-b+sqrt(d))/(2*a),(-b-sqrt(d))/(2*a)];
disp(['x1=',num2str(x(1)),'x2=',num2str(x(2))]);
简单进行一些说明，input类似于scanf，单引号内的内容是为了更好的交互，需要用户人为输入才能进行变量的赋值。
matlab中的神奇之处在于变量类型可变以及向量和矩阵的快速赋值，如上代码中x变量其实就是一个一维数组（一维向量），利用中括号[...]就可以实现，调用某一个元素利用小括号（这一点与C不一样），同理多维向量（矩阵，立方体...），一般来说向量中变量类型都一致，不一致系统也会进行转换。
disp函数用于在命令行中显示变量内容，num2str表示将变量转化为字符数组。


简单说明一些基础，input类似于scanf，引号内的字符串（居然是单引号）

• 判断结构
  if 语法类似于 VB，与 C 不太一样。

%单判断
if 表达式
  执行程序块
end

%两个选择
if 表达式
  执行程序块1
else 
  执行程序块2
end

%多判断
if 表达式1
  执行程序块1
elseif 表达式2
  执行程序块2
elseif 表达式3
  执行程序块3
else 
  执行程序块4
end
注：程序块不需要加大括号，利用end进行分割，有点类似Pasical

• 分支语句（注意与 C 的区别）

switch 开关表达式
  case 表达式1
     语句段1
  case {表达式2，表达式3，...,表达式m}
     语句段2
  ...
  otherwise
     语句段n
end

• 循环结构(相比于 C 更加简介)

for k = 1:2:100 %变量k从1开始，步长为2到100结束
  for n = 1:100 %变量n从1开始，默认步长为1到100结束循环
    A(k,n) = 1/k+n+1; %A是矩阵，指定元素进行赋值
  end
end

k = 0;
while k<100
  k = k + 2;
end

注：由于背诵语法细节属于浪费时间，我们大可不必纠结细枝末节，必要时可以查看 matlab 的手册，这里罗列最基础也是最常用的知识点。

0x02 矩阵及其运算

在 MATLAB 中,数组是由一组实数和复数排成的矩形阵列。数组是一个二维的“矩形”,习惯上用“行”和“列”来表示二维。从某种意义上来说，标量其实就是一个 1 行 1 列的数组，如果一个数组只有一维的“行”或“列”，那么称之为向量。三维数组可以认为是由若干相同的二维数组堆叠而成，直到更高的维数也是如此。

矩阵创建

• 直接输入法

• 载入外部数据文件

• 利用内置函数创建特殊矩阵

• 利用 M 文件创建和保存矩阵

• 利用冒号生成矩阵

矩阵寻访

• 单值寻访（双下标或者单下标）

• 多元素访问（善于利用冒号）

矩阵拼接

包括水平拼接和垂直拼接

• 基本拼接

注：矩阵拼接必须要满足对应行列相等，否则会报错

• 拼接函数

• cat 函数

• repmat 函数（分块矩阵合体）

矩阵运算

• 矩阵的加减运算

• 矩阵乘法

• 矩阵转置
  m*n 的矩阵 A 的转置矩阵 B 为 n*m 的矩阵，且满足 A(i,j) == B(j,i)
  表达式：B = A'
• 矩阵除法（我测，线代）

常用的是左除 A\B，表示 AX=B 的解，要求 A 和 B 的行数相同。

• 矩阵按位乘除
  大小相等的矩阵才可以进行按位乘除运算。

和 C 完全不一样！！！！！！！！

按位乘除即矩阵对应位置的元素进行相应运算，运算符是 .* .\(按位左除) ./(按位右除)，下图为例子

0x03 自定义函数

在 matlab 中也有类似于自定义函数的功能，不过是依托于函数文件实现的。

具体方法：

1. 在当前目录创建一个函数，命名为 myfunction.m
2. 在函数文件中进行编写

function[s] = myfunction(m)
...

s 为实参，m 为形参，不需要额外的 return。
3. 在主脚本文件中调用函数即可。

0x04 补充

• 分号的使用
  允许行末没有分号，但是会将表达式的值显示在命令行中。加了分号就不会显示了，可以防止冗余变量的占用命令行显示。我们只需要将需要显示的利用 disp 函数显示即可。
• 常用清空函数

>> clc       %清除命令窗口的内容，对工作环境中的全部变量无任何影响
>> close     %关闭当前的Figure窗口
>> close all %关闭所有的Figure窗口
>> clear     %清除工作空间的所有变量
>> clear all %清除工作空间的所有变量，函数，和MEX文件

figure 窗口指的是图形窗口。

0x10 图像处理基础

0x11 数字图像处理概述

图像术语

0x12 Matlab 图像处理基础

文件格式与类型

图像类型包括二值图像、灰度图像、索引图像和 RGB 图像。

• 二值图像：1-Bit Image
  每一个像素点占用一位，即 0，1 。其中 0 是黑，1 是白。（理解：无光是黑，所有颜色汇聚是白）
  存储二值图像：二维矩阵，矩阵元素的值为 0 或 1

• 灰度图像：8-Bit Gray-Level Image
  每一个像素点占用 8 位（1 字节），分辨率为 256，表示颜色的明亮程度。
  存储灰度图像：二维矩阵，矩阵元素的值为 0~255

拓展：打印机打印只有打和不打的区别，即 01，是如何打印出灰度图片的呢？

本质问题：如何利用二值图像表现出灰度图像？

事实上，打印机确实只能打 0 和 1。但是一张灰度图片的像素点，打印机可以用多个像素来表示，即利用一个大的像素块来表示一个灰度值，像素块中 0 越多，整体的灰度就越大。主要的思想是利用空间换取灰度，该算法称为 Dithering（抖动）。

• RGB 图像
  每一个像素占用 3 个字节，分别表示 Red，Green，Blue 的颜色分量。分辨率为 256 *256 *256。
  存储 RGB 图像：三维矩阵，矩阵中每一个元素（x，y，i）取值为 0~255，x，y 是坐标，i：1 2 3 表示红、绿、蓝。
• 索引图像：8 位彩色图像
  每一个像素占用一个字节，表示一个索引值，分辨率为 256。在呈现图片时，系统会根据索引值从色图中选择相应的 RGB 颜色。

对比一下 RGB 图像和索引图像：

发现后者颜色显然没有前者丰富。

如何将 RGB 图像转化为索引图象？

本质：$256256256$ 个桶内元素分配到 256 个桶中，且尽量均衡。

方法一：近似法

R:16 48 80 112 144 176 208 240

G:16 48 80 112 144 176 208 240

B:32 96 160 224

由于人眼对红和绿更加敏感，我们将一个字节拆成 3+3+2 位，红绿占三位，蓝占一位。一共有 8 种红，8 种绿，4 种蓝，共计 $8*8*4=256$ 种颜色。我们按照 256 的等分设置分量，如上所示，对于 RGB 中的一个像素对应的三元组[R,G,B]，我们选择它们对应最接近的进行近似。例如[17,49,89]我们可以近似于[16,48,96]。

方法二：Median-cut algorithm

上述方法存在一定的弊端，如果 Blue 分量接近 32 的几乎没有，那么这一个近似值就相当于浪费了。我们要做到每一个桶内的元素数量尽可能地均衡，于是就产生了中位切分算法。

算法流程：

• 首先得到图片 Red 分量的直方统计图，得到中位数 R1。将 Red 分量小于 R1 的标记为 0，大于等于 R1 的标记为 1。（2）
• 得到所有标记为 0 的 Green 分量的直方统计图，得到中位数 G1。将 Gree 分量小于 G1 的在原先标记后面加 0，大于等于 G1 的在原先标记后面加 1。同理所有标记为 1 的 Green 分量。（4）
• 得到所有标记为 00 的 Blue 分量的直方统计图，得到中位数 B1。将 Blue 分量小于 B1 的在原先标记后面加 0，大于等于 B1 的在原先标记后面加 1。同理标记为 01，10，11 的 Blue 分量。（8）
• 接着红（16），绿（32），蓝（64）；再红（128），再绿（256）。
• 最后得到 256 个不同的位串，原 RGB 图片的每一个 RGB 分量都对应于一个位串中，这样都尽可能保证了每个位串所表示的 RGB 分量均衡。
• 至于每一个位串颜色怎么取，我们可以取平均颜色或者别的方法。

图像处理基本函数

• 文件查询与读取

>> info = imfinfo("2.png")

info = 

  包含以下字段的 struct:

           Filename: 'C:\Users\87632\Desktop\matlab实验\上课实践\2.png'
        FileModDate: '04-Jul-2023 10:11:56'
           FileSize: 15847
             Format: 'jpg'
      FormatVersion: ''
              Width: 397
             Height: 334
           BitDepth: 24
          ColorType: 'truecolor'
    FormatSignature: ''
    NumberOfSamples: 3
       CodingMethod: 'Huffman'
      CodingProcess: 'Progressive'
            Comment: {}

• 图像文件读入

常用的是前两个。根据读入的不同类型（RGB，灰度，二值，索引）产生的矩阵不同，RGB 三维，灰度、二值都是二维，索引需要加入调色板。

• 图像文件写入

常用的也是前两个，fmt 依然可以省略，文件后缀名不能忽略。

• 图像显示
  imshow 函数,一些特殊用法可以查阅帮助手册，一般就是用来展示图象。
  其中 RGB 图像，灰度图像，二值图像用法一样，直接 imshow(myimage);如果是索引图象还需要加入相应的调色板 map，即 imshow(myimage,map)
  注：imshow 函数显示图像的时候，如果矩阵类型是 double，则默认值在[0,1]，超过 1 的就是白色；矩阵类型是 uint8，则默认值在[0,255]。所以进行矩阵转化过程中，为了防止截断通常利用 double 强类型转化，但是最终还是需要转化为 uint8。

图像类型转换

• gray2ind 函数

将灰度图像或二值图像转化为索引图象。

• im2bw 函数

通过设定亮度阈值，将灰度，真彩，索引图像转换为二值图像。

• ind2gray 函数

将索引图像转化为灰度图像。

• ind2rgb 函数

将索引图像转化为真彩图像。

• rgb2gray 函数

将真彩色图像转化为灰度图像或将彩色图像转化为灰色色图。

• rgb2ind 函数

将真彩色图像转化为索引图象。

综合例子

rgbImage = imread("RGB.png"); %读入真彩色图像
grayImage = rgb2gray(rgbImage); %真彩色图像转化为灰度图像
[ind1Image,map1] = rgb2ind(rgbImage,8); %真彩色图像转化为索引图象，色图至多8种颜色
binImage = im2bw(rgbImage); %真彩色图像转化为二值图像，同理灰度、索引
[ind2Image,map2] = gray2ind(grayImage,8);%灰度图像转化为索引图象，色图至多8种颜色
I = ind2gray(ind1Image,map1); %索引图象转化为灰度图像
M = ind2rgb(ind1Image,map1); %索引图象转化为真彩色图像

subplot(3,3,1); imshow(rgbImage); xlabel("原真彩色图像") ;
subplot(3,3,2); imshow(grayImage); xlabel("真彩色转灰度");
subplot(3,3,3); imshow(ind1Image,map1); xlabel("真彩色转索引");
subplot(3,3,4); imshow(binImage); xlabel("真彩色转二值");
subplot(3,3,5); imshow(ind2Image,map2); xlabel("灰度转索引");
subplot(3,3,6); imshow(I); xlabel("索引转灰度");
subplot(3,3,7); imshow(M); xlabel("索引转真彩色");

0x13 Matlab 图像绘制

二维绘图

基本步骤：

• 数据准备
• 设置当前绘图区
• 绘图
• 设置坐标轴和网格线属性
• 标注图形
• 保存和导出图形

例一

clc;clear;close all;
%(1) 数据准备
x = 0:pi/100:3*pi;
y1 = sin(x);
y2 = sin(3*x);
y3 = sin(5*x);
%(2)设置当前绘图区
figure;
%(3)绘图
plot(x,y1,x,y2,x,y3);%一一对应关系
%(4)设置坐标轴和网格线属性
axis([0 7 -1.5 1.5]);%xmin xmax ymin ymax
grid on; %格子线
%(5)标注图形
xlabel('x');
ylabel('y');
title('演示绘图基本步骤');
legend('sin(x)','sin(3x)','sin(5x)');
%右上角的曲线标注，与plot顺序一致

例 2：（作业）

%main.m
clc;clear;close all;
rgbimage = imread("RGB.png");
grayimage = rgbToGray(rgbimage);
imwrite(grayimage,"GRAY.png");
subplot(2,2,1); imshow(grayimage); xlabel("灰度图像");
binimage = rgbToBin(rgbimage);
imwrite(binimage,"BIN.png");
subplot(2,2,2); imshow(binimage); xlabel("二值图像");
subplot(2,2,3); StaticsGray(grayimage); xlabel("灰度图像统计");
subplot(2,2,4);  StaticsBin(binimage); xlabel("二值图像统计");

%rgbToBin.m
function [s] = rgbToBin(m)
   [r,c] = size(m);
   s1 = rgbToGray(m);  %得到灰度图像
   s = s1 > 128; 
end

%rgbToGray.m
function [s] = rgbToGray(m)
  
   s(:,:) = m(:,:,1) * 0.299 + m(:,:,2)*0.587 + m(:,:,3)*0.114;

end %好神奇的写法！！！

%StaticsGray.m
function [] = StaticsGray(m)
  x = 0:255;
  [r,c] = size(m);
  y=zeros(1,256);
  for a = 1:r
      for b = 1:c
          y(m(a,b)+1) = y(m(a,b)+1)+1; %下标从1开始，我们人为+1
      end
  end
  bar(x,y);
  axis([0 255 0 7000]);
end

%StaticsBin.m
function [] = StaticsBin(m)
  x=0:1;
  y = zeros(1,2);
  [r,c]=size(m);
  total = sum(sum(m));
  y(2) = total;
  y(1) = r*c - total;
  bar(x,y);
  axis([-1 2 0 200000]);
end

figure； 呈现一个图形显示窗口。

subplot(m,n,p); 将图像窗口分为 m*n 格，在第 p 格(顺序是从顶向下，从左到右)进行图像呈现，如 plot,imshow,bar 等

注意：如果主脚本文件中有 figure 或者 subplot，此时呈现的图形窗口是主窗口。如果函数文件所要呈现的图形要放于主窗口中，函数文件不需要再 figure，否则会新生成一个 figure 窗口进行覆盖。

plot(x,y); 主要用于绘制曲线，要求 x，y 向量长度一致，保证一一对应。

bar(x,y); 主要用于绘制柱状图，要求 x，y 向量长度一致，保证一一对应。

axis([xmin xmax ymin ymax]); 设置坐标轴，分别是 x 轴最小值、最大值，y 轴最小值、最大值。

sum 函数可以用于计数：如果 A 是一维向量，sum(A)返回的是元素之和；如果 A 是二维矩阵，sum(A)返回的是一个一维向量，对应于每一列的元素之和。因此小技巧 sum(sum(A))可以用于求矩阵元素和。

矩阵可以进行按位比较，A 是二维矩阵，执行 B = (A > 5)，则 B(i,j)=1 当且仅当 A(i,j)>5，否则为 0。不过好像不能用于相等比较。

注意：matlab 中数组下标是从 1 开始而不是 0，需要特别注意！！！另外，尝试多用冒号语法进行一般矩阵的操作，减少循环的使用！（虽然也不是不行）

*​三维绘制

0x20 颜色转化模型及点运算

0x21 颜色模型基础

RGB 模型

HSV 模型

HSI 颜色空间

0x22 颜色模型转换

• RGB 模型和 HSV 模型的转换

• RGB 和 HSI 模型空间的转换

0x23 图像的直方图特征

灰度直方图反映的是一幅图像中各灰度级像素出现的频率之间的关系。以灰度级为横坐标，纵坐标为灰度级的频率，绘制频率同灰度级的关系图就是灰度直方图。

直方图绘制

• bar 函数
  0x13 例 2 中就是利用 bar 函数对灰度图进行统计。
• imhist 函数
  imhist 函数可以直接对灰度图进行频率统计并生成灰度直方图。
  用法：imhist(A);或 imhist(A,n) 。其中 A 是灰度图，n 是灰度级数，可以理解为直方图显示的柱形条数有 n 个（均匀）。

直方图性质

• 灰度直方图只能反映图像的灰度分布情况，而不能反映图像像素的位置，即所有的空间信息全部丢失。
• 一幅图像对应唯一的灰度直方图，反之不成立。不同的图像可对应相同的直方图。如：

0x24 图像的点运算

图像的点运算本质上是一种映射关系。B 是原始图像，A 是生成图像，且 A(x,y)=f(B(x,y))，则 A 是 B 关于函数 f 的映射图像。其中 f 函数称为 GST 函数。根据映射关系不同，点运算可以分为线性点运算，非线性点运算和直方图修正。

线性点运算

对比度指的是一幅图像中明暗区域最亮的白和最暗的黑之间不同亮度层级的测量，差异范围越大代表对比越大，差异范围越小代表对比越小。粗浅地理解就是像素之间颜色的差距。

非线性点运算

直方图均衡化

数字图像处理中，灰度直方图是最简单且最有用的工具。可以说，对图像的分析与观察直到形成一种有效的处理方法，都离不开直方图。

所谓直方图均衡化，就是把原始图像的直方图变换成均匀分布的形式，增加像素灰度值的动态范围，从而达到增强图像整体对比度的效果。

-直方图均衡化不改变灰度出现的次数，改变的是出现次数所对应的灰度级，以避免改变图像的信息结构。

-直方图均衡化力图使等长区间内出现的像素数接近相等。
Matlab 图像处理工具箱提供了 histeq 函数用于实现图像的均衡化处理。其调用格式：

直方图规定化

同样可以利用 histeq 函数实现图像的均衡化处理。同样是力图使对应灰度级接近于用户指定的规定化向量 。

0x25 综合例子

clc;clear;close all;
rgbimage = imread("RGB.png");
subplot(4,4,1); imshow(rgbimage); xlabel("RGB图像");
hsvimage= rgb2hsv(rgbimage);    %RGB模型转化成HSV模型
subplot(4,4,2); imshow(hsvimage); xlabel("HSV图像");
B1 = rgbimage + 50;             %线性点运算1
subplot(4,4,3); imshow(B1); xlabel("灰度增强");
B2 = rgbimage * 2.0;            %线性点运算2
subplot(4,4,4); imshow(B2); xlabel("对比度增强");
B3 = rgbimage * 0.75;           %线性点运算3
subplot(4,4,5); imshow(B3); xlabel("对比度减弱");
B4 = -double(rgbimage) + 255;   %线性点运算4
subplot(4,4,6); imshow(uint8(B4)); xlabel("图像求补");

%非线性点运算
x = 1:255; y = x+x.*(255-x)/255;
subplot(4,4,7); plot(x,y); xlabel("函数的曲线图");
B5 = double(rgbimage) + 0.006 *double(rgbimage).*(255-double(rgbimage));
subplot(4,4,8); imshow(uint8(B5)); xlabel("非线性处理效果");

Gray = rgb2gray(rgbimage); imwrite(Gray,"Gray.tif");
subplot(4,4,9); imshow(Gray); xlabel("灰度图");

%直方图均衡化
B6 = histeq(Gray); 
subplot(4,4,10); imshow(B6); xlabel("均衡化图像");
subplot(4,4,13); imhist(Gray,64); title("原始图像直方图");
subplot(4,4,14); imhist(B6,64); title("均衡化图像直方图");

%直方图规定化 hgram人为规定,hgram(i) = i-1; i:1~256 
hgram = 0:255; J = histeq(Gray,hgram);
subplot(4,4,11); imshow(J); xlabel("直方图规定化");
subplot(4,4,15); imhist(J,64); title("直方图规定化直方统计");

注意：在进行图像运算过程中，由于很多图像像素值是有范围限制的，因此计算结果会存在溢出的现象。在本例的图像求补以及非线性点运算中把图像数据转换为 double 正是避免这种情况的产生。最后需要将矩阵再转化成 uint8。

0x30 图像的基本运算

0x31 图像的代数运算

图像加法

常用于图片的叠加，例如渐变效果或者半透明叠加效果。

• imadd 函数

两个 uint8 数据相加，数值大于 255 时截取（即取 255）

clear all;
I = imread("beauity.png"); I = imresize(I,[384,512]);
J = imread("animal.png");  J = imresize(J,[384,512]);
K = imadd(I,J);
D = imadd(I,50);
subplot(2,2,1); imshow(I); xlabel("美女图像");
subplot(2,2,2); imshow(J); xlabel("野兽图像");
subplot(2,2,3); imshow(K); xlabel("两个图像叠加");
subplot(2,2,4); imshow(D); xlabel("加一个常数");

• imnoise 函数

图像减法

常用于检测运动目标、图像变换等。

• imsubtract 函数

两个 uint8 数据相减，小于 0 截取（即等于 0）

clear all;
I = imread("beauity.png"); I = imresize(I,[384,512]);
J = imread("animal.png");  J = imresize(J,[384,512]);
K = imsubtract(I,J);
D = imsubtract(I,50);
subplot(2,2,1); imshow(I); xlabel("美女图像");
subplot(2,2,2); imshow(J); xlabel("野兽图像");
subplot(2,2,3); imshow(K); xlabel("两个图像相减");
subplot(2,2,4); imshow(D); xlabel("减一个常数");

图像乘法

• immultiply 函数

clear all; close all;
I = imread("beauity.png"); I = imresize(I,[384,512]);
I16 = uint16(I);
I16 = immultiply(I16,I16);
K = immultiply(I,0.5);
subplot(2,2,1); imshow(I); xlabel("美女图像");
subplot(2,2,2); imshow(I16); xlabel("自乘图像");
subplot(2,2,3); imshow(K); xlabel("乘上一个常数");

图像除法

与乘法类似，函数是 imdivide，用于两个图像相除或者除以常数。

图像求补

• imcomplment 函数

clear all; close all;
I = imread("beauity.png"); I = imresize(I,[384,512]);
K = imcomplement(I);
subplot(2,2,1); imshow(I); xlabel("美女图像");
subplot(2,2,2); imshow(K); xlabel("补图");

应用

• 渐变效果
  我们可以利用图像的加法和乘法来实现两个图片转化的渐变效果。不同重叠效果只是两个图片乘上不同系数后再相加。
  例如过渡帧数为 3，图像 A 转化到图像 B.

$1.0 * A + 0.0 B -> 0.75 A + 0.25 * B -> 0.50 * A + 0.50 * B ->0.25 * A +0.75B -> 0.0A + 1.0 *B$

• 给定一张图片和一组掩膜（Musk）图像，实现抠图效果
  掩膜图像是二值图像，我们可以采用图像乘法，因为掩膜图像中黑色区域像素值为 0，做乘法后依旧是 0。白色区域像素值为 1，做乘法后原始图片像素值不变。因而可以达到抠图的效果。

• 如何在视频场景中检测正在运动的物体？
  首先我们获得没有运动物体的背景图片，用新图片减去背景图片即可抠出运动物体。

0x32 图像的几何运算

图像缩放

• 线性插值方法
  图像的缩放其实是会减少或者增加像素点的，那么这些像素点的颜色值我们该如何确定？

若以一维数组的扩充为例，原长度为 4 的数组可以表示成一个折线图。要扩展成长度为 9 的数组，我们将长度 4 进行 9 等分，每个等分点对应于某条折线上的数值，利用数学知识进行计算即可。

那么二维的情况又该怎么办呢？

假设原始图像为 A，x 轴方向缩放系数为 sx，y 轴方向缩放系数为 sy，新图像为 B。则 B 上(x0,y0)位置的像素点需要照应于原图像的像素点。我们让 x0 和 y0 除以对应轴方向的缩放系数得到 x1 和 y1，此时存在 3 种情况。

1. x1 和 y1 都为整数点，那么 B(x0, y0, :)=A(x1, y1, :);
2. 其中一个是整数点，那么对应在线上，按照权重比例插值。如 x1 是整数点，y1 是小数，那么取点 P(x1,floor(y1))，Q(x1,ceil(y1))，显然(x1,y1)夹在 PQ 之间，则 B(x0, y0, :) = A(x1, floor(y1), :) * (ceil(y1) - y1) + A(x1, ceil(y1), :) * (y1 - floor(y1))
3. 两个都是小数点，此时(x1,y1)夹在一个边长为 1 的正方形中。

我们利用上面加权的相似方法求出 R2 和 R1 点颜色值（即 y 是整数），然后再对 R2 和 R1 进行加权求出 P 点颜色值。
具体代码：

%缩放后的坐标(i,j),缩放系数分别是sx和sy，原图像m
function [s] = F(i,j,sx,sy,m)  
  [w,h,l] = size(m);
  dx = min(j/sx + 1,h);
  dy = min(i/sy + 1,w); %反映射到原图像m坐标,防止出现0和越界
  Ax = floor(dx); Ay = floor(dy); %左上角
  Bx = ceil(dx); By = floor(dy);  %右上角
  Cx = floor(dx); Cy = ceil(dy);  %左下角
  Dx = ceil(dx); Dy = ceil(dy);   %右下角

  if(dx == ceil(dx) & dy == ceil (dy)) %在一个点上
      s = m(Ay,Ax,:);
  elseif(dx == ceil(dx)) 
      s = m(Ay,Ax,:) * (Cy - dy) + m(Cy,Cx,:) * (dy - Ay);
  elseif(dy == ceil(dy)) 
      s = m(Ay,Ax,:) * (Bx - dx) + m(By,Bx,:) * (dx - Ax);
  else
      R1 = m(Ay,Ax,:) * (Bx - dx) + m(By,Bx,:) * (dx - Ax);
      R2 = m(Cy,Cx,:) * (Dx - dx) + m(Dy,Dx,:) * (dx - Cx);
      s = R1 * (ceil(dy) - dy) + R2 * (dy  - floor(dy));
  end
end

需要注意的是数组边界问题，不能为 0 也不能越界。因此我们要做一些处理。最后得到的矩阵是 double 类型的，要转成 uint8，否则依托白色。

图像旋转

angle 的单位是度

图像剪切

图像平移

在实际 matlab 中实现平移其实就是简单的赋值。

clear all; close all;
I = imread("beauity.png"); I = imresize(I,[384,512]);
J = zeros(size(I));
H = size(I);
dx = 50; dy =50;
J(dx+1:H(1),dy+1:H(2),1:H(3)) = I(1:H(1)-dx,1:H(2)-dy,1:H(3));
J = uint8(J); %矩阵进行计算后都会变成double，还原成uint8
subplot(1,2,1); imshow(I); xlabel("原始图像");
subplot(1,2,2); imshow(J); xlabel("平移后的图像");

镜像变换

实际 matlab 中还是利用矩阵的赋值操作实现镜像变换。

0x33 图像的仿射变换

两个向量空间之间的仿射变换（或仿射映射）是由一个线性变换接上一个平移组成的。仿射变换可以理解为对坐标进行缩放、旋转、平移后取得的新坐标的值，或者是经过对坐标的缩放、旋转、平移后原坐标在新坐标领域中的值，可以表示为f(x) = Ax+b,其中 A 是变换矩阵，b 是平移矩阵。二维空间中，A 可以按照 4 个步骤分解：尺寸、伸缩、扭曲和旋转。

尺度与伸缩变换

扭曲与旋转变换

综合仿射变换

仿射变换本质

仿射变换求法

0x34 图像的逻辑运算

clear all; close all;
I= imread("RGB.png");
subplot(2,3,1); imshow(I); xlabel("原始图像");
J = imdivide(I,2);
K1=bitand(I,J);   %位与
subplot(2,3,2); imshow(K1); xlabel("位与运算");
K2=bitcmp(I,'uint8'); %位补 即2^8-I
subplot(2,3,3); imshow(K2); xlabel("位补运算");
K3=bitor(I,J);    %位或
subplot(2,3,4); imshow(K3); xlabel("位或运算");
K4=bitxor(I,J);   %位异或
subplot(2,3,5); imshow(K4); xlabel("位异或运算");
K5=bitshift(I,2); %位移位 即 *2^k
subplot(2,3,6); imshow(K5); xlabel("位移运算");

0x40 无损图像压缩

0x41 信息论

物理中有熵的概念，用于描述一个系统的有序性。信息也有类似的概念，我们称之为信息熵。

• 信源（information source）
  需要进行处理的信息。
• 码表（alphabet）
  S={s1,s2,s3...,sn}，信源中所有单位元素都出自码表中。例如只含小写英文字母的信源的码表就是 26 个英文字母。
• 信息熵（entropy of information source）
  用于形容信源信息量的量，信息熵越大，信息量越大，信息越无序。

信息熵的计算方式：
\eta = \sum_{i=1}^{n}P(i)log_{2}\frac{1}{P(i)} = - \sum_{i=1}^{n}P(2)log_{2}P(i)
其中 P(i)表示字母 si 在信源中出现的频率。

一般来说，其他条件相同的情况下，P(i)越均匀，信息熵越大。

而数学家已经证明，信息熵本质上代表了用于表达信源所需要的最小位数（上取整），这就有了之后压缩的想法诞生。

0x42 差分编码

我们可以直接根据图片的信息熵进行编码，但容易看出，这种效率不一定高，如 0x41 中的例子，信息熵为 8，和原始图片的位数相同，这样也就失去了压缩的意义。

差分编码是一种在不改变信源（或者说我们有方法重新得到信源）的前提下，使得信息熵尽可能地小的编码方式。差分编码又称作预测编码（Predictive coding）。原始图像是 A，差分图像是 B，我们人为规定一种预测方式，得到预测矩阵 C，则 B(i,j)=A(i,j) - C(i,j)

常见的预测方式有：

1. 列差分：C(i,j) = A(i,max(j-1,1)); % 取 max 是因为 C 的第一列就是 A 的第一列
2. 行差分：C(i,j) = A(max(i-1,1),j); % 同理，取 max 是因为 C 的第一行就是 A 的第一行
3. 对角差分：C(i,j) = A(i-1,j-1); % 此时初始值应该是第一行和第一列
4. C(i,j) = [A(i-1,j) + A(i,j-1)]/2;
5. C(i,j) = A(i-1,j) + A(i,j-1) - A(i-1,j-1);
6. C(i,j) = A(i-1,j) + [A(i,j-1) - A(i-1,j-1)]/2;
7. C(i,j) = A(i,j-1) + [A(i-1,j) - A(i-1,j-1)]/2;

除非特殊构造的矩阵，一般情况下 B 矩阵的信息熵会比 A 矩阵的信息熵要小。（一般反映在直方图上是比较窄）

灰度图做一遍差分即可，RGB 图像对 R、G、B 分别做一遍。

注：差分编码实际上得到的是一个差分矩阵，尚未进行真正意义上的压缩操作，其目的是降低信息熵，为后续的压缩做准备。

0x43 无损图像压缩

有损压缩和无损压缩的区别在于：能否根据压缩后的图片还原成原始图片。

差分编码得到差分矩阵，我们可以利用 Huffman 编码对差分矩阵进行编码，得到二进制文件，该二进制文件即是压缩文件。

本质上进行 Huffman 编码的过程才是压缩的过程。

Huffman 编码不再赘述。

• Matlab 自带的 Huffman 编码函数

需安装 Communications Toolbox

1. 对字符串进行 Huffman 编码

2. 对图片进行 Huffman 编码

注意：如果要对差分矩阵进行 Huffman 编码的话，码表应该是-255~255。

基于差分编码的无损压缩和解码

1. 读入原始图像 A
2. 对 A 进行差分编码得到差分矩阵 B
3. 对差分矩阵 B 进行 Huffman 编码
   1. 明确码表取值
   2. 求取码表中各元素在矩阵 B 中的频率
   3. 利用 huffmandict 函数得出字典（即构建出 Huffman tree，元素对应于位串），保存字典
   4. 利用 huffmanenco 函数把差分矩阵 B 压缩成二进制一维矩阵 C，将 C 保存成二进制文件
4. 解码
   1. 读取字典和二进制文件
   2. 利用 huffmandeco 函数解码，得到差分矩阵 B
   3. 根据预测方式还原成原始矩阵 A

注：上述提到的 huffman 函数在上面都有例子。

0x50 基于 Haar 小波变换的图像处理

0x51 Haar 小波变换

一维 Haar 小波变换

• 设有一维信号\{ x1,x2\},平均信号a=\frac{x1+x2}{2},细节信号d= \frac{x1-x2}{2}

• 一维信号可以表示成\{a,d\},且原信号还原：x1=a+d;x2=a-d

• 当 x1 和 x2 非常接近的时候，一维信号\{x1,x2\}可近似用\{a\}表示，可实现信号压缩

  • a 可以看成信号的整体信息
  • d 可以看成原信号用 a 表示时丢失的细节信息

• 设有一维多元素信号\{x1,x2,x3,x4\}

• $\begin{align} &a_{1,0}=(x1+x2)/2 \
  &a_{1,1}=(x3+x4)/2 \

\end{align}$ $\begin{align} &d_{1,0}=(x1-x2)/2 \
&d_{1,1}=(x3-x4)/2 \

\end{align}$

• 信号可以表示成\{a_{1,0},a_{1,1},d_{1,0},d_{1,1}\},丢失细节信号可以压缩为\{a_{1,0},a_{1,1}\}
• a_{0,0}=(a_{1,0}+a_{1,1})/2 d_{0,0}=(a_{1,0}-a_{1,1})/2
• 信号可进一步表示为\{a_{0,0},d_{0,0}\},丢失细节信号可以压缩为\{a_{0,0}\}，可以知道a_{0,0}=(x1+x2+x3+x4)/4
• 最终得到离散 Haar 小波信号为\{a_{0,0},d_{0,0},d_{1,0},d_{1,1}\}
  • a：原始信号的低分辨率的近似表达，称为低频信息
  • d：原始信号的低分辨率的近似表达下丢失的细节信息，称为高频信息

• 函数的重建：小波信号还原

• 正变换计算方法 1：

• 正变换计算方法 2：

• 反变换计算方法 1：

• 反变换计算方法 2：

二维 Haar 小波变换

我们现在已经知道了如何对一维信号进行 Harr 小波变换，二维也是类似，将每一行每一列都看成一维的信号处理。

先对所有行进行一次 Harr 小波变换（不用做到底），此时图片变成左右两个部分，左边为图片I_{L}，右边为图片I_{H}

再对I_{L}的所有列进行一次 Harr 小波变换，此时图片变成了上下两个部分，上部分为图片I_{LL}，下边为图片I_{LH}

同理对图片I_{H}也进行类似操作得到两个图片I_{HL},I_{HH}

• 二级小波分解示意图

第二层操作方法和第一层一样。我们称左上角的四个图为 2 级小波，除去这个四个图的为 1 级小波。同理三级小波分解以及 n 级小波分解。

0x52 基于 Harr 小波变换的信号去噪

• 一般噪声特点： （1）高频信号（细节） （2）幅度（数值）小，使用阈值
• 去噪声过程：
  根据需求对原信号进行 1 级或 2 级或 n 级小波分解。
  1 级小波：去除原始信号高频成分（细节）中幅度小于阈值部分
  2 级小波：设定 2 个阈值，称“阈值 1”和“阈值 2”。
  去除 1 级噪声：去除 2 级小波细节分解中小于“阈值 1”部分。
  去除 2 级噪声：去除 2 级小波细节分解中小于“阈值 2”部分。

（去除：赋值为 0）

• 恢复： 对去噪后的小波进行重构，得到去噪后的信号。

0x53 基于 Harr 小波变换的图像压缩

设定了阈值进行滤波就对原始信号发生了改变，这就是有损压缩。阈值越大，压缩比例越大，图像越失真。

通常做三级小波变换效果会比较好。

• 标准分解方法：

• 非标准分解方法：

• 图像的八带分解：

• 图像压缩及解压流程：

1. 读入原始图像 A
2. 进行小波分解，得到矩阵 B
3. 有损压缩：设定阈值，滤波，得到矩阵 C； 无损压缩：C=B
4. 对矩阵 C 进行 Huffman 编码，得到字典和二进制文件并保存
5. 根据字典和二进制文件解压出矩阵 C，矩阵 C 进行小波重构得到原始图像 A（无损压缩）或 有损图像 D（有损压缩）。

因为高频信号逼近于 0（分解层数越多，越逼近），进行小波分解会使得矩阵中接近 0 的像素点数变多，导致信息熵减少。
若进行滤波操作（有损压缩）还会进一步导致信息熵减少。

0x60 图像增强

图像增强是一种基本的图像处理技术，可以通过突出图像中的某些信息，同时抑制或去除某些不需要的信息来提高图像的质量。例如，强化图像的高频分量，可使图像中物体的轮廓更为清晰。图像增强不以图像保真为原则。

其主要目的有两个：

1. 改善图像的视觉效果，提高清晰度。
2. 使图像变得更有利于计算机的处理，如通过锐化处理突出图像的边缘轮廓线，以便进行各种特征分析。

根据增强处理的作用域不同，图像增强可以划分为空域增强方法和频域增强方法两大类。空域增强是指直接在像素域进行处理，而频域增强是指在图像的变换域内对图像进行处理，然后反变换后即可得到增强以后的图像。此外，还有一些特殊的图像增强，称为彩色增强，包括伪彩色增强和真彩色增强。针对不同的图像应采用不同的图像增强方法，或者同时采用几种适当的增强算法进行实验，从中选出视觉效果好、计算量适中的算法。因此，图像增强技术大多数属于试探式和面向问题的。

0x61 灰度增强变换

本质上就是灰度图像点函数的线性变换和非线性变换，详细可见 0x24。

灰度变换函数

• imadjust 函数

• brighten 函数

brighten 函数直接对图床进行操作，在展示完图片代码的后面加上即可。

0x62 空域滤波增强

空域滤波原理

平滑空间滤波器

- 图像在传输过程中，由于传输信道、采样系统质量较差，或受各种干扰的影响，而造成图像毛糙，此时，就需对图像进行平滑处理。平滑可以抑制高频成分，但也使图像变得模糊。

- 平滑空间滤波器用于模糊处理和减少噪声。模糊处理经常用于预处理，例如，在提取大的目标之前去除图像中一些琐碎的细节，桥接直线或曲线的缝隙。

• 线性平滑滤波器（均值滤波器）

• 非线性平滑滤波器

一般来说，中值滤波器适合用于椒盐噪声（黑白噪声）的去除，且模糊程度比线性平滑滤波器低。

• 综合例子

                                    添加椒盐噪声。

                                          添加椒盐噪声。

• medfilt2 函数——中值滤波函数

• filter2 函数——二维线性滤波器函数

• fspecial 函数——创建预定义的二维滤波器——得到一个掩膜矩阵

锐化空间滤波器

离散序列: f(1),f(2),f(3)...,f(n),...

一阶微分: f(2)-f(1),f(3)-f(2),...f(n+1)-f(n),...

二阶微分(一阶微分再进行一次微分): f(3)+f(1)-2f(2),... ,f(n+2)+f(n)-2f(n+1)...

• 拉普拉斯算子

使用拉普拉斯掩膜对图像进行处理本质上是强化了像素之间的差值，得到的是强化后的差值矩阵。

其中(b)——imshow(uint8(abs(f_lps)));相当于画出 f_lps 的幅值。
(c)——imshow(f_lps,[])，所谓**规约化就是最小值映射到 0，最大值映射到 255，**​ **[**​ ]​ ->0​ ~​255

实际操作过程我们可以自己写拉普拉斯算子对原矩阵进行滤波得到 f_lps，f_lps 取绝对值得到拉普拉斯图像，f_lps 规约化得到矩阵 f_lps_G。
f​ _​lps与原始矩阵相加即得增强了的图像。也可以直接利用 filter2 函数实现。

imshow(a,[])本质：
---
b = (a - min(a(:)))./(max(a(:))-min(a(:))); %a 是 double 型矩阵
imshow(uint8(b));

• 使用 filter2 函数进行线性锐化滤波

• 手写 LaplaceFilter 函数

clear all; close all;
I = imread("moon2.tif"); 
I = double(I);
subplot(1,4,1); imshow(I,[]); xlabel("原始灰度图像");
J = LaplaceFilter(I);
subplot(1,4,2); imshow(uint8(abs(J))); xlabel("拉普拉斯图像");
subplot(1,4,3); imshow(J,[]); xlabel("规约化的拉普拉斯图像");
M = I + J;
subplot(1,4,4); imshow(uint8(M)); xlabel("增强后的图像");

%%  3*3的对角拉普拉斯算子
function [s] = LaplaceFilter(m)

dx = [-1,0,1,-1,1,-1,0,1];
dy = [-1,-1,-1,0,0,1,1,1];

[r,c] = size(m);

s =  8 * m;

for y = 1:r
    for x = 1:c
        for k = 1:8
            yy = y + dy(k);
            xx = x + dx(k);
            if(xx <= 0 || yy <= 0 || xx > c || yy >r)
                continue;
            end
            s(y,x) = s(y,x) - m(yy,xx);
        end
    end
end
end

0x70 图像分割与区域处理

0x71 边缘检测

 边缘检测技术对于处理数字图像非常重要，因为边缘是所要提取目标和背景的边界线，提取出边缘才能将目标和背景区分开来。在图像中，边界表明一个特征的终结和另一个特征区域的开始，边界所分开区域的内部特征或属性是一致的， 而不同区域的特征或属性是不同的。边缘的检测正是利用物体和背景在某种图像特性上的差异来实现的，这些差异包括灰度、颜色或纹理特征。边缘检测实际上就是检测图像特性发生变化的位置。经典的边缘提取方法是考察图像的每个像素在某个邻域内灰度的变化，利用边缘邻近一阶或二阶方向导数变化规律，用简单的方法检测边缘，这种方法称为边缘检测局部算子法。​边缘检测的基本思想是通过​检测每个像素元及其邻域的状态，以决定该像素元是否位于一个物体的边界上。如果每一个像素元位于一个物体的边界上，则其相邻像素元灰度值的变化就比较大。假如可以应用某种算法检测出这种变化并进行量化表示，那么就可以确定物体的边界。常用的边缘检测算子主要有罗伯特(Roberts) 边缘算子、索贝尔(Sobel)边缘算子、Prewitt 边缘算子、拉普拉斯(Laplacian)边缘算子、高斯-拉普拉斯(Laplacian ofGaussian)边缘算子和坎尼(Canny) 边缘算子，下面分别介绍。

注意各算子的适用情况以及卷积核

梯度算子

Roberts 算子

Sobel 算子

Prewitt 算子

拉普拉斯(Laplacian)边缘算子

在锐化空间滤波器中我们提到了拉普拉斯算子，其目的是突出图像细节，这些细节很大程度上就有可能是边界。算法过程一样，详细可以看 0x62 锐化空间滤波器。

高斯-拉普拉斯（LOG）算子

• 高斯平滑滤波器

可以使用 fspecial 函数获得高斯滤波器，可以用于去噪。利用 edge 中零交叉检测法和高斯掩膜也可以判断边界。

• 高斯-拉普拉斯（LOG）算子

Canny 算子

总结

0x72 直线提取与边界追踪

Hough 变换提取直线

本质：空间变换

• Hough 变换

• 直角坐标系下 Hough 变换提取直线

• 极坐标系下 Hough 变换提取直线

减少了参数的取值范围。

• Matlab 实现方式

Matlab 实现直线提取内部本质是极坐标系 Hough 变换

在 Matlab 中通过 Hough 变换在图像中检测直线需要三个步骤。

1. 利用 hough()函数执行 Hough 变换，得到 Hough 矩阵
2. 利用 houghpeaks()函数在 Hough 矩阵中寻找峰值点（前 k 大）
3. 利用 houghlines()函数在前两步结果的基础上得到原二值图像中的直线信息。

对于一个彩色或者灰度图像，先进行边缘检测处理得到二值图像 BW。

上述各函数调用方法（具体可以参看帮助手册）：

[H,theta,rho] = hough(BW);

[H,theta,rho] = hough(BW,Name,Value);

H 是变换得到的 Hough 矩阵（即上述提到的统计矩阵，只不过是在极坐标系下的，行表示距离 rho，列表示角度 theta），theta 和 rho 分别对应于 Hough 矩阵每一列和每一行的\theta和\rho组成的向量。或者可以理解成\theta和\rho的取值范围，整数域。

peaks = houghpeaks(H,numpeaks,Name,Value);

其中 H 是由 hough()函数得到的 Hough 函数，numpeaks 是要寻找的峰值数目，默认为 1；peaks 是一个 Q*2 的矩阵，每行的两个元素分别为某一峰值点在 Hough 矩阵中的行、列索引。Q=numpeaks

lines = houghlines(BW,theta,rho,peaks,Name,Value);

BW 是边缘检测后的二值图像；theta 和 rho 是 Hough 矩阵每一列和每一行的\theta和\rho值组成的向量，由 hough 函数返回；peaks 是一个包含峰值点信息的 Q*2 的矩阵，由 houghpeaks 函数返回；lines 是一个结构体数组，数组长度是找到的直线条数。

lines 数组：point1 表示线段的起点（行,列），point2 表示线段的终点(行、列)，theta 和 rho 表示对应直线在极坐标系下的两个参数。

• 显示 Hough 变换矩阵（统计矩阵）

• 示例：寻找图片最长的三条直线

clear all; close all;
I = imread("RGB.png");
I = imrotate(I,60,'crop'); 
rotI = rgb2gray(I);
BW = edge(rotI,'canny');
[H,T,R] = hough(BW);
P  = houghpeaks(H,5,'threshold',ceil(0.3*max(H(:))));
lines = houghlines(BW,T,R,P,'FillGap',5,'MinLength',7);
figure, imshow(I), hold on

%线段两端点知道可以求模长，按此排序
for i = 1:length(lines)-1
   for j = 1:length(lines)-i
       len_1 = norm(lines(j).point1-lines(j).point2);
       len_2 = norm(lines(j+1).point1-lines(j+1).point2);
       if(len_1<len_2)
           temp = lines(j);
           lines(j) = lines(j+1);
           lines(j+1) = temp;
       end
   end
end
for k = 1:3
   xy = [lines(k).point1; lines(k).point2];
   plot(xy(:,1),xy(:,2),'LineWidth',2,'Color','green');

   % Plot beginnings and ends of lines
   plot(xy(1,1),xy(1,2),'x','LineWidth',2,'Color','yellow');
   plot(xy(2,1),xy(2,2),'x','LineWidth',2,'Color','red');
end

边界跟踪

• bwtraceboundary 函数——边界跟踪

• bwboundaries 函数——边界提取

[B,L]=bwbundaries(BW,conn); %conn：跟踪对象时遵守的连通性

BW 是边缘检测处理后的二值图像。L 是标签矩阵，所有属于区域 k 的点在 L 中对应值为 k（类似于染色）。B 是边界元胞数组，B{i}是一个向量，表示区域 i 的所有边界点的坐标。（B 中不包括背景区域，而 L 中包含区域 0 即背景，需要特别注意）
   利用 matlab 自带的边界跟踪函数标记出面积最大的 3 个区域
利用 L 标签矩阵，可以获得区域 k 的像素点数，按此进行排序。注意背景不要算上。

clear all; close all;
I = imread('rice.png');                 %导入图像
figure(1),
subplot(1,3,1);
imshow(I),title('原始图像')
BW = im2bw(I, graythresh(I));           %生成二值图像
subplot(1,3,2);
imshow(BW),title('二值图像')
[B,L] = bwboundaries(BW,'noholes'); %提取边界，并返回边界元胞数组B 和区域标志数组L
subplot(1,3,3);
imshow(label2rgb(L, @jet, [.5 .5 .5]))  %以不同的颜色标志不同的区域
title('彩色标记图像')
hold on
area_ids = unique(L); %离散化处理，L中元素0~k-1对应到1~k，k表示区域数，包括背景，背景看作区域0，即area_ids(1)=0
for i = 1:length(area_ids)
    area{i} = zeros([1,2]);
    area{i}(1) = length(find(L(:)==area_ids(i))); %获得区域i-1的像素点数
    area{i}(2) = i;
end %此时area{1}是背景信息
for i =1:length(area_ids)-1
   for j =1:length(area_ids)-i
       if(area{j}(1)<area{j+1}(1))
           temp = area{j};
           area{j} = area{j+1};
           area{j+1} = temp;
       end
   end
end
for k = 2:4
    boundary = B{area{k}(2)-1};
    plot(boundary(:,2), boundary(:,1), 'w', 'LineWidth', 2)   %在图像上叠画边界
end

放大后确实是最大的三个区域：