目 录试验 1 实验 2 实验 3 斐波那契数列 ..................................................................... 1 图像轮廓线提取技术 ......................................................... 8 RGB 向量空间中的图像分割技术 .................................... 16试验 4 图像的伪装技术 ................................................................. 23 实验 5 用矩阵变换实现对图像形状及颜色畸变的校正........... 34试验 6 Galton 钉板实验 ................................................................ 54 试验 7 迭代与分形 ....................................................................... 63试验 8 微分方程求解 ..................................................................... 71 试验 9 怎样计算平面图形的面积 ................................................. 81 试验 10 数值微积分 ..................................................................... 92试验 11 数值仿真 ........................................................................... 99 附录 MATLAB 软件使用简介 ......................................................... 109试验 1 斐波那契数列试验 1斐波那契数列一、 实验目的与要求 1、认识 Fibonacci 数列，体验发现其通项公式的过程； 2、了解 matlab 软件中进行数据显示与数据拟合的方式； 3、掌握 matlab 软件中 plot, polyfit 等函数的基本用法； 4、提高对数据进行分析与处理的能力。 二、 问题描述 某人养了一对兔，一个月后生育了一对小兔。假设小兔一个月后就可以长 大成熟，而每对成熟的兔每月都将生育一对小兔，且兔子不会死亡。问：一年后 共有多少对兔子？ 三、 问题分析 这个问题，最早由意大利数学家斐波那契(Fibonacci)，于 1202 年在其著 作《珠算原理》中提出。根据问题的假设，兔子的总数目是如下数列： 1，1，2，3，5，8，13，21，34，55，89，144，233，? 问题的答案就是此数列的第 12 项，即一年后共有 144 对兔子。 这个数列通常被称为斐波那契(Fibonacci)数列，研究这个问题就是研究 Fibonacci 数列。把这个问题作更深入的研究，我们会问：第 n 个月后，总共有 多少对兔子？即 Fibonacci 数列的第 n 项是多少？这就需要我们探素 Fibonacci 数列的通项公式。根据问题的描述，我们知道第 n+2 个月后兔子的对数，等于 第 n+1 个月后兔子的对数（表示原来就有的老兔子对数），加上第 n 个月后兔 子的对数（表示生育出来的新兔子对数）。这样就得到关于 Fibonacci 数列的一 个递推公式： Fn? 2 ? Fn ?1 ? Fn 利用 matlab 软件的数据可视化功能将这些数据显示成平面曲线的形式后， 我们可以观察到 Fibonacci 数列的变化规律； 通过 matlab 软件的数据拟合功能， 我们可以大概知道 Fibonacci 数列的函数关系式，结合上面的递推公式，就可以 推导出来 Fibonacci 数列的通项公式。 四、 背景知识介绍 1、数据的可视化。 将离散的数据： F1 , F2 , F3 , F4 ,?, Fn ,? ， 看成平面坐标系里的点： (1, F1 ), (2, F2 ), (3, F3 ), (4, F4 ),?, (n, Fn ),? ， 利用 matlab 软件的 plot 函数在平面坐标系里划出一个点列，就可以实现离1试验 1 斐波那契数列散数据的可视化。plot 函数的基本使用格式为：plot(y)，其中参数 y 表示竖坐标， 即需要显示的数据。 例1 y=1:20;y=y.^3;plot(y) 2、数据的拟合。 数据拟合就是寻找一个目标函数，作为被拟合数据的近似函数关系。目标 函数的类型，可以是多项式、指数函数等。作数据拟合，首先需要估计目标函数 的类型， 这一点可以通过数据可视化来观察得到，而一阶多项式是最常见的目标 函数，此时称为线性回归。确定拟合系数的原则是最小二乘法，即所有误差的平 方和取最小值。在 matlab 软件中以多项式为目标函数作数据拟合的函数是 polyfit，它的基本使用格式为：polyfit (x,y,n)。 其中(x,y)是被拟合的数据，即平面上的一个点列，而 n 是事先确定的多项 式的阶数。 例2 x=[1,3,4,5,6,7,8,9,10];y=[10,5,4,2,1,1,2,3,4]; polyfit (x,y,2) 结果： 0.2676t 2 - 3.6053t + 13.4597 3、数列的通项公式。 寻找一个整标函数，使得它在 n 处的函数值，等于数列的第 n 项的值，这个 函数就是数列的通项公式。 4、黄金分割。 把一条线段分割为两部分，使其中一部分与全长之比等于另一部分与这部 分之比（如下图） 。其比值是一个无理数 ( 5 ? 1) ? 2 ，取其前三位数字的近似值 是 0.618。由于按此比例设计的造型十分协调美观，因此称之为黄金分割。AM:AB=MB:AM五、 实验过程 本试验将 Fibonacci 数列的有限项， 看成是待处理的数据。 首先利用 matlab 软件的可视化功能， 将这些数据显示在平面坐标系中， 观察其图形类似什么曲线， 结论是：指数函数的曲线。进一步，利用指数函数与对数函数的互逆关系，将原 有数据取对数，再观察其曲线形状是否类似直线，以验证原来的观察是否正确。 通过观察到的目标函数， 然后利用 matlab 软件的数据拟合功能， 得到 Fibonacci 数列通项公式的近似关系。最后，从近似关系出发，推导出来 Fibonacci 数列的 通项公式。 1、观察数据的大概函数关系。 为了研究 Fibonacci 数列的变化规律，我们取此数列的前 30 项来观察。利 用 Matlab 软件的数据可视化功能，将这些数据显示在平面坐标系中，观察其中 蕴涵的函数关系。具体的实现流程为：（1）定义数组 fn；（2）显示数组 fn。 具体的代码如下： function plotfibo(n) %定义函数显示 Fibonacci 数列前 n 项2试验 1 斐波那契数列fn=[1,1];%将数列的前两项放到数组 fn 中for i=3:n %fn 的第 3 项到第 n 项 fn=[fn,fn(i-2)+fn(i-1)]; %将第 i 项添加到数组 fn 中 end %循环结束 plot(fn) %将装有数列前 n 项的数组显示出来 这个函数的调用方式是：plotfibo(30)，显示出来的图像为图 1，经观察， 觉得曲线的形状象指数函数的曲线，其数据无限增大。可以改变参数 n 的值，反 复观察。图 1 n=30图 2 n=50图 3 n=500图 4 n=10002、进一步验证上一步得到的结论。 经过上一步的观察，觉得这些数据应该是指数函数的形式。为了进一步验 证这个结论是否正确， 可以利用指数函数与对数函数的互逆关系。如果这些数据 确实是指数函数的形式，则经过取对数后应该是一个线性关系，即一阶多项式， 从图形上看应该象一条直线。因此，再利用 Matlab 软件的数据可视化功能，将 这些数据取对数后显示在平面坐标系中，观察它是否象一条直线。具体的实现流 程为：（1）定义数组 fn；（2）数组 fn 取对数；（3）显示数组 fn。具体的代 码如下： function plotlnfibo(n) %显示取对数后的前 n 项 fn=[1,1]; %将数列的前两项放到数组 fn 中3试验 1 斐波那契数列for i=3:n%fn 的第 3 项到第 n 项fn=[fn,fn(i-2)+fn(i-1)]; %将第 i 项添加到数组 fn 中 end %循环结束 fn=log(fn) %将原来的数据取对数 plot(fn) %将装有数列前 n 项的数组显示出来 这个函数的调用方式是： plotlnfibo(30)， 显示出来的图像为图 5， 经观察， 觉得它确实象一条直线。可以改变参数 n 的值，反复观察。图 5 n=30图 6 n=50图 7 n=500图 8 n=10003、获得数据的近似关系式。 经过以上第一步的观察，确定 Fibonacci 数列的数据是指数函数的关系， 第二步验证了第一步得到的结论， 因此我们认为 Fibonacci 数列的数据关系就是 指数函数，取对数后就是线性函数，即一阶多项式。利用 Matlab 软件的数据拟 合功能，通过取对数后的数据，用一阶多项式拟合出它的函数关系式，可以得到 Fibonacci 数列通项公式的一个近似表达式。具体的实现流程为：（1）定义数 组 fn； （2）数组 fn 取对数； （3）用一阶多项式拟合数组 fn。具体的代码如下： function fitlnfibo(n) %根据取对数后的数据，拟合出线性表达式 fn=[1,1]; %将数列的前两项放到数组 fn 中 for i=3:n %fn 的第 3 项到第 n 项 fn=[fn,fn(i-2)+fn(i-1)]; %将第 i 项添加到数组 fn 中4试验 1 斐波那契数列end%循环结束xn=1:n; %定义横坐标 fn=log(fn) %将原来的数据取对数 polyfit(xn,fn,1) %拟合装有数列前 n 项的数组 这个函数的调用方式是：fitlnfibo(30)，运行后返回结果是：0.4799， -0.7768。这两个数据就是一阶多项式的系数，即： log( Fn ) ? ?0.9n 为了提高精度，可以加大 n 的值。取 n＝1000 时得到： log( Fn ) ? ?0.2n 从上面的表达式，可以得到数列通项公式的近似： Fn ? 0.4476 ? 1.6180n 4、推导 Fibonacci 数列的通项公式。 通过以上的观察和分析，我们知道 Fibonacci 数列的数据大概是指数函数 的关系。因此，我们猜测它的通项公式具有形式： Fn ? k ? r n 。将这个表达式代 入递推公式 Fn? 2 ? Fn ?1 ? Fn 中，得到： k ? r n?2 ? k ? r n?1 ? k ? r n 。经过简化得到：r2 ? r ?1 这是一个一元二次的代数方程，其两个根形式如下： r ? (1 ? 5) ? 2 考虑到 Fibonacci 数列的数据无限增大，我们取 r ? (1 ? 5) ? 2 ，于是得到通项公式如下：Fn ? k ? [(1 ? 5) ? 2]n上面的公式对吗？我们可以来验证。取 n=1 和 n=2 代入上面的公式中计算， 显然得不到 F1 ? 1, F2 ? 1 ，因此它不是 Fibonacci 数列的通项公式。 但这个公式并非一无是处，我们可以来考虑这个公式与 Fibonacci 数列到 底相差多少。因此，我们引入以下一个数列： Tn ? Fn ? k ? [(1 ? 5) ? 2]n 可以验证，这个新的数列也满足同样的递推公式： Tn ? 2 ? Tn ?1 ? Tn ，因此我 们猜测它同样是指数函数的形式，可以假设其表达式为： Tn ? ? ? r n ，代入递推 公式后，同样可以得到： r 2 ? r ? 1 。这里的 r 显然不同于上面的 r，故这个 r 取 值为： r ? (1 ? 5) ? 2 ，从而得到： Tn ? ? ? [(1 ? 5) ? 2]n 。故有： Fn ? k ? [(1 ? 5) ? 2]n ? ? ? [(1 ? 5) ? 2]n 由条件 F1 ? 1, F2 ? 1 确定其中的常数，得到： Fn ? {[(1 ? 5) ? 2]n ? [(1 ? 5) ? 2]n } ? 5 可以证明，它确实就是 Fibonacci 数列的通项公式。 这个公式是法国数学家比内(Binet)在 1843 年发现的，称为比内公式。 六、 结论与应用5试验 1 斐波那契数列1、Fibonacci 数列的阶。 通过以上试验过程，我们得到了 Fibonacci 数列的通项公式，它类似一个 指数函数，当 n 无限增大时，其数据也无限增大，变化的阶是：{[(1 ? 5) ? 2]n } ? 5在 Fibonacci 数列的通项公式中，出现了 (1 ? 5) ? 2 和 (1 ? 5) ? 2 等无理 数，而它们运算以后的结果确是正整数，多么奇妙啊。 2、Fibonacci 数列与黄金分割数的关系。 上面的两个无理数之间，存在着这样的关系： (1 ? 5) ? 2 ? [( 5 ? 1) ? 2]?1 而 G ? ( 5 ? 1) ? 2 就是著名的黄金分割数。因此，Fibonacci 数列的通项公 式又可以写成如下形式：Fn ? [G ? n ? (?1) n ?1 G n ] ? 5 可以验证，Fibonacci 数列与黄金分割数之间还有如下的关系：F2 F4 F F F F ? ? ? ? 2 n ? ? ? G ? ? ? 2 n ?1 ? ? ? 3 ? 1 F3 F5 F2 n ?1 F2 n ? 2 F4 F2Limn ??Fn 5 ?1 ? ?G Fn ?1 23、Fibonacci 数列通项公式的其它形式。 Fibonacci 数列的通项公式还有其它形式，比如：1 1 Fn ?1 ? 0 0?1 1 1 00 ?1 1 00 0? ?0 0 0 10 0 0 1?1 ? 0 ?? ? ? ? ? ? ?4、自然界中的 Fibonacci 数列。 Fibonacci 数列与自然界中的许多现象有着紧密的联系，比如：植物的分枝 生长、向日葵种子的排列、花瓣的数量、晶体的结构等。花瓣的数量，一般都是 Fibonacci 数6试验 1 斐波那契数列以斐波那契螺旋方式的排列如果顺时针与逆时针螺旋的数目，是斐波那契数列中相邻的 2 项，可称其 为斐波那契螺旋，也被称作黄金螺旋。这样的布局能使植物的生长疏密得当、最 充分地利用阳光和空气。 5、应用。 Fibonacci 数列在纯粹数学、运筹优化、 计算机科学等领域具有重大的应用 价值。 本试验研究的是 Fibonacci 数列的变化规律，而数列本质上就是一些数据。 因此，对于一般的数据（比如：从调查中获得的数据、从试验中获得的数据）, 我们也可以参照这样的方式，来分析其中蕴涵的规律。利用 Matlab 软件的数据 可视化功能和数据拟合功能，可以极大地方便我们进行数据处理分析。 七、 练习 1、讨论数列 an?1 ? an ? 1 ? an , a1 ? 1的变化规律。 （1）在平面坐标系中画出数列变化的折线图； （2）观察图形，你认为数列的极限是什么； （3）观察图形，寻找恰当的函数去拟合这个数列； ? 1 ? n 变化规律。 2、讨论调和级数的 n ?1 （1）画出部分和数列 S n 变化的折线图，观察变化规律； （2）引入数列： n ? S 2 n ? S n ，作图观察其变化，猜测是否有极限； H （3）引入数列： n ? S n ，作图观察其变化，寻找恰当的函数拟合； G （4）调和级数的部分和数列的变化规律是什么？2? ?7试验 2 图像轮廓线提取技术实验 2一、实验目的与要求：图像轮廓线提取技术1、 能熟练应用 matlab 去分析问题、解决问题； 2、 熟悉对 matlab 的图像处理的功能，掌握基本的图像处理的若干命令； 3、 在应用 matlab 进行图像处理方面具备一定的编程能力。 4、 掌握 figure， imread， image， colormap， imshow， imwrite， subplot， title，rgb2gray，imfinfo 等语句的基本使用方法。5、 掌握图像轮廓线提取的简单方法并上机实现。 6、 掌握 matlab 自带的一些常用边界检测算子的使用， 提高对复杂图像处理的能力。 二、问题描述 “图像轮廓线提取”是数字图像处理中对图像进行处理和分析之前的一项非 常重要的工作。指的是从原始图像中，以手动或自动的方法，将图片中的人物、 动物、植物或者其他任何对象的（特征）轮廓线提取出来，使之成为一幅独立的 黑白线条图。从而达到将物体与背景分开，物体与物体分开的效果。提取轮廓线 被应用于许多方面，例如人脸检测和跟踪。它结合了认知科学、图象处理、计算 机图形学、机器视觉和模式识别等多个研究领域。 三、问题分析 既然“图像轮廓线提取”的黑白线条图所在位置往往是图像中两区域交界位 置，则可以通过图像特征（如形状、颜色、纹理等）变化情况来检测两区域交界 处。 最简单的方法就是采用阈值检测法，即将当前检测点的特征与周围点的特征 进行比较， 若发现有较大的差异， 则认为当前检测点属于两区域的交界点， 否则， 认为同一区域内的点。 四、背景知识介绍 首先介绍几种基本的图像格式，再介绍一下 matlab 中常见的图像处理命令 及其用法。 ? 常见图像格式1、二值图 单色图像则是带有颜色的图像中比较简单的格式， 它一般由黑色区域和白色 区域组成，可以用一个比特表示一个像素， “1”表示黑色， “0”表示白色，当然试验 2 图像轮廓线提取技术也可以倒过来表示，这种图像称之为二值图像。 2、灰度图 我们也可以用 8 个比特（一个字节）表示一个像素，相当于把黑和白等分为 256 个级别， “0”表示为黑， “255”表示为白，该字节的数值表示相应像素值的 灰度值或亮度值，数值越接近“0” ，对应像素点越黑，相反，则对应像素点越白， 此种图像我们一般称之为灰度图像。单色图像和灰度图像又统称为黑白图像。 3、彩色图 与黑白图像对应，就存在着彩色图像，这种图像要复杂一些，表示图像时， 常用的图像彩色模式有 RGB 模式、CMYK 模式和 HIS 模式，一般情况下我们只 使用 RGB 模式，R 对应红色，G 对应绿色，B 对应蓝色，它们统称为三基色， 这三中色彩的不同搭配， 就可以搭配成各种现实中的色彩，此时彩色图像的每一 个像素都需要 3 个样本组成的一组数据表示， 其中每个样本用于表示该像素的一 个基本颜色。基于三刺激理论（Tristimulus Theory） ，我们的眼睛通过三种可见 光对视网膜的锥状细胞的刺激来感受颜色。这些光的波长为 630nm（红色） 、 530nm（绿色）和 450nm（蓝色）时的刺激达到高峰。通过光源中的强度比较， 我们感受到光的颜色。这种视觉理论是使用三种颜色基色：红、绿和蓝在视频监 视器上显示彩色的基础， 称为 RGB 颜色模型。 以一个常见的例子来说明： Windows 环境下主要的图像格式之一，BMP 灰度图像，以其格式简单，适应性强而倍受 欢迎。这种文件格式就是每一个像素用 8bit 表示，显示出来的图像是黑白效果， 最黑的像素的灰度（也叫作亮度）值为“0” ，最白的像素的灰度值为“255” ，整 个图像各个像素的灰度值随机的分布在“0”到“255”的区间中，越黑的像素， 其灰度值越接近于“0” ，越白（既越亮）的像素，其灰度值越接近于“255” ；与 此对应的是在该文件类型中的颜色表项的各个 RGB 分量值是相等的，并且颜色 表项的数目是 256 个。每一个像素颜色由其红、绿、蓝三色的强度值联合决定。 4、常见的颜色值如下： 颜色 红 绿 蓝 黄 紫 青 白 黑 灰 红 R 255 0 0 255 255 0 255 0 128 绿 G 0 255 0 255 0 255 255 0 128 蓝 B 0 0 255 0 255 255 255 0 128 ? 常用图像处理命令1、图像文件的读取 i=imread(文件名,’图像文件格式’)。将文件名指定的图像文件读入数组 I，I 的数据类型为无符号 8 位整数（uint8） 。对于灰度图像，I 是一个二维数组；对 于真彩色图像，I 是一个三维数组（n×m×3） 。文件名是指定图像文件名称的字 符串，’图像文件格式’是指定图像文件格式的字符串。 如：i=imread('C:\Documents and Settings\ My Documents\My Pictures\a','jpg'); 或直接在第一个参数上加扩展名： i=imread('C:\Documents and Settings\ My Documents\My Pictures\a.jpg'); 注意：后面最好加上“； ”号，因为图像若大则需时间较长。 2、图像的显示 image(i):将矩阵 I 作为图像显示，I 中的每一个元素在图像中表示为一个颜 色块，I 可以是 n×m 阶矩阵或者 n×m×3 的数组，其元素为 double 或 uint8 数9试验 2 图像轮廓线提取技术据类型。 imshow:是在图像显示中最常用，也是功能最强的图像显示函数，常用的格 式为：imshow(i,n),使用 n 个灰度等级显示图像 I。如果默认 n，则使用 256 级或 64 级灰度显示图像。 3、显示灰度直方图 imhist(i),显示灰度图像 I 的灰度直方图。其中 I 是一个二维矩阵，imhist 用 柱状图形式显示图像中每个灰度出现的频率， 即矩阵中每一个元素对应的数值在 整个矩阵中出现的频率。 如：先获取 r 颜色的二维矩阵：j=i(:,:,1); 再运行：imhist(j); 4、获取像素点信息 pixval,在使用 image 或 imshow 显示图像后，在执行 pixval 指令，便可以读 取光标所指向的像素点的坐标和灰度值；如果按住鼠标，在图像中拖曳，便可以 显示光标所移动的距离。 再次运行 pixval 则是关闭获取像素点信息。 5、求矩阵的规模 D=size(i),函数 size 用来求矩阵的规模。如果 I 是 n×m 的矩阵，则 D 是包含 两个元素的向量： D=[n, m ]； 同样道理， 如果 I 是 n×m×3 的矩阵， D=[n,m,3]。 则 如果我们事先知道矩阵的维数，也可以直接用这样的方法求得矩阵的规 模:[n,m]=size(i)或[n, m, k]=size(i)等。 6、数据类型转换 i=uint8(X),这个函数实现将矩阵 X 中的所有元素转换为无符号 8 位整形数 据，把转换的结果保存在矩阵 I 中。这个函数对所有类型的数据都适用，转换后 数值都落在 0 到 255 之间。 i=double(X),这个函数与 uint8（）相似，实现的功能是将矩阵 X 中的所有元 素转换为 double 类型（16 位双精度类型） 。 五、实验过程 根据分析， 图像中位于轮廓线上的点，它与其相邻的点的灰度值差有一定的 跳跃，故通过值的对比，就可以将那些边缘点提取出来。 考虑到图像边框上的像素点对于轮廓线获取的结果影响不大， 且为了代码的 简单实用， 首先去掉图像的边框上的所有像素点， 去掉灰度值矩阵的第一行、 即： 第一列、最后一行和最后一列。这样做的好处是显而易见的，因为可以保证每一 个待比较的点其周围都有 8 个供比较的点，如图示：在算法中，将中心被检测点依次与其上下、左右、左上、右下，和右上、左 下 8 个点做比较， 若差值大于规定值， 则该检测点就是轮廓线上的点， 反之不是。 算法关键的地方是对 IMREAD 生成的矩阵进行了非线性变换。10试验 2 图像轮廓线提取技术算法流程如下： 首先，将这个矩阵从 8 位无符号型变成双精度的数值存储； 然后才能进行下一步的将其每个值都除以 255，这样就把每个数值都调整为 0～1 之间的小数； 再用正弦 sin 函数将这些值做变换，并取求得的结果放大 40 倍。这样一来， 本来是在 0～255 之间的正整数， 就被变换到 0～1 之间的小数， 然后再离散成 0～ 40 之间的双精度数值。其目的是为了在进行比较灰度值的时候，方便自定义各 种差值。 最后，用户给定一个差值，根据这个值来比较检测点与其周围 8 个点的灰度 值，若大于给定差值，则认为检测点位于轮廓线上，否则，不在轮廓线上。 灰度图的源代码（.m 文件）如下： function graydrawout(pix,n) %灰度图的轮廓线提取 A=imread(pix); %读取指定的灰度图 [a,b]=size(A); %a,b 分别等于矩阵 A 的行数和列数 B=double(A); %将矩阵 A 变为双精度矩阵 D=40*sin(1/255*B); %将矩阵 B 进行线性变换 T=A; %新建与 A 同等大小矩阵 for p=2:a-1 %处理图片边框内的像素点 for q=2:b-1 if (D(p,q)-D(p,q+1))&n|(D(p,q)-D(p,q-1))&n| (D(p,q)-D(p+1,q))&n|(D(p,q)-D(p-1,q))&n| (D(p,q)-D(p-1,q+1))&n|(D(p,q)-D(p+1,q-1))&n| (D(p,q)-D(p-1,q-1))&n|(D(p,q)-D(p+1,q+1))&n T(p,q)=0; %置边界点为黑色 else T(p,q)=255; %置非边界点为白色 subplot(2,1,1); %将窗口分割为两行一列，下图显示于第一行 image(A); %显示原图像 title('灰度图原图'); %图释 %保持图片显示比例 subplot(2,1,2); %下图显示于第二行 image(T); %显示提取轮廓线后的图片 title('提取轮廓线'); %图释 %保持图片显示比例 注释：pix－为要提取轮廓线的灰度图名（带路径） ，由单引号括住。 n－自定义的灰度值差值，超过该值就是轮廓线上的点，反之不然。 这是一个最重要的参数，通过调节它的值，修整轮廓线的效果，取范 围在 0～40 之间的任何有理数。 A－MATLAB 读取原图片后返回的数据矩阵，2 维（M N ） T－新建的与 A 矩阵同等行列数的矩阵，待放入比较后结果，代码的主体11试验 2 图像轮廓线提取技术是 if－else－end 部分，其原理就是将中心点与其周围的 8 个点依次比 较，发现有一个差值大于自定义的值时，就判断其为轮廓线上的点并 将之置为黑色；若其与周围 8 个点的比较值都小于自定义值时，则其 不在轮廓线上，置白色。 灰度图提取轮廓线示例（n=20） ：彩色图提取图像轮廓线 彩色图片经过 IMREAD 函数读取后， 由于它的每一个像素点都是由红、 绿、 蓝三色的强度值一起定义其颜色的，所以，生成的是一个 3 维的（M N 3） 矩阵。矩阵平面（：： ，，1）代表对应像素点的红色强度值，矩阵平面（：： ，，2） 是绿色强度值，矩阵平面（：： ，，3）则是蓝色强度值。在每一个矩阵平面中强度 值的范围是[0,255]的正整数。位于轮廓线上的点，它的红、绿、蓝三色的强度值 与其周围的点必然有一定的差值，也正式利用这些差值的比较，可以将那些位于 边缘的点提取出来。 在算法上，彩色图与灰度图不同的是，彩图的每一个单色矩阵可以单独用 来判定该点是否位于轮廓线上。其算法与灰度图的原理非常相似。 在彩色图的算法中，加入了对判定矩阵的选择，即提供使用者自由选择用 图片的哪一个单色矩阵进行轮廓线提取。因为根据图片本身的特点，选择更适合 于进行提取的矩阵， 其效果是十分不同的，在下面的讨论中将会再次提到这个问 题。 算法关键的地方同样是对 IMREAD 生成的矩阵进行了线性变换。首先，将 这个矩阵从 8 位无符号型变成双精度的数值存储； 然后才能进行下一步的将其每 个值都除以 255，这样就把每个数值都调整为 0～1 之间的小数；再用正弦 sin 函 数将这些值做变换，并取求得的结果放大 40 倍。这样一来，本来是在 0～255 之间的正整数，就被变换到 0～1 之间的小数，然后再离散成 0～40 之间的双精 度数值。 其目的是为了在进行比较单色强度值的时候，方便自定义各种带小数的 差值。彩色图轮廓线提取的源代码（.m 文件）如下： function colordrawout(pix,n) %彩色图片轮廓线提取函数 A=imread(pix); %读取指定彩色图片 B=A(:,:,1); %红色强度值矩阵 C=A(:,:,2); %绿色强度值矩阵 D=A(:,:,3); %蓝色强度值矩阵 for i=1:3 %依次从三个矩阵中提取轮廓线 if i==1 %从红色矩阵提取 E=B; else if i==2 %从绿色矩阵提取12试验 2 图像轮廓线提取技术E=C; else E=D;%从蓝色矩阵提取 H=double(E); %将选择的矩阵变为双精度矩阵 F=40*sin(1/255*H); %进行非线性变换 [k,j]=size(B); % k,j 分别为矩阵 D 的行数和列数 T=A; for p=2:k-1 for q=2:j-1 if (F(p,q)-F(p,q+1))&n|(F(p,q)-F(p,q-1))&n| (F(p,q)-F(p+1,q))&n|(F(p,q)-F(p-1,q))&n| (F(p,q)-F(p-1,q+1))&n|(F(p,q)-F(p+1,q-1))&n| (F(p,q)-F(p-1,q-1))&n|(F(p,q)-F(p+1,q+1))&n T(p,q,1)=0;T(p,q,2)=0;T(p,q,3)=0; %置边界点黑色 else T(p,q,1)=255;T(p,q,2)=255;T(p,q,3)=255;%置非边界点白色 subplot(2,2,i+1); %将窗口分割为两行两列，下图显示于第 i＋1 位置 image(T); %显示轮廓线 title(i); %图释 %保持图片显示比例 subplot(2,2,1); %下图显示于第 1 位置 image(A); %显示原彩色图片 title('彩色图原图'); %图释 %保持图片显示比例 注释：pix－为要提取轮廓线的灰度图名（带路径） ，由单引号括住。 n－自定义的强度值值差值，超过该值就是轮廓线上的点，反之不然。 这是一个最重要的参数，通过调节它的值，修整轮廓线的效果 范围为 0～40 之间的任何有理数。 A－MATLAB 读取原图片后返回的数据矩阵，3 维（M N 3） T－新建的与 A 矩阵同等行列数的矩阵，待放入比较后结果。 代码中加入的部分是对单色强度值矩阵的选择。从比较中可以看出，选择 彩色图片中最大程度的颜色矩阵，可以提高提取轮廓线的效果。譬如：图片以红 色为主，就选择 1－红色矩阵。 代码的主体同样是 if－else－end 部分。其原理仍是将中心点与其周围的 8 个点依次比较， 发现有一个差值大于自定义的值时，就判断其为轮廓线上的点并 将 T 矩阵中对应点置为黑色；若其与周围 8 个点的比较值都小于自定义值时， 则其不在轮廓线上，T 矩阵对应点置白色。13试验 2 图像轮廓线提取技术彩色图提取轮廓线示例：注：n 值为 1.6 时， 1）由红色矩阵提取的轮廓线；2）由绿色矩阵提取的轮廓线； 3）由蓝色矩阵提取的轮廓线。 六、结论与应用 在前面的算法中， 在逐行逐列扫描比较各个象素点是否属于轮廓线上的点的 过程中， 事先预设除去边框的所有点， 这是为了得到比较精简实用的算法。 否则， 算法会相对复杂， 首先需将图像的边框引入，在算法中便多了四条特殊的线和四 个特殊点，就是上下左右的边线及图像的四个顶点。拿四个顶点而言，可供它们 进行比较的参考点只有三个， 而且对于每个点而言， 三个参考点的位置都不相同， 这样就造成了算法的累赘冗余。同理，边线也是，边线上的点只有五个参考点， 也是各边方向不同，需要大量相似的代码。 实验证明，提取轮廓线时，加上边框与不加边框的区别非常不明显，故采用 后者的精简算法。 在灰度图和彩色图的算法中，对读取生成的图片数据矩阵进行 了线性变换， 为的是更好的调节自定义强度值的差值， 从而得到令人满意的结果。 选择三角函数 sin，用其进行变换后，可以保证结果值都落到 0～1 之间，从而控 制扩大后的数值都在 0～40 的范围。其实这个函数也可以使用其他离散函数，譬 如 log，它的离散作用更加明显，可以将值都变换到[0，＋ ? )区间上。) 由于噪声和模糊的存在， 轮廓线可能会变宽或在某些点处发生间断。在算法 中，发现图片出现断点，可以通过减小 n 值来提高轮廓线的精度从而减少断点。 但是，轮廓线变粗的同时，弊处是也使原有的清晰的其他线条变得更加粗，甚至 出现了模糊的散点。为此，只能寻找一个合适的 n 值，使断点不至于太多，轮廓 线也不会太黑糊。14试验 2 图像轮廓线提取技术在 MATLAB 的图像处理中，导数算子具有突出灰度变化的作用，对图像运 用导数算子， 灰度变化较大的点处算得的值比较高，因此可将这些导数值作为相 应点的边界强度，通过设置门限的方法，提取边界点集。 一阶导数与是最简单的导数算子， 它们分别求出了灰度在 x 和 y 方向上的变 化率，而方向α 上的灰度变化率可以用相应公式进行计算;对于数字图像，应该 采用差分运算代替求导，差分公式参考相关教材。 函数 f 在某点的方向导数取得最大值的方向是， 方向导数的最大值是称为梯 度模。利用梯度模算子来检测边缘是一种很好的方法，它不仅具有位移不变性， 还具有各向同性。为了运算简便，实际中采用梯度模的近似形式。另外，还有一 些常用的算子，如 Roberts 算子和 Sobel 算子。 由于 Sobel 算子是滤波算子的形 式，用于提取边缘。也可以利用快速卷积函数，简单有效，因此应用很广泛。 拉普拉斯高斯（loG）算法是一种二阶边缘检测方法。它通过寻找图像灰度 值中二阶微分中的过零点（Zero Crossing）来检测边缘点。其原理为，灰度级变 形成的边缘经过微风算子形成一个单峰函数，峰值位置对应边缘点；对单峰函数 进行微分，则峰值处的微分值为 0，峰值两侧符号相反，而原先的极值点对英语 二阶微分中的过零点，通过检测过零点即可将图像的边缘提取出来。 MATLAB 的图像处理工具箱中提供的 edge 函数可以实现检测边缘的功能， 其语法格式如下： 1、采用 Sobel 算子进行边缘检测 BW = edge(I,'sobel') 2、指定算子方向的 Sobel 算子边缘检测 BW = edge(I,'sobel',direction) 其中： direction=‘horizontal’，为水平方向； direction=‘vertical’，为垂直方向； direction=‘both’，为水平和垂直两个方向。 3、Roberts 算子 BW = edge(I,'roberts') 4、拉普拉斯高斯算子 BW = edge(I,'log') 七、练习 1、任意选取一个灰度图像和彩色图像，对算法中若干关键语句中的参数进行调 整，得出不同的实验结果，并对这些结果进行分析。 2、根据自己所学知识，提出自己的轮廓线提取方法，与简单阈值法进行比较分 析。 3、练习 matlab 自带算子的检测结果。 4、写出图像轮廓线提取的实验报告。15试验 3RGB 向量空间中的图像分割技术实验 3RGB 向量空间中的图像分割技术一、实验目的与要求： 1、能熟练应用 matlab 去分析问题、解决问题； 2、熟悉对 matlab 的图像处理的功能，掌握基本的图像处理的若干命令； 3、在应用 matlab 进行图像处理方面具备一定的编程能力。 4、掌握 figure，imread，image，colormap，imshow，imwrite，subplot， title，rgb2gray，imfinfo 等语句的基本使用方法。 5、掌握图像分割技术中的简单方法并上机实现。 6、进一步提高对复杂图像处理的能力。 二、问题描述 “图像分割”简单地说，就是要把图像中有意义的区域与背景分离开，依据 的原则是同一区域具有“同质性”“区域”是图像中相邻的具有类似性质的点组 。 成的集合。 图像分割是在图像预处理的基础上对信息进行组织与加工，它是在对图像自 动识别与理解之前必不可少的步骤， 其目标就是将图像划分为具有区域强相关性 的各组成部分。 三、问题分析 “图像分割”模拟提取图像中景物的信息。因此，可以通过图像中物体所具 备的特征（如大小、形状、颜色、纹理等）来指导分割，并利用这些特征之间的 相互信息，进行更复杂的分割。 常见的分割方法有：阈值法、类间最大距离法、区域增长法等。本实验重点 实现前面两种方法。 四、背景知识介绍 重点介绍两种简单的分割方法，简单阈值法和类间最大距离法。 1、简单阈值法 在理想状态下， 背景与对象之间的灰度值相差很大，且同一个对象具有基本 相同的灰度值。体现在图像的灰度直方图上，就是直方图呈明显的双峰分布，两16试验 3RGB 向量空间中的图像分割技术类灰度级之间无交叠。 因此， 可以认为处于直方图谷底的灰度值就是目标对象和背景之间的分界线。称 该灰度为“阈值” ，选取了阈值以后，根据实际情况，把灰度高于该阈值的像素 点置为白色， 低于该阈值的像素点置为黑色。这样处理之后就可以把目标图像和 背景分开了。 设 f(x, y)表示原图像，g(x, y)表示分割后的图像，T 为选定的灰度阈值，分 割算法表示为：?255 , g ( x, y ) ? ? ? 0, f ( x, y ) ? T f ( x, y ) ? T从表达式可以看出，算法的关键在于阈值 T 的确定。在理想状态下，我们 可以取直方图中处于谷底的点作为阈值；或者通过多次调试，确定效果最优的灰 度值作为阈值。 2、类间最大距离法 前面介绍的简单阈值法有一个最大的不足，就是阈值的选取不够智能化。于 是我们想寻求一种智能的分割方法。类间最大距离法就是其中一种方法。 显然，一幅灰度图像的目标区域跟背景区域有较大的灰度区别，才能使我们 分辨出目标和背景。 因此可以设计一种算法，以不同区域之间的距离达到最大为 目标，用运筹学中的规划方法，得到实现目标的最优阈值。现在，关键的问题是 如何定义区域之间的“距离”。可以采用不同的方法来衡量距离，这里采取一种 相对简单的方法，就是把两类区域像素灰度均值的方差作为度量的距离。 ? 灰度均值的计算： （1）ui = Σ f(x,y) /N;其中：f(x,y)表示二维图像(x,y)处像素点的灰度值，N 表示像素点的个数。 ? 两区域间相对距离的计算：给定一个初始阈值 Th0，将图像分为 C1 和 C2 两区域，对应区域的灰度均值 分别为 U1、 U2。则两区域的相对距离为： S = (U2 C Th0)×(Th0 C U1)/(U2 C U1)^2 ? 选择最佳的阈值 （2）17试验 3RGB 向量空间中的图像分割技术给定阈值变化的区间[a,b],用优化函数 fmin（）获得区间范围内的最佳阈值 Th，使得图像按照该阈值 Th 分为 C1 和 C2 两类后，满足:S Th ? max { S T }T ?[ a ,b ]五、实验过程 1、简单阈值分割实验 实验流程图：源代码（.m 文件） ： function I = smpthrshd(A , T) %简单阀值法的源代码，为了提高函数的可重用型，把原图像作为参数 I = A; %I 是目标图像，与原图像大小相同 [n , m] = size(I); %读出原图像的大小 for i = 1 : n for j = 1 : m %两重循环，实现扫描矩阵的每一个像素点18试验 3RGB 向量空间中的图像分割技术if A(i , j) &= T I(i , j) = 0; elseif A(i , j) & T I(i , j ) = 255;%以下代码实现： %当像素点的灰度小于阈值则置为黑色 %当像素点的灰度大于阈值则置为白色 imshow(I); %显示分割后的图像 title('经过简单阀值分割后的图像'); 2、类间最大距离法分割实验 算法流程图如下：算法源代码（.m 文件）如下： %――――――子函数 1，计算灰度均值――――――% function [U1, U2] = getaverage(I, Th) %计算两个区域的灰度均值 %阈值 Th 把图像分成两个区域：灰度高于 Th 的区域和灰度低于 Th 的区域 [n,m] = size(I); sum1 = 0; %设置灰度和的初始值 sum2 = 0; num1 = 0; %设置每个区域像素点个数的初值 num2 = 0; for i = 1 : n for j = 1 : m if I(i , j) &= Th %两重循环，实现扫描整个矩阵 %如果灰度小于或等于阈值，则属于区域 119试验 3RGB 向量空间中的图像分割技术num1 = num1 + 1; %区域 1 的像素个数加 1 sum1 = sum1 + I(i , j);%区域 1 灰度之和加上当前像素点的灰度 else %否则，该像素点属于区域 2 sum2 = sum2 + I(i , j);%区域 2 灰度之和加上当前像素点的灰度 num2 = n * m - num1; %区域 2 的像素个数等于所有像素数减区域 1 的像素数 if num1 ~= 0 %如果区域 1 不为空 U1 = sum1 / num1; %区域 1 的灰度均值等于灰度和除以区域 1 的像素个数 else U1 = 0; %如果区域 1 为空，则设置其灰度均值为 0 if num2 ~= 0 %如果区域 2 不为空 U2 = sum2 / num2; %区域 2 的灰度均值等于灰度和除以区域 2 的像素个数 else U2 = 0; %如果区域 2 为空，则设置其灰度均值为 0 %――――――子函数 2 获取类间距离――――――% function S = get_S(Th, I) % 获得距离的函数 [U1,U2] = getaverage(I,Th); %调用函数 getaverage，返回两区域的灰度均值 if (U1 - U2) ~= 0 %如果两个区域的灰度均值不相等 S = - (U1 - Th) * (Th - U2)/(U1 - U2)^2; %用式 2.2 计算距离 else S = 0; %如果灰度均值相等，则距离为 0 %――――――主函数，用类间最大距离法进行分割――――――% function I = maxsp(A) %利用类间最大距离法实现分割，其中用到两个函数:getaverage,get_S A=double(A); %因为要对矩阵元素进行加减操作，先转成双精度类型 %下面用优化函数获取使距离达到最大的那个阈值 %[Th,option] = fmin('get_S',0,255,[1,1.e-12],A); %for matlab 6.5.1 Th = fminbnd(@(Th) get_S(Th,A),0,255); %for matlab 7.0.1 %得到最优的阈值后，下面的操作与简单阈值法相同 I = A; [n, m] = size(I); for i = 1 : n20试验 3RGB 向量空间中的图像分割技术for j = 1 : m if A(i , j) &= Th I(i , j) = 0; elseif A(i , j) & Th I(i , j) = 255; imshow(I); title('经过类间最大距离法分割后的图像'); 六、结论与应用 1、简单阈值分割实验 下面以肝癌的 CT 图片为例进行实验。图 1 肝癌的 CT 图片（左侧两箭头指向的区域为癌块）图 2 原图的灰度直方图灰度直方图的分析：从图２可以看出直方图出现了两个波峰，对照原图， 我可以发现图１中根据像素基本上可以分为两个区域。 这两个区域分别具有不同 的亮度，恰好对应着直方图中的两个波峰。因此，理论上我们可以取 100 到 230 之间的任何一个灰度作为阈值， 但为了分割后图像的效果较为美观，阈值最好取 在 100 到 170 之间。图 3 取 T=100 后的分割效果21试验 3RGB 向量空间中的图像分割技术简单阈值法只考虑直方图的灰度信息，没有利用图像的其他信息，如图像 中包含物体的大小， 光源的位置等。在图像直方图为单峰或者图像较复杂的情况 下无效。 而且这种方法的一个最大的不足就是对阈值的选取不够智能化，必须根 据直方图，通过人工的判断，并经过多次尝试才能得到好的分割效果。虽然对于 简单的图片， 这种算法分割的效果比较能让人满意，但如果有大量的图片需要处 理或者人工无法干预的情况下该怎么办呢？有没有一种更智能一些的方法呢？ 我们可以先不要求有那么高的处理效果，但要求智能，以减少工作量。 2、类间最大距离法分割实验 把图 1 作为原图，采用类间最大距离法进行分割，得到最优的阈值为 133.5159， 分割效果图如下:图 4 阈值为 89.5448 的分割效果图类间最大距离法实现了智能选取阈值的功能，分割效果不错。对于一些特殊 的图像，相信根据实际情况给出新的距离的定义，也可以得到很好的分割效果。 但这种方法依然没有解决简单阈值法中提到的一个问题， 就是同样只考虑直方图 的灰度信息， 没有考虑到图像的其他有用信息， 如图像的大小和光源的位置等等！ 如果能把这些有用的信息考虑进去，对图像的分割将可以达到更好的效果。 七、练习 1、任意选取一个灰度图像，对算法中若干关键语句中的参数进行调整，得出 不同的实验结果，并对这些结果进行分析。 2、根据自己所学知识，提出自己的阈值选取方法，与类间最大聚类法的分割 效果进行比较分析。 3、写出图像分割的实验报告。22试验 4图像的伪装技术试验 4 图像的伪装技术一、实验目的与要求 让学生了解图像伪装的意义以及相应的算法在实际生活中的价值。 通过实验 能够掌握基本的图像加密和伪装算法， 并能够对具有保密价值的图形图像作出合 理巧妙的伪装处理。 二、问题描述 对一幅具有保密价值的数字图像进行伪装与加密。 三、问题分析 在已知各种图像格式的前提下， 利用现有的线性代数及高等数学的相关知识 实现对一幅具体图片的伪装转换，利用各种线性或非线性变换产生一串钥匙链， 并对此链进行二次加密， 从而实现对图像的保密。 另一方面还可以通过矩阵变换， 将原图变换为其他不相干的图片的方法实现对原图的伪装。 四、背景知识介绍 图像的伪装技术是随着计算机的发展和人们对各种特殊图片的保密需求而 发展起来的一门新型技术。数字图像具有信息量大、信息表达直观的特点，人们 易于从图像图形来感知世界。 随着计算机性能和网络带宽的迅速提高，数字图像 及视频不再难以存储、传输和处理。数字图像的自相关性会造成大量的冗余，而 图像的自相关性和人眼的分辨能力决定了即使图像有所变动， 人的视觉系统仍然 可以提取出数字图像的特征， 这就是以数字图像为载体并在其中实现信息隐藏的 前提和基础。 保护图像信息安全最容易想到的方法是直接采用加密技术。但是现有的保密 系统大多是为保护文本信息而设计的，直接用它们来保护图像信息并不合适，因 为图像的数据量一般比文本的数据量大许多， 如果直接采用传统的保密系统则加 解密所需的时间将会大大增加。此外，在解密文本数据(密文)时得到的文本必须 和原来的文本(明文)相同，但是图像则不需要满足这个条件，因为图像最终是由 人的肉眼来观察的，只要解密出的图像的失真在主观视觉上可以接受的范围内， 即可认为是达到了要求。 所以实现数字图像的安全性应该不同于文本数据，对数 字图像信息可以设计出效率更高的安全保密方法。 随着多媒体技术和计算机网络的飞速发展，数字产品的版权和信息安全问题 日益倍受关注。 数字水印技术和信息伪装技术正是适应这一要求发展起来的。前 者是为了保护数字产品的版权， 而要对数字产品加载所有者的水印信息，以便在 产品的版权产生纠纷时作为证据。而后者是指把不便于公开的信息，通过伪装算 法加载到可以公开的信息载体中去，在信息的接收端利用解伪算法，恢复出被伪 装的信息。 前者是为了保护载体数字产品的版权问题，而后者则是为了被伪装信 息的安全与保密。二者似乎没有本质的差别，而只不过是使用的目的不同而已。 但是，从技术要求来看，二者又有很大的差别，一般的，数字水印的技术要求较 高，要求所加载的水印不仅要具有安全性，还要具有很强的稳健性.而后者的要23试验 4图像的伪装技术求则更加偏重于安全与保密性. 目前，数字水印技术已经成为人们研究的重点，已经出笼的水印算法很多， 基本上可以分为空间域、频域变换域等几种。但是对信息伪装技术的研究较少。 数字水印技术可以用于信息伪装，但所伪装的信息较少，伪装之后的冗余太大， 因而有必要专门研究信息伪装技术，以供保护信息的安全和保密使用。 自从信息伪装和数字水印的概念提出以来， 信息伪装技术和数字水印技术都 得到了迅速发展。 但用于保密通信的信息伪装技术远没有数字水印技术发展得迅 速。 由于图像在信息的表达中扮演了重要的角色，人们也将会越来越多地使用图 像.因此，如何保证图像信息的安全已成为人们关注的重要问题。 近年来，人们已提出了几种保护图像信息的密码系统。它们都是将要保密的 图像加密变成一幅经过编码的密图(Cipher image)如果密图在传送过程中被非法 拦截， 当拦截者试图查看密图内容时，他只能看到一幅由一堆乱码所构成的毫无 意义的图像。但另一方面，这也容易引起拦截者的怀疑，从而激发起了拦截者破 解密图信息的兴趣。 为了弥补这个缺点，我们着眼于研究基于信息伪装的保护图 像信息的方法。 其基木思想是经过编码的密图仍是有意义的，因此不易引起拦截 者的注意，从而大大减少了遭受拦截者试图破译的可能性。图像经过伪装后，多 了一层保护，仿佛是自然界中生物的保护色，增加了图像的保密性和安全性. 本实验提出了一种频域的信息伪装技术，该技术隐藏效果好，人眼系统不能 发觉隐藏了秘密图像的载体图像与原载体图像之间的差别；隐藏的内容是安全 的；可以有效地抵抗 JEPG 等图像压缩编码的破坏，稳健性很好。此外，该算法 还可以用于数字水印， 而且可以将水印图像由二值图像扩展到灰度图像；可以在 同一幅载体图像中隐藏多幅秘密图像。 五、实验过程 1、基于信息隐藏的图像加密 设 ci（carrier image）是 m×n 大小的载体图像，si(secret image)为想要隐藏 的 p×q 大小的秘密图像，p=m/2，q=n/2。系数 a 将对隐藏后的载体图像 ci 的质 量以及由其恢复的新秘密图像的质量产生重大的影响，a 越大，恢复的新秘密图 像的质量越好；相反，a 越小，图像 ci 的质量越好。一般的，可以将 a 的值取为 0.01-0.1 之间。通过本算法，我们要把图像 si 隐藏到图像 ci 中。算法分以下几步 进行： 1. 0.1=&a 2. 生成 8×8 的离散余弦变换矩阵 c 3. 读入载体图像 ci 和秘密图像 si 4. ci 的第 1 个平面=&ci1，si 的第 1 个平面=&si1 5. 对 si1 的每一个行向量进行低通滤波得 si2，间隔采样后得 si3 6. 对 si3 的每一个列向量进行低通滤波得 si4，间隔采样后得 si5 7. 将 ci1 分解为（m/8）×（n/8）个 8×8 大小的方块 bci1 8. 对每个 bci1 进行 DCT 变换，即 bci2＝c'×bci1×c 9. 将 si5 分解为（m/8）×（n/8）个 2×2 大小的方块 bsi1 10. 对于每一个 bsi1 及其对应的 bci2 按下图所示进行数据隐藏，即 si＝a×ti，其 中 i＝1，2，3，424试验 4图像的伪装技术s1 s2 t1 t2 t3 t4 bsi1 s3 s4bci211. 再对 bci2 进行逆 DCT 变换，即 bci3=c×bci2×c' 12. 将各个 bci3 合并成图像的一个面 13. 对载体图像和秘密图像的其它两个面进行同样的处理 14. 再将三个面合并为一个整图，即得隐藏有秘密图像得载体图像? 1 ? 2 ? 2 ? cos 1 ? ? 2n n? ? ? n ?1 ? ?cos 2n ? 1 2 3 cos ? 2n cos 3(n ? 1) ? 2n ? 1 ? 2 ? 2n ? 1 ? ? cos ? ? 2n ? ? (2n ? 1)( n ? 1) ? ? cos ?? 2n ? ?c?源代码（.m 文件） clc clear a=0.1; n=8; %生成 8×8 的离散余弦变换矩阵 for i=1:n-1 for j=1:2:2*n-1 c(i+1,(j+1)/2)=i*j*pi/(2*n); end end c(1,:)=0; c=cos(c); c(1,:)=sqrt(1/2); c=sqrt(2/n)*c; ci=imread('ci.jpg'); %读入载体图像和秘密图像 figure(1) imshow(ci) %显示载体图像 ci=double(ci); si=imread('si.jpg');25试验 4图像的伪装技术figure(2) imshow(si) %显示秘密图像 si=double(si); for k=1:3 %将图像分成三个面来处理 si1=si(:,:,k); ci1=ci(:,:,k); [p,q]=size(si1); for j=1:q-1 %对 si1 进行行低滤通波 si2(:,j)=si1(:,j)+si1(:,j+1); end si2(:,q)=si1(:,q)+si1(:,1); for j=1:2:q %间隔采样 si3(:,(j+1)/2)=si2(:,j)/2; end for i=1:p-1 %对 si3 进行列低滤通波 si4(i,:)=si3(i,:)+si3(i+1,:); end si4(p,:)=si3(p,:)+si3(1,:); for i=1:2:p %间隔采样 si5((i+1)/2,:)=si4(i,:)/2; end [m,n]=size(ci1); for i=1:m/8 %对载体图像和秘密图像进行分块 for j=1:n/8 bci1=ci1((i-1)*8+1:(i-1)*8+8,(j-1)*8+1:(j-1)*8+8); bci2=c'*bci1*c; %对载体图像的每一块进行 DCT 转换 bsi1=si5((i-1)*2+1:(i-1)*2+2,(j-1)*2+1:(j-1)*2+2); bci2(2,6)=a*bsi1(1,1); %数据隐藏 bci2(3,7)=a*bsi1(1,2); bci2(5,3)=a*bsi1(2,1); bci2(6,4)=a*bsi1(2,2); bci3=c*bci2*c'; %对载体图像的每一块进行逆 DCT 转换 %将载体图像的每一块合并为一个整图，即为隐藏有秘密图像的新图像 ci1((i-1)*8+1:(i-1)*8+8,(j-1)*8+1:(j-1)*8+8)=bci3; end end si6(:,:,k)=si5; %si6 是经过压缩后的秘密图像 tci(:,:,k)=ci1; %将图像的三个面合并 end tci=uint8(tci); figure(3) imshow(tci) %显示已经隐藏了秘密图像的载体图像 imwrite(tci,'tci.jpg','Quality',100)26试验 4图像的伪装技术信息提取： 设图像 tci 为已经隐藏了秘密图像的载体图像。现将所隐藏的秘密图像从 tc 工中提取出来。其过程为上述隐藏算法的逆运算： 1. 0.1=&a 2. 生成 8×8 的离散余弦变换矩阵 c 3. 读入载体图像 tci 4. 将 tci 的第一个面分解为（m/8）×（n/8）个 8×8 大小的方块 bci1 5. 对每个 btci1 进行 DCT 变换，即 btci2＝c'×btci1×c 6. 对每一个 btci2，按下式得到 btsi1，即 si＝（1/a）×ti，其中 si 和 ti 如上图所 示 7. 将各个 btsi1 合并成图像的一个面 wsi1 8. 对 wsi1 进行逆 DWT 变换的 wsi3 9. 对载体图像 tci 的其它两个面进行同样的处理得到秘密图像的其它两个面 10. 将这三个面合并即得隐藏在载体图像中的秘密图像 源代码（.m 文件） %提取图像 clc clear a=0.1; n=8; %生成 8×8 的离散余弦变换矩阵 for i=1:n-1 for j=1:2:2*n-1 c(i+1,(j+1)/2)=i*j*pi/(2*n); end end c(1,:)=0; c=cos(c); c(1,:)=sqrt(1/2); c=sqrt(2/n)*c; tci=double(imread('tci.jpg')); for k=1:3 %将图像分成三个面来处理 [m,n]=size(tci(:,:,k)); for i=1:m/8 for j=1:n/8 btci1=double(tci((i-1)*8+1:(i-1)*8+8,(j-1)*8+1:(j-1)*8+8,k)); btci2=c'*btci1*c; %对已经隐藏了秘密图像的载体图像的每一块进行 DCT 转换 btsi1(1,1)=btci2(2,6)/a; %提取数据 btsi1(1,2)=btci2(3,7)/a; btsi1(2,1)=btci2(5,3)/a; btsi1(2,2)=btci2(6,4)/a; %将秘密图像的每一块合并为一个整图 wsi1((i-1)*2+1:(i-1)*2+2,(j-1)*2+1:(j-1)*2+2,k)=btsi1; end27试验 4图像的伪装技术end [p,q]=size(wsi1(:,:,k)); for i=1:2*p %对提取到的秘密图像进行逆 DWT 变换 wsi2(i,:,k)=wsi1(round(i/2),:,k); end for j=1:2*q wsi3(:,j,k)=wsi2(:,round(j/2),k); end wsi(:,:,k)=wsi3(:,:,k); end wsi=uint8(wsi); imshow(wsi) %显示提取到的秘密图像图1图2图3图4图5图628试验 4图像的伪装技术2、基于置乱算法的图像加密技术 图像的置乱算法的核心是将图像中像素的位置互相置换。步骤如下： (1)对原图像取一个固定模框，模框中像素位置排列如下所示； (2)建立一个宽度与原图像模框不同的置乱模框, 在置乱模框中把原图模框中的 像素位置信息按横向次序填入； (3)将置乱模框中的像素位置信息按纵向次序读取，填回原图像中就得到信息次 序的混乱。1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 置乱模框1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 原图像1 4 7 10 13 16 2 5 8 11 14 3 6 9 12 15 结果模框其中置乱模框的宽度是密匙的特征参数之一。其长度则随原图像模框中像素 的多少而变化。置乱模框除上面给出的方法之外，还可以选取其它的形状，如菱 形，星形，三角形等。图像经过多次上述的“置乱”(取不同的特征参数)后，信 息变得杂乱无章，不知道“密匙”的人几乎没有办法恢复图像的原来面貌。 源代码（.m 文件） clc clear si=double(imread('si.jpg')); %si=[1 2 3 4;5 6 7 8;9 10 11 12;13 width=30; for k=1:3 si1=si(:,:,k); [p,q]=size(si1); s=1; i2=1; j2=1; for i=1:p for j=1:q si2(i2,j2)=si1(i,j); if j2==width j2=1; i2=i2+1; else j2=j2+1; end s=s+1; end2914 15 16];试验 4图像的伪装技术end [p2,q2]=size(si2); if p*q~=p2*q2 si2(p2,q2-(p2*q2-p*q)+1:q2)=-1; end s=1; i2=1; j2=1; for j=1:q2 for i=1:p2 if si2(i,j)~=-1 si3(i2,j2)=si2(i,j); if j2==q j2=1; i2=i2+1; else j2=j2+1; end s=s+1; end end end nsi(:,:,k)=si3; end imshow(uint8(nsi));图 7 原图 按幻方的图像置乱变换 以自然数 1，2，?，n2 为元素的 n 阶矩阵图8置乱后的效果? a11 ?a A ? ? 21 ? ? ? ?a n1na12 a 22 ? an 2? a1n ? ? a2n ? ? ? ? ? ? ? a nn ?i ?1 ( j ?1, 2,?, n )满足j ?1 ( i ?1, 2,?, n )? a ij ?? a ij ? ? a ij ?i? jnni ? j ? n ?1?anij?c30试验 4图像的伪装技术其中 c ?n(n 2 ? 1) 2则称 A 为标准幻方。 假定数字图像相应于 n 阶数字矩阵 B。 对取定的 n 阶幻方 A， B 与 A 按行 将 列作一一对应。 我们把 A 中的元素 1 移至 2 的位置， 将元素 2 移至 3 的位置， ?。 一般说来，对任何 m ? ? ,2,?, n 2 ? 1?，它从其在 A 中的位置移至 m＋1 在 A 中所 1 在的位置。若 m ? n 2 ，则将 n 2 移至 1 在 A 中所在的位置。经过这样的置换之后， 矩阵 A 转换为矩阵 A1 。记 A1 ? EA 。对 A1 来说，可以重复上述置换得矩阵A2 ? EA1 ，继而 A3 ? EA2 ， A4 ? EA3 等等，这便是一系列的置乱变换，经过 n2步， An 2 ? A 。一般说，A1，A2，?， An 2 ?1 并不满足式（1） ，它们不是幻方。 对数字图像矩阵 B，注意 B 与 A 元素之间的对应关系，随 A 转换为 A1 而把 B 中对应象素信息（灰度，RGB）做相应的移置，产生数字图像矩阵 B1，记为EB ? B1 ,一般地，有 E m B ? Bm 。例如，我们考虑一个由 4×4 象素组成的图像 B，它可以被看作是一个 4 阶 方阵?b11 b12 ?b b B ? ? 21 22 ?b31 b32 ? ?b41 b42b13 b23 b33 b43b14 ? b24 ? ? b34 ? ? b44 ?在如下 4 阶幻方 A16 9 2 7 3 13 8 1 6 125 11 10 4 14 15置乱作用下，得到新的 4 阶方阵 B1?b43 ?b B1 ? EB ? ? 41 ?b24 ? ?b13b44 b23 b33 b14b12 b31 b21 b42b34 ? b32 ? ? b22 ? ? b11 ?源代码（.m 文件） clc31试验 4图像的伪装技术clear si1=double(imread('smallsi2.jpg')); imshow(uint8(si1)); [p,q]=size(si1(:,:,1)); if p~=q error('图像的维数必须是 M×N×3，并且 M=N'); end magic=magic(p); for n=1:1 for k=1:3 for i1=1:p for j1=1:q for i2=1:p for j2=1:q temp=magic(i1,j1)-1; if temp==0 temp=p*q; end if magic(i2,j2)==temp si2(i1,j1,k)=si1(i2,j2,k); end end if j2~= end end end end end si1=si2; end imshow(uint8(si1)); 六、结论及应用 通过实验我们可以了解到， 对数字图像的伪装与加密其实就是数学上的矩阵 变换。我们可以根据实际要求，针对不同图像所要求的保密等级，设计不同的加 密深度。事实上，由于矩阵变换的任意性和复杂性，我们可以随意定义一个变换 规则，基于这个规则，可以实现任意复杂的加密算法。尽管从理论上讲任何加密 得法都有被破解的可能性， 但由于算法的任意性，要实现对一个任意抽象的算法 的还原， 可能会耗费大量的工作，而这种工作在时间和精力上的开销对于破解一 幅图片可能已经没有任何意义。 而且我们还可以采用图像的伪装技术，使得对图 像的还原出现不同的还原路径， 即使别人从截获的数据恢复出若干个图像，由于32试验 4图像的伪装技术无从判断原始图像的真实性， 也会使得这种破解过程变得没有实用价值。从这个 角度来看， 对图像的加密和伪装二者的结合，当算法具有较高复杂度和任意性的 前提下， 如果不事先知道对图像变换的顺序和所用的变换函数类，要实现对原始 图像的还原，基本是不可能的。 七、练习 1、找一幅需要保密的图片，通过定义一个变换规则，对其进行加密； 2、将一幅需要保密的图片，通过矩阵变换变为众人司空见惯的另一幅图片， 以实现对图像的伪装。33试验 5用矩阵变换实现对图像形状及颜色畸变的校正实验 5用矩阵变换实现对图像形状及颜色畸变的校正一、实验目的与要求 让学生了解数字图像的数学表达及相关概念，通过实验让学生加深对数学 在相关学科的应用价值的认识， 培养学生的实际操作能力，并引导他们建立基础 学科在处理具体问题时方法上联系。 二、问题描述 对于在颜色或形状上发生畸变的图像，通过数学的方法实现校正。 三、问题分析 先由教师讲授数字图像的基本概念（包括图像的数学化、采样、量化、灰 度、各种数学图像的文件格式、表色系、颜色映像等） ，再通过具体的实例给学 生示范对于在颜色或形状上发生畸变的图像如何通过数学的方法实现校正的过 程。 最后让学生动手完成对某些特殊畸变的图像的校正，写出数学原理和实验报 告。 四、背景知识介绍 1、数字图像的数值描述及分类 图像是对客观存在物体的一种相似性的生动模仿与描述， 是物体的一种不完 全的不精确的描述。数字图像是用一个数字阵列来表示的图像。数字阵列中的每 个数字，表示数字图像的一个最小单位，称为像素。采样是将空域上或时域上连 续的图像变换成离散采样点（像素）集合的一种操作。 对一幅图像采样后， 若每行像素为 M 个，每列像素为 N 个，则图像大小为 M ? N 个像素。例如，一幅 640 ? 480 的图像，就表示这幅连续图像在长、宽方向上分 别分成 640 个和 480 个像素。显然，想要得到更加清晰的图像质量，就要提高图 像的采样像素点数，即使用更多的像素点来表示该图像。 客观世界是三维的，从客观场景中所拍摄到的图像是二维信息。因此，一幅 图像可以定义为一个二维函数 f(x,y),其中 x,y 是空间坐标。对任何一对空间坐标 (x,y)上的幅值 f(x,y)，成为表示图像在该点上的强度或灰度，或简称为像素值。 因为矩阵是二维结构的数据，同时量化值取整数，因此，一幅数字图像可以用一 个整数矩阵来表示。矩阵的元素位置(i,j)，就对应于数字图像上的一个像素点的 位置。矩阵元素的值 f(i,j)就是对应像素点上的像素值。 值得注意的是矩阵中元素 f(i,j)的坐标含义是 i 为行坐标，j 是列坐标。而像 素 f(x,y )的坐标含义一般指直角坐标系中的坐标，两者的差异如下图：试验 5用矩阵变换实现对图像形状及颜色畸变的校正0 矩阵元素 f (i,j) 行坐标(i)列坐标(j)纵坐标(y)像素 f(x,y) 0 图 1.1 矩阵坐标系与直角坐标系 横坐标(x)对应于不同的场景内容，数字图像可以大致分为二值图像，灰度图像，彩色 图像三类。 1）二值图像 它是指每个像素不是黑就是白，其灰度值没有中间过度的图像。二值图像对 画面的细节信息比较粗略， 适合于文字信息图像的描述。 它的矩阵取值非常简单， 即 f(i,j)=0(黑),或 f(i,j)=1(白)，除此之外没有其他的取值。当然，0 和 1 表示黑或 白都只是人定义的，可以人为地反过来定义。这种图像具有数据量小的优点。 2）灰度图像 它是指每个像素的信息由一个量化后的灰度级来描述的数字图像， 灰度图像 中不包含彩色信息。 标准灰度图像中每个像素的灰度有一个字节表示，灰度级数 位 256 级，每个像素可以是 0~255（从纯黑到纯白）之间的任何一个值。值越接 近 0 就越黑，越接近 255 就越白。 3）彩色图像 常用的图像彩色模式有 RGB 模式、CMYK 模式和 HIS 模式，一般情况下只 使用 RGB 模式。它是根据三基色成像原理来实现对自然界中的色彩描述的。这 一原理认为，自然界中的所有颜色都可以由红，绿，蓝（R,G,B）三基色组合而 成。如果三种基色的灰度分别用一个字节（8bit）表示，则三基色之间不同的灰 度组合可以形成不同的颜色。 2、数字图像质量决定因素 数字图像的效果与以下几个评价参数有关。 1）图像分辨率 即采样所获得的图像总像素的多少。 2）采样密度 即在图像上单位长度所包含的采样点数。 采样密度的倒数是 像素间距。 3）采样频率 即一秒钟采样的次数。它反映了采样点之间的间隔大小，采 样频率越高，丢失的信息越少，图像的质量越好。 4）扫描分辨率 表示一台扫描仪输入图像的细微程度， 指每英寸扫描所得的 点，单位是 DPI（DotPerInch） 。数字越大，表示被扫描的图像转化为数字 化图像越逼真，扫描仪质量也越好。 3、彩色空间 1） 三基色原理 近代的三色学说研究认为,人眼的视网膜中存在着三种锥体细胞,它们包含不 同的色素,对光的吸收和反射特性不同,对于不同的光就有不同的颜色感觉.研究 发现,第一种锥体细胞专门感受红光,第二和第三种锥体细胞则分别感受绿光和蓝35试验 5用矩阵变换实现对图像形状及颜色畸变的校正光.它们三者共同作用,使人们产生了不同的颜色感觉. 这三种色光以不同比例混 合,几乎可以得到自然界中的一切色光,混合色域最大;而且这三种色光具有独立 性,其中一种原色不能由另外的原色光混合而成,由此,称红,绿,蓝为色光三原色. 为 了 统 一 认 识 ,1931 年 国 际 照 明 委 员 会 (CIE) 规 定 了 三 原 色 的 波 长 nmR0.700=λ,nmG1.546=λ,nmB8.435=λ。 2） 彩色的基本特征 色调(hue) 色调又称为色相,是当人眼看到一种或多种波长的光时所产生 的彩色感觉,它反映颜色的种类,是决定颜色的基本特性.色调用红,橙,黄,绿,青,蓝, 靛,紫等术语来刻画. 不透明物体的色调是指该物体在日光的照射下,所反射的各 光谱成分作用于人眼的综合效果;透明物体的色调则是透过该物体的光谱综合作 用的效果。 饱和度(saturation) 饱和度是指颜色的纯度,即色彩含有某种单色光的纯净 程度,它可用来区别颜色的深浅程度.对于同一色调的彩色光,饱和度越深颜色越 鲜明或说越纯,例如鲜红色饱和度高,而粉红色的饱和度低.完全饱和的颜色是指 没有渗入白光所呈现的颜色,例如仅由单一波长组成的光谱色就是完全饱和的颜 色。 亮度(brightness) 亮度是视觉系统对可见物体辐射或者发光多少的感知属 性.亮度是光作用于人眼时所引起的明亮程度的感觉,它与被观察物体的发光强度 有关.由于其强度的不同,看起来可能会亮一些或暗一些.对于同一物体,照射光越 强,反射光也越强,感觉越亮;对于不同的物体在相同照射情况下,反射越强者看起 来越亮。 通常把色调和饱和度通称为色度.亮度是用来表示某彩色光的明亮程度,而 色度则表示颜色的类别与深浅程度。 五、实验过程 1、图像畸变介绍 从数字图像处理的观点来考察畸变校正, 实际上是一个图像恢复的过程, 是对一幅退化了的图像进行恢复。在图像处理中，图像质量的改善和校正技术， 也就是图像复原，当初是在处理从人造卫星发送回来的劣质图像的过程中发展、 完善的。目前，图像畸变校正的应用领域越来越广，几乎所有涉及应用扫描和成 像的领域都需要畸变校正。 图像在生成和传送的过程中， 很可能会产生畸变， 如： 偏色、模糊、几何失真、几何倾斜等等。前几种失真主要是体现在显示器上，而 后一种失真则多与图像集角度有关。不正确的显影，打印、扫描，抓拍受反射光 线的影响等方式，都会使图像产生偏色现像。模糊、几何畸变主要是在仪器采集 图片过程中产生，大多是因机器故障或操作不当影响导致，如在医学成像方面。 而几何空间失真广泛存在于各种实际工程应用中，尤其是在遥感、遥测等领域。 2、图像畸变校正过程所用到的重要工具 灰度直方图是关于灰度级分布的函数，是对图象中灰度级分布的统计。灰度 直方图是将数字图象中的所有像素，按照灰度值的大小，统计其所出现的频度。 通常，灰度直方图的横坐标表示灰度值，纵坐标为想像素个数。直方图上的一个 点的含义是， 图像存在的等于某个灰度值的像素个数的多少。这样通过灰度直方 图就可以对图像的某些整体效果进行描述。从数学上讲，图像的灰度直方图是图 像各灰度值统计特征与图像灰度值出现的频率。从图形上来讲，它是一个一维曲36试验 5用矩阵变换实现对图像形状及颜色畸变的校正线，表征了图像的最基本的统计特征。 作为表征图像特征的信息而在图像处理中起着重要的作用。由于直方图反映 了图像的灰度分布状况， 所以从对图像的观察与分析， 到对图像处理结果的评价， 灰度直方图都可以说是最简单、最有效的工具。 3、图像颜色畸变校正介绍 图像颜色畸变现象可以是由摄像器材导致， 也可以是由于真实环境本身就偏 色导致，还有的是由于图像放置过久氧化、老化导致。无论其产生的原因如何， 其校正方法都是类似的。 图像颜色畸变校正在社会生活、工作中应用十分广泛。小到家庭生活图像处 理，大到医学成像应用、罪犯识别和国防侦察，都离不开它。 由于灰度图比较简单，因此本文跳过灰度图，直接研究彩色图的颜色畸变校 正，但无论是灰度图还是彩色图，其校正原理都是一样的，程序实现上只须对程 序进行小小的调整。 如果用 Matlab 显示颜色畸变的图像 RGB 基色直方图 ，发现相对正常图像， 颜色畸变的图像的直方图的三种基色的直方图中至少有一个直方图的像素明显 集中集中在一处，或则集中在 0 处或则集中在 255 处，而另一部分有空缺，或则 集中在中间而两边空，因此通过调整该直方图的像素点的像素值在区间[0，255] 上的分布来解决图像颜色畸变问题。 如果直方图中像素集中在 0 一边则说明该基 色偏暗， 如果集中在 255 处则说明该基色偏亮。下图是一有颜色畸变的图像的基 色 B 的直方图。图1基色 B 的直方图很明显几乎所有像素点都集中在区间[a，b]上，这是偏暗的情况。那么要做 的是把代表基色 B 的矩阵的数据拉伸，使得区间[a，b]扩大为区间[a，c]。只要 做以下处理即可得到以上目的。对每一个 x ，x 在[a，b]上，x *（c-a）/（b-a） ， 而所有的 y，y 在区间[b，c]上，y=c，c=255。其算法流程图为：37试验 5用矩阵变换实现对图像形状及颜色畸变的校正开始读取偏色图像到矩阵 a取矩阵的行列数保存 在 m、n设定拉伸系数 fr1 =& i1=& J拉伸处理 bm=b(i,j,n)*N bm&255 ? Y将 bm 值付给矩阵 b：b(i,j,n)=i +1=&I j&n? Y N i&n? N 显示处理后的图像 Y j+1=&j结束MATLAB 中的算法实现如下： function dealcolor(pic,n,d) a=imread(pic); %提取指定图像到矩阵 a b=double(a); %将矩阵 a 的数据转化为 double 型38试验 5用矩阵变换实现对图像形状及颜色畸变的校正[m,n]=size(b(:,:,n)); %取图像矩阵的行列数 fr=255/d; %设定拉伸系数 for i=1:m % 二重循环 对矩阵内的每一个数据进行处理 for j=1:n bm=b(i,j,n)* % 拉伸处理 if bm&255 %将所有值大于 255 的点都设为 255 bm=255; end b(i,j,n)= end end c=uint8(b); %将矩阵 b 转化为 8 个字节的整型数据 image(c); %显示处理过的图像 注释：a ――― 要处理的图像矩阵，是一个三维矩阵 pic ――― 要处理的图像的路径 n ――― 要处理的第几个基色矩阵，1、2、3 分别代表 R、G、B d ――― 向量，它的值是要拉伸的像素值中的最大值，0~255 之间 4、图像颜色畸变校正实例 对于已经发生颜色失真的图像（下图） ，通过数学变换进行校正。图 2 发生了颜色畸变的图像 可以看出，图像明显偏黄。下面研究它的 RGB 直方图，确定拉伸系数。39试验 5用矩阵变换实现对图像形状及颜色畸变的校正图 3图 2 的 RGB 直方图由于 R、 的直方图则在 0 到 255 之间分布较为均匀， G 而在基色 B 的直方图 中，像素值的分布明显是集中在 140 以下的部分，而基本没有多少点在 140 到 255 之间，这说明基色 B 偏暗。像素值在 0~255 之间分布偏向一边是导致这幅图 像颜色畸变的原因，因此从基色 B 着手处理。 既然基色 B 直方图中像素分布不均匀，那么就应该把像素值在 0 到 140 之 间的点的像素值拉伸以至它们均匀分布在 0 到 255 上。因此做以下的处理。将代 表基色 B 的矩阵的每一个小于 140 的元素的值都乘以系数 255/140 ，显然这样 可以使区间[0，140] 拉伸为[0，255]，以达到需要的效果。因此处理函数的输入 参数为： dealcolor(‘color.jpg’,3,140); 处理后得到的图片如下：图4校正后的图像40试验 5用矩阵变换实现对图像形状及颜色畸变的校正很明显图像整体色质得到了很好的改善。 在这里参数 d 的值可以不断地变换，对比处理后的图像画质，选定一个最 好的值作为参数 d 的值来校正图像。经过测试，发现 d 的值在区间[130，160] 得到的校正图像基本没什么大的差别。 以上讨论的是单通道（单个基色直方图）向上调整（向 255 拉伸）的情况，而双 通道、 三通道需要调整的情况跟单通道的调整方法是一样的，我们只需要研究各 个直方图，拿出需要的参数值，多次调用函数 function dealcolor(pic,n,d) 即可。 而对于需要向两边调整和向下调整（向 0 拉伸）的情况则只须先将需要拉伸的区 间平移到 0 处， 再进行向上调整。这些只需要对算法做一个小小的改动就可以实 现。 function dealcolor(pic,n,d) % d 是一个包含两个元素的向量, d(1)表示拉伸区间的下界 d(2)表示拉伸区间的上界 a=imread(pic); b=double(a); [m,n]=size(b(:,:,n)); fr=255/(d(2)-d(1)); for i=1:m % 二重循环 对矩阵内的每一个数据进行处理 for j=1:n bb= b(i,j,n)-d(1); %将每个元素值减去 d(1)向 0 靠拢 if bb&0 %将所有处理前值小于 d(1)的点的处理后的值设置为 0 bb=0; end bm=bb* % 拉伸处理 if bm&255 %将所有处理后值大于 255 的点都设为 255 bm=255; end b(i,j,n)= end end c=uint8(b); %将矩阵 b 转化为 8 个字节的整型数据 image(c); 5、图像模糊校正算法介绍 图像模糊校正的应用与偏色校正的应用一样广泛而重要， 在很多领域它们具 有同样的重要地位：如医学成像、罪犯识别、国防侦察等等。虽然图像模糊校正 最先是由处理间谍卫星图像而发展起来的， 但是随着摄像设备特别是随着摄像监 控技术在社会管理中的广泛应用， 警察可以轻易从监控设备得到的犯罪罪犯图像 资料。但由于技术上的缺陷，得到的图像资料通常都是模糊的，如果不进行相关 的处理很难得到罪犯的特征，追缉工作难度也会增加。 图像的灰度变化情况可以表现为一曲线。当读入一个图像后，灰度变化就转变成 了矩阵数据的变化。 反映数据变化的数学手段可以采用微分算子。从数学的微分 含义来看， “一阶微分”是描述“数据的变化率”“二阶微分”是描述“数据变 ， 化率的变化率” 。在感应灰度变化方面，二阶微分比一阶微分更具敏感性，尤其 是对斜坡渐变的细节。因此采用二阶微分算子来处理。41试验 5用矩阵变换实现对图像形状及颜色畸变的校正最简单的各向同性微分算子是拉普拉斯微分算子。设原图为 f(x,y)，一个二 维的拉普拉斯微分算子定义为： ? 2 f ??2 f ?2 f ? ?x 2 ?y 2将它展开就得到? 2 f ? 4 f ( x, y) ? f ( x ? 1, y) ? f ( x ? 1, y) ? f ( x, y ? 1) ? f ( x, y ? 1) 。写成图像处理运? 0 ?1 0 ? 算模版的形式就是 L ? ?? 1 4 ? 1? 。 ? ? ? 0 ?1 0 ? ? ?算法的流程为：42试验 5用矩阵变换实现对图像形状及颜色畸变的校正开始 读入图像到矩 阵a 取矩阵的长宽存放到 h,w 中生成一个与矩阵 a 相同维数的矩阵 g1=&Y1=&X1=&J调用函数 pick_tem(),接收 返回值为相应输出图像点 的像素值 J+1=&JJ&k+1? X+1=&XYYX&w? Y+1=&Y Y&h?Y显示处理过的图像 g结束设处理后的图像为 g(x,y)，则 g ( x, y) ? f ( x, y) ? ? 2 f ( x, y) 。用模板表示则是? 0 ?1 0 ? L0 ? ? ? 1 5 ? 1? ，模板中心点就是要处理的像素点。如果连对角线方向都考 ? ? ? 0 ?1 0 ? ? ? ?? 1 ? 1 ? 1? 虑的话，模板的形式就表现为 L1 ? ?? 1 9 ? 1? 。 ? ? ?? 1 ? 1 ? 1? ? ?43试验 5用矩阵变换实现对图像形状及颜色畸变的校正由于模板是一个 3 阶矩阵， 所以模板处理不了图像矩阵边缘的点。将这些处 理不了的点像素值都设置为 255。 在 Matlab 中的算法实现为： function pick_tem(b,j,tem_n) % b 是一个三维矩阵，j 是代表第几层矩阵 %tem_n 代表模板号 0 代表采用第一种模板，1 代表采用第二种模板 if tem_n==0 %第一种模板 %返回二阶微分处理后的值 pic_tem=5*b(y,x,j)-(b(y,x-1,j)+b(y+1,x,j)+b(y-1,x,j)+b(y,x+1,j)); else if tem_n==1 then %第二种模板 %返回二阶微分处理后的值 pic_tem=9*b(y,x,j)-(b(y-1,x-1,j)+b(y-1,x,j)+b(y-1,x+1,j)+b(y,x-1,j)+b(y+1,x+1,j)+b(y +1,x-1,j)+b(y+1,x,j)+b(y+1,x+1,j)); end end function faintness(pic,tem_n) %pic 代表处理的图片的路径，tem_n 代表采用第几种模板 a=imread(pic); %导入图像，并把数据存放到三维矩阵 a 中 b=double(a); %将数据转化为双精度 [h,w,k]=size(b); %取矩阵的维数 g=zeros(h,w,k)+255; %生成一个与原图像矩阵有相同维数的矩阵 g for y=2:h-1 %剔除图像矩阵最外一层的点，处理可以被模板包含的点 for x=2:w-1 for j=1:k g(y,x,j)=pick_tem(b,j,tem_n); %采用二阶微分算子处理，返回处理点（x,y）后的像素值 end end end image(uint8(g));%显示处理后的图像 6、图像模糊校正实例 下面的图像是一幅因调焦不准确而导致模糊了的图像。根据上述算法，用 Matlab 对该图像进行校正。44试验 5用矩阵变换实现对图像形状及颜色畸变的校正图5需要处理的照片先采用模板 L0 处理上图，函数 faintnesss 输入参数如下： faintness（’ faintness.jpg’,0）; 处理的结果为：图6用 L0 处理后的图像再采用模板 L1 处理，函数 faintnesss 输入参数如下 faintness（’ faintness.jpg’,1）;45试验 5用矩阵变换实现对图像形状及颜色畸变的校正处理的结果为：图7用模板 L1 处理后的图像图 6 和图 7 都比原图 5 要清晰得多， 说明采用二阶微分算子来处理模糊是正 确有效的。对比图 6 和图 7，发现图 7 的线条轮廓比图 6 要清晰，但是相对图 6， 图 7 的噪音也加强了很多。 这就是采用二阶微分处理模糊图像的一个弊端――加 强噪音。 以上所研究的都是整幅图像模糊需要校正的情况，而局部模糊的图像校正 跟全局模糊校正都是一样的基理。只须将需要校正的区域提取出来，把提取出来 的图像区域作为一个图像整体看待，再用全局处理的方法来处理。那么可以根据 这个原理把函数程序 faintness 一般化到局部可选情况。 function faintness(pic,tem_n,x1,y1,x2,y2) %（x1,y1）表示局部区域左上角点 %（x2,y2）表示局部区域右下角点 a=imread(pic); b=double(a); [h,w,k]=size(b); g=b; % 替代 g=zeros(h,w,k)+255; for y=y1:y2 for x=x1:x2 for j=1:k g(y,x,j)=pick_tem(b,j,tem_n); end end end image(uint8(g)); 7、图像几何畸变校正算法介绍46试验 5用矩阵变换实现对图像形状及颜色畸变的校正几何畸变广泛存在于各种实际工程应用中，尤其是在遥感、遥测等领域。目 前, 利用数字图像技术进行畸变校正主要应用于两个领域: 一是医用内窥镜图像 的校正, 目的是测量图像上病变区域的大小; 另一应用是机器视觉领域, 是为了 提高它的视觉定位精度, 正确完成动作。几何畸变有桶形畸变、枕形畸变、几何 倾斜等等。这里选择桶形畸变的校正作为研究的对象。 对于畸变的光学系统， 畸变空间中的直线在像空间中一般不再是直线，而只 有通过对称中心的直线是例外。因此在进行桶形畸变校正时须先找出对称中心， 再进行通用的几何畸变校正过程。 桶形畸变校正一般步骤： （1）找出畸变图对称中心，将畸变图代表的地址空间关系转换为以对称中心为 原点的空间关系。 （2）空间变换：对输入图像（畸变图）上像素重新排列以恢复原空间关系。也 就是利用地址映射关系为校正图空间上的每一个点找到它们在畸变图空间上的 对应点。 （3） 灰度插值： 对空间变换后的像素赋予相应的灰度值以恢复原位置的灰度值。 几何畸变的校正要使用几何（坐标）变换，包括平行移动、旋转、扩大缩小等简 单的变换。 在这里， 先在直角坐标系下研究地址映射关系， 在程序中则使用矩阵坐标系。 设[f(x,y)]是原图，[f(u,v)]是发生畸变后的图像。 畸变校正的基本思想是，找出由(u,v)? (x,y)的坐标变换 T? (地址映射),然后 令： f ( x, y) ? f (T? (u, v)) ? 为参数向量。若畸变只是简单的纵横比的改变和倾 斜，那么仿射变换可以校正这种畸变。取： T? : ? ? ?A, B, C , D, E, F ?x ? ?Au ? Bx ? C ?y ? [ Du ? Ev ? F ]?x?表示去最接近 x 的整数若能得到参数 ? 的估计，问题就可以解决了。 但一般的畸变都不只是简单的纵横比的改变和倾斜， 通常遇到的都是空间扭 曲型几何畸变，俗语讲就是橡胶层面拉伸。它是曲线畸变，这里采用二次多项式 来模拟它,可表示为：u ? a0 ? a1 x ? a 2 y ? a3 x 2 ? a 4 xy ? a5 y 2 (1) v ? b0 ? b1 x ? b2 y ? b3 x 2 ? b4 xy ? b5 y 2 (2)同样，只要可以取得参数 (ai , bi ) 的估计，畸变函数便可知，那么原则上，可以通 过上式多项式变换来获得修正的空间扭曲映射。 (1)和(2)都是一个有 6 个参数的二元二次方程，那么只要在畸变图和校正图 上各取六对对应点 （其中从校正图上的点是估计值）就可以通过解方程组得到参 数 (ai , bi ) 的估计， 理论上， 取的对应点对数越多得到的参数 (ai , bi ) 估计就越精确。 设取 m 对对应点，用向量来表示为U t ? [u1 ,u 2 ,..., um ] V t ? [v1 , v2 ,..., v m ]47试验 5用矩阵变换实现对图像形状及颜色畸变的校正?1 x1 ? ?1 x 2 A ? ?. . ? ?. . ?1 x m ?y1 y2 . . ymx1 x2 . .2 2x1 y1 x2 y 2 . .xm2xm y m2 y1 ? 2? y2 ? . ? ? . ? 2 ym ? ?其系数 a t ? [a0 , a1 ,..., a m ]b t ? [b0 , b1 , . . bm ] .,假设所取的 m 对对应点组成的矩阵 A 可逆，也就是说 m 对对应点线性不相 关，则容易计算得到系数 (ai , bi ) ： a ? A ?Ub ? A ?V由于采用灰度插值， 所以在校正的处理过程是对校正图上每一点映射到畸变 图， 然后通过灰度插值来得到这一点的灰度值。 因此校正所采用的是地址逆映射：u ? 1, x, y, x 2 , xy, y 2 * a??v ? 1, x, y, x 2 , xy, y 2 * b??(?)由地址映射( ? )计算得到的（u,v）可能是非整数，而畸变图[f(u,v)]是数字图 像， 其像素值仅在坐标为整数处有定义，所以在非整数处的像素值要用其周围一 些整数处的像素值来计算，这叫灰度插值。灰度插值有邻近插值法、双线性插值 法等等。 邻近插值法得到的图像通常都会出现锯齿现象，而双线性插值法得到的 图像比邻近插值法得到的图像精确得多。通常情况下，双线性插值的精确度已经 可以满足一般图像处理的要求， 并不需要更高精确度的灰度插值。因此在这里将 采用双线性插值法。 双线性插值利用(u,v)周围的四个最邻近像素的灰度值，根据下面方法来计 算(u,v)处的灰度值。 设(u,v)四个邻近像素点为 ABCD 坐标分别为(i,j), (i+1,j), (i,j+1),(i+1,j+1)(i,j+1) C E (u,v) D FA (i,j)B (i+1,j)设 ? ? u ?i? ?v? j首先先计算出 E,F 两处的灰度值，f(E)和 f(F) f(E)= ? [f(C)-f(A)]+f(A) f(F)= ? [f(D)-f(B)]+f(B)48试验 5用矩阵变换实现对图像形状及颜色畸变的校正再计算(u,v)： f(u,v)= ? [f(F)-f(E)]+f(E) 此 f(u,v)值代表的就是校正后图像中(x,y)处的灰度值。 算法的流程为：49试验 5用矩阵变换实现对图像形状及颜色畸变的校正开始读入图像到矩阵 a取矩阵的长宽存放到 h,w 中生成一个与矩阵 a 相同维数的矩阵 sp 存放校正图信息求解系数估计 a0 、b01=&I1=&J构造逆向映射多项式 x=[1,j-og(1),i-og(2),(j-og(1))^2,(i-og(2))*(j-og(1)),(i-og(2))^2];用逆向映射求理想图点在失真图中的映射 u=x*a0; v=x*b0;N点（u,v）在畸变图中Y对 u,v 取整并计算参数 arf bta1=&kJ+1=&J做双线性插值YK&4? J+1=&JYJ&w+1? I+1=&IYI&h+1结束50试验 5用矩阵变换实现对图像形状及颜色畸变的校正在 Matlab 中的算法实现为： function gmodify(pic,uv,gm,og) %pic 表示要处理的图像的路径文件名 %uv 是一个二维矩阵,uv(:,1)代表上面提到的 U t ,uv(:,2)表示 V t %gm 是一个二维矩阵,gm(j,:)代表在校正图空间上与 uv(j,:)一一应的点 %og 代表对称中心，它是一个二维向量 a=imread(pic); b=double(a); n=size(gm(:,1)); for k=1:n%转换到以对称点为原点的空间关系并构造矩阵 A A(k)=[1,gm(k,1)-og(1),gm(k,2)-og(2),(gm(k,1)-og(1)^2, (gm(k,1)-og(1))*(gm(k,2)-og(2)),(gm(k,2)-og(2))^2]; end [h,w]=size(b(:,:,1)); sp=zeros(h,w,3)+255; a0=pinv(A)* uv(:,2); %计算上面提到的地址映射的系数估计 a b0=pinv(A)* uv(:,1); %计算上面中提到的地址映射的系数估计 b for i=1:h %从理想图像矩阵出发处理 for j=1:w x=[1,j-og(1),i-og(2),(j-og(1))^2,(i-og(2))*(j-og(1)),(i-og(2))^2]; u=x*a0+og(2); % 逆向映射（j，i）到畸变图像矩阵（v，u） v=x*b0+og(1); if (u&1)&&(u&w)&&(v&1)&&(v&h) %处理在图像大小范围内的像素点 uu=floor(u); %对 u 取整 vv=floor(v); %对 v 取整 arf=u- %计算上面提到的 ? bta=v- %计算上面提到的 ?for k=1:3 %进行灰度双线性插值 ft1=(1-bta)*b(vv,uu,k)+bta*b(vv+1,uu,k); ft2=(1-bta)*b(vv,uu+1,k)+bta*b(vv+1,uu+1,k); sp(i,j,k)=(1-arf)*ft1+arf*ft2; end end end end image(uint8(sp)); %显示校正图像 8、图像几何畸变校正实例 下面是一幅已经发生严重几何变形的图片，根据上面描述的算法，对该图片 进行几何畸变校正。51试验 5用矩阵变换实现对图像形状及颜色畸变的校正图8畸变图这里选择网状图是因为网状图更容易分析，校正效果也更容易观察。 观察畸变图， 发现从上数下第五条横线和从左数起第六条竖线都是直线，因 此选用这两条直线的交叉点作为对称中心进行畸变校正。然后，从畸变图中选出 4 个线性不相关的点，并估计它们在校正图上的对应点。因此得到的三个输入参 数的值分别如下： pic=’tt.jpg’; uv=[144,26;26,198;144,416;308,241;19,29;303,30;307,325]; gm=[144,29;29,198;144,412;304,241;15,26;306,26;307,327]; og=[145，241]; gmodify(pic,uv,gm,og); 得到的校正图如下：图 9校正后的图片52试验 5用矩阵变换实现对图像形状及颜色畸变的校正虽然经过校正得到的图像仍有一定程度的畸变，但是畸变程度小多了。 对几何畸变图像的校正计算, 不可能完全准确的恢复原来的图像。 这是因为影响 畸变校正精度的因素是多方面的。其中对校正算法精度影响较大的几种误差有: 对称中心估算误差、系数估计误差、插值误差和灰度校正误差等。 该校正方法关键在于对称中心的估计和地址映射的系数估计，而系数估计的 关键是在于所取对应点对的精确程度。显然，校正后的图像仍有一定程度的畸变 是由于对应点对不够精确。 校正图上的对应点的估计是很困难的，通常都需要对 程序运行的结果图进行分析，然后不断进行调整。 六、结论及应用 通过实验我们可以了解到，我们学习了数字图像的基本概念，包括图像的 数学化、采样、量化、灰度、各种数学图像的文件格式、表色系、颜色映像等。 再通过具体的实例给学生示范了对于在颜色或形状上发生畸变的图像如何通过 数学的方法实现校正的过程。 七、练习 寻找一份自己以前的已经发生颜色或形状畸变的照片，按照上面实验的原 理，做出有效的校正，要求写出实验报告。53试验 6Galton 钉板实验试验 6 Galton 钉板实验一、实验目的与要求 1、复习概率论中随机变量、概率分布、二项分布、均值和分布函数等概念。 2、理解 Galton 钉板实验中小球落入格子所服从的规律。 3、了解 Matlab 软件中进行动画演示的命令。 4、掌握 Matlab 软件中计算二项分布概率、产生二项分布随机数的命令，掌 握计算离散型随机变量数学期望的方式。 5、掌握 Matlab 软件中进行随机模拟的方法。 二、问题描述 所有现象的“因”和“果”，即“条件”和“结果”之间在客观上都存在着 一定的规律，这种规律通常可以分成两类：一是确定性的规律，另一类是非确定 性规律。 对于确定性的系统， 当已知条件是充分时， 那么实验的结果也是确定的， 即在每一次试验进行以前，可以预见试验产生的结果。但若条件不充分时，就无 法预测试验的结果，这就产生了“因果律的缺失”的随机现象。随机现象在实践 中是大量遇到的，如掷骰子。虽然无法由“因”预测“果”，但是当进行大量重 复试验时，因果之间仍会呈现一种统计规律。概率方法建立在“重复试验”的基 础上，统计规律只有在大量重复后才会