MATLAB的图像处理介绍.ppt

关于matlab中图像处理的简单介绍徐健2010年7月16日内容主要分为以下四个篇章：一.介绍应用matlab处理图像问题（以一个图像隐写的简单例子为说明线索）二.应用matlab解决三维血管重建中的图像处理问题三.图像处理练习（应用matlab解决双目定位问题中的图像处理）四.附录与说明(包括：函数命令库)一.介绍应用matlab处理图像问题（以一个图像隐藏的简单例子为说明线索）1.图像与数字图像简介2.图像的不同类型及在matlab中的显示3.常用的几个图像处理命令4.一个简单的图像隐写介绍5.图像处理中一些注意问题图像与数字图像简介n图像：利用各种系统观测客观世界获得的且可以直接或间接感知的视觉实体。信息其他(20%)视觉(60%)听觉(20%)n图像:n（1）模拟图像：光学图像、模拟电视图像等。处理速度快，但精度和灵活度差。n（2）数字图像：数码相机、数字电视等。是将连续的模拟图像经过离散化处理后得到的计算机可以识别及处理的点阵图像。数字图像像素数字图像的优点：（1）精度高：目前计算机可将模拟图像转化成高精度数字图像（2）处理方便：数字图像是一组数据，可利用计算机对其处理（3）重复性好：数字图像可比模拟图像有正常的保质时间数字图像已经应用到各个领域，无处不在。那么对数字图像的处理主要有以下方面：（1）图像变换：傅立叶变换，小波变换等。（2）图像增强与复原：突出图像信息，抗干扰。（3）图像压缩编码：简化图像利于传输等。（4）图像分割：提取图像中的有意义的特征。（5）图像分析：对图像中的信息进行各种分析。（6）图像识别：提取图像中的信息进行判别。（7）图像隐藏：对图像加入水印进行信息伪装。图像的不同类型及在matlab中的显示n1.数字图像数字图像（按纪录方式分）：（按纪录方式分）：n（1）矢量图像：利用数学的矢量方式纪录图像内容。以线条和色块为主，容易放大、缩小或旋转，且不易失真，精确度高，可以绘制3D图像。但是不易做成色彩丰富的图像。n（2）位图图像：将图像中每一个像素点转换成一个数据。如果以8位记录，可以表现出256种颜色（），所以色彩丰富。通常有：16色，256色，增强16位和真彩色24位（）.但随着颜色数和分辨率的提高，存储空间大，且较易失真。用数码相机和扫描仪获得的图像都属于位图。3）象素：是图像在计算机显示中的度量单位，可以变化，可大可小。4）分辨率：是用于度量图像在显示器中清晰程度的一个参数，分辨率越高，图像越清晰。分辨率是与象素相关的，即单位长度上的象素数就是分辨率。由此可知，分辨率越高，象素的几何尺寸就越小。5）图像文件的大小：指一幅图像在计算机中保存时所占用的磁盘空间，其大小与所用的颜色模式有关。灰度图像中的每一个灰度象素只占用一个字节（8位），RGB图像中红、绿、蓝各占用一个字节。另外，图像文件的大小也直接与其分辨率有关，原因是当分辨率增加时，一幅图像所包含的象素量急剧增加。6）句柄：就是对象的代号或标志，它能使计算机方便地找到所需要的对象并加以相应的操作。MATLAB中的句柄图形对象包括轴、文本、菜单、控制框、图像等。2.几种常见的几种常见的MATLAB 图像文件格式简介图像文件格式简介：A）BMP格式。即位图文件，整幅图可视为一个数字矩阵。它包括1、4、8、24位非压缩图像，8位RLE（行程编码）图像。文件内容包含文件头、位图信息数据块和图像数据。选择BMP格式保存一幅灰度模式图像时，可选择以Windows格式保存。而且在选中4位或8位位图时，还可选压缩（RLE）项，在用RLE方式压缩保存后图像将毫无损失。这是用得最广的图像格式之一。B）TIFF格式。处理1、4、8、24位非压缩图像，1、4、8、24位packbit压缩图像，1位CCITT压缩图像等。文件内容包括：文件头、参数指针表与参数域、参数数据表和图像数据四部分。是一种用途广泛的文件格式，其特点是可移植性好，几乎所有的扫描仪及在Windows、Macintosh平台上常用的版面设计软件都支持TIFF文件格式。但图像文件结构比较复杂，不压缩时文件比较大。C）JPEG格式。是一种联合图像专家组的图像压缩格式，是目前所用对静止灰度或彩色图像的压缩标准。它实际上定义了3种编码系统：a基于DCT有损编码基本系统，用于绝大多数压缩场合；b用于高压缩比、高精度或渐进重建应用的扩展编码系统；c用于无失真应用场合的无损系统。JPEG没有规定文件格式、图像分辨率或所用的彩色空间模型，这使它适用于MATLAB。D）PCX格式。可处理1、4、8、16、24位等图像数据。文件内容包括文件头、图像数据、扩展调色板数据。E）XWD格式。1、8位Zpixmaps,Xybitmaps,1位XYPixmaps。F）TGA格式。处理1、4、8、16、24位非压缩图像和行程编码图像。文件包由5个固定长度字段和3个可变长度字段组成。G）HDF格式。有8位，24位光栅图像数据集。3.MATLAB图像文件类型：图像文件类型：根据数据矩阵和图像象素颜色匹配关系，MATLAB中图像可分为：索引图像、灰度图像、二值图像和RGB图像。1）索引图像：它的数据信息包括一个数据矩阵和一个双精度色图矩阵，它的数据矩阵中的值直接指定该点的颜色为色图矩阵中的某一种。色图矩阵中，每一行表示一种颜色，每行有三个数据，分别表示该种颜色中红、绿、蓝的比例情况，所有元素值都在0，1内。数据矩阵数据矩阵（图片）（图片）像素点对应双精度色图矩阵双精度色图矩阵（色彩）（色彩）2）灰度图像：数据矩阵中的元素值一般都在0,1或0,255之间，灰度图像根据这些数据利用线性插值来和色图中的颜色种类匹配。灰度图像读入matlab中是一个二维的平面矩阵，其中行与列的乘积代表其图片中像素点的个数。注意：灰度图像一般看起来是一副黑白图像，但是色彩明暗度较二值图像更为丰富。因为每一个像素点的取值在0,1或0,255之间。.像素点取值：0,1或0,255mnm*n个像素，如1024*7683）二值图像：数据矩阵中的元素值只是0或1。读入matlab也是一个二维矩阵。注意：二值图像读入matlab中也是一个二维的平面矩阵，但像素点取值只限于0，1。4）RGB图像：图像中每个象素的颜色用三个数据来存储，分别指定红、绿、蓝三原色在象素颜色中的比例关系，组成一个三维数组，读入matlab后是一个三维的矩阵。注意：美术教科书中称红、黄、蓝为三原色，讲的是绘画颜料的使用。一般电视光色等光色是红、绿、蓝。RGB图像就是采用红、绿、蓝作为三原色的，其中R为红色，G为绿色，B为蓝色。上图是一个2048*1536大小的图像，其中这个三维矩阵的第一维就是上图中第一层代表红色数值，第二维为第二层代表绿色数值，第三维为第三层代表蓝色数值。也可以这样理解：将索引图像中的数据矩阵中每一个像素点直接加载上色图矩阵中对应的颜色值。图像在图像（x1,y1）点的RGB值是（r,g,b）n4.matlab中图形图像的读入中图形图像的读入n在matlab中利用函数imread将图像读成一个矩阵的形式。其主要格式如下：nA=imread(filename,fmt)nX,map=imread(filename,fmt)n=imread(filename)n=imread(URL,)n=imread(,idx)(CUR,ICO,andTIFFonly)n=imread(,frames,idx)(GIFonly)n=imread(,ref)(HDFonly)n=imread(,BackgroundColor,BG)(PNGonly)nA,map,alpha=imread()(ICO,CUR,andPNGonly)n其中最基本的格式是：X,map=imread(filename,fmt)存储图像数据的矩阵名图像调色板图像文件名文件格式常用n一个例子：图像Imread读入Matlab中的矩阵A（一个三维矩阵）Matlab窗口Imread命令读入的A是一个三维的数据为uint8位的矩阵命令窗口n对A（1536*2048*3，uint8）的解释如下图：通过左图的表示，这样这个三维矩阵A就可以表示成一个彩色矩阵，也就是一张数字图像可以在matlab中读成一个矩阵A值得注意的是数据类型，上面记录的数据是uint8型，关于数据类型，有如下内容：在图像（x1,y1）点的RGB值是（r,g,b）且数据为uint8位5.MATLAB中图像的存储运算中图像的存储运算将一幅图片读入到MATLAB中，其数值一般都采用double型（64位）存储和运算。但为了节省存储空间，MATLAB提供了特殊的数据类型uint8（8位无符号整数），以此方式存储的图像称为8位型像。函数image能够直接显示8位图像，但8位型数据和double型数据在image中意义不一样。对于灰度图像，uint8表示范围0,255，double型表示范围0,1。可见，double型和uint8型灰度图像不一样，二者转换格式为：I8=uint8(round(I64*255);I64=double(I8)/255;反之，imread根据文件中的图像种类做不同的处理。当文件中的图像为灰度图像时，imread把图像存入一个8位矩阵中，把色图矩阵转换为双精度矩阵，矩阵中每个元素值在0，1内；当为RGB图像时，imread把数据存入到一个8位RGB矩阵中。还有另一种格式是unit16型，在这里就不做更多的介绍了，但是要注意的是，有时针对一幅图像进行相应处理时，需要对数据类型进行转换，这里就需要注意很多匹配的问题。n前面介绍了图像的读入，下面看看输出：n6.matlab中图形图像的输出中图形图像的输出n（1）以图像文件的形式输出，应用函数imwriten格式是：nimwrite(A,filename,fmt)nimwrite(X,map,filename,fmt)nimwrite(,filename)nimwrite(,Param1,Vall,Param2,Val2,)n备注：1）命令中各参数含义可以参照前面imread命令，也可在matlab命令窗口输入：helpimwrite，来获得提示信息。n2）imwrite获得的文件存放在matlab的work文件夹中。一个图片文件矩阵Aimwriten（2）以图像的形式输出，函数为image(imshow)nimage的格式是：nimage(C)nimage(x,y,C)nimage(,PropertyName,PropertyValue,)nimage(PropertyName,PropertyValue,)nhandle=image()n其中，x,y分别表示图像显示位置的左上角坐标，C表示所需显示的图像。函数imagesc与image函数类似，但是可以自动标度输入数据。n关于image中各参数的意义可以使用helpimage查询Matlab输出的图像矩阵AImage/imshown一个image的例子：得到图像输出窗口点击DataCursor按钮并将光标移动到图像的一个位置可以得到此处的数据信息，如右图所示：nimshow的格式是：nimshow(I,n)nimshow(I,lowhigh)nimshow(BW)nimshow(X,map)nimshow(RGB)nimshow(,display_option)nimshow(x,y,A,)nimshowfilenamenh=imshow()n关于imshow的参数说明可以参照image，以及在命令窗口可以得到相应的例子。imshow和image:图像的显示是最为重要的，用imshow和image都可以显示图像，但是有一定的区别。用的不对，可能出错或得到一张空白图或者是彩色图显示成颗粒状、反相黑白图等等。image是用来显示附标图像，即显示的图像上有x,y坐标轴的显示，可以看到图像的像素大小，而imshow只是显示图像。它们都可以用subplot来定位图像显示的位置。image显示imshow显示显示真彩色图像像素如果是uint8类型，要求数据范围为0-255，如是double型，则其数据范围为0-1。uint16数据类型与uint8类似，取值范围为0-65536。常用的几个图像处理命令介绍nimfinfo函数可以观察图像信息文件此处没有分号得到图像信息Colorbar:Displaycolorbar(MATLABfunction)显示彩条Imagesc:Scaledataanddisplayasimage(MATLABfunction)缩放数据并显示图像Immovie:Makemoviefrommultiframeindexedimage由多帧图像制作图像Imtool:DisplayimageintheImageViewer在imagetool中显示图像Montage:Displaymultipleimageframesasrectangularmontage将多个图像帧显示为矩形蒙太奇Subimage:Displaymultipleimagesinsinglefigure在单个图形中显示多个图像（图像处理函数详见附表）特别注意一些转换函数：（1）数据转换函数：）数据转换函数：A=unit8(B)、B=double(A)或或B=A/255前面已经介绍。前面已经介绍。Im2double，im2unit8,im2unit16分别是将图像矩阵数值转换为分别是将图像矩阵数值转换为double型，型，unit8型，型，unit16型。型。（2）图像类型转换函数：）图像类型转换函数：索引图像灰度图像二值图像RGB图像Gray2ind Ind2gray Int2rgb Rgb2gray Rgb2ind 三层变一层三层变一层一层变三层一层变三层（方法研究）（方法研究）Im2bw 涉及到阈值选择，即多少作为黑白图像的分界线，一个阈值选择函数：graythresh。一个转换函数例如：W=im2bw(W,graythresh(W)一个简单的图像叠加隐写n1.图像隐藏是指媒体信息的相互隐藏，常见的有数字水印和图像的信息伪装等。n2.把一幅图像中隐藏一个信息，并可以通过一种方式将图像信息解码出来。A（原始图像）B（需要隐藏的信息）C（加入信息的图像）要求：A与C尽量相近，没有差别n3.介绍的目的：体会matlab中关于imread和image的用法。n4.采用一种原始的图像叠加的简单思想处理：（1）隐写算法：n1）将原始图像A和隐写信息B读入matlabn2）将隐写信息B扩充到与A有相同大小的点阵n3）将矩阵A和B中的数据转换成uint8型n4）对矩阵B的数据做如下处理：B(i,j,k)/100;ni,j分别是点阵的大小，k=1,2,3；n5）将矩阵A和B做叠加处理得到带有隐写信息的图像Cn（2）解码算法是隐写算法的逆过程xx叠加算法的具体图形表达5.算法可成立的支持理由：对于隐写信息B由于将象素点的RGB值减小100倍，再叠加到原始图像A中，这样基本看不出隐写信息；在解码中将处理后带有隐写信息的图像C减去A并将象素点扩大100倍可以得到隐写信息B。6.算法的实现结果：（1）隐写结果：+n（2）解码结果n7.算法中图像叠加和图像相减语句可以采用编程实现，也可采用函数imadd和imsubtract。n8.这种简单算法的问题：n1）由于采用uint8型（0-255）所以当数据叠加超过255时只以255计算，这样可能在恢复隐写信息时造成信息失真，特别对于白色图像（255）相加。n2）由于在前面隐写时像素值缩小了100倍在uint8型下省去了小数点后的数，故而得到恢复图像的像素值都为0，100，200。所以有很大的局限性。n9.可以想到的解决方向：n1）采用double型来保持小数点后面的值。n2）传输处理后图像的同时传输一个记录数据的矩阵。n3）隐写算法（已经成为一个独立的研究方向）一个较为正规的数字水印程序rgb2grayim2bw二维余弦变换（二维余弦变换（DCT2）在低频部分加载水印信息在低频部分加载水印信息公开图像公开图像灰度图像灰度图像水印图像水印图像二值图像二值图像含有水印图像含有水印图像利用二维余弦反变换利用二维余弦反变换（IDCT2）得到水印信息图像得到水印信息图像具体程序具体程序图像处理中一些注意问题当前路径存放变量默认路径改变路径各种图片和编写的M文件和程序放在默认路径n图像保存中的失真问题：图像保存中的失真问题：A=imread(山水照片.jpg);imshow(A);得到：大小为大小为2048*1536用此处保存得到图片再读用此处保存得到图片再读入和入和imshow得到得到保存白边，大保存白边，大小为小为335*470解决方案：程序中直接采用解决方案：程序中直接采用imwrite写入磁盘，这样大写入磁盘，这样大小就不会改变。一般小就不会改变。一般bmp格格式较不易失真。式较不易失真。最常碰到的问题有：图像读入imread，应是上文提过的MATLAB支持的7种格式之一，显示图像用imshow(h)语句，h为图像句柄；输出图像若需要永久保存，可用imwrite(h,map,filename.bmp,bmp),写入存储器。注意在该语句前要设置调色板，即map=(gray(256)。要注意图像格式的转化。不同的图像格式对应不同的处理方式，如果处理与格式不符，将引起错误。尽管MATLAB允许未定义使用数组，但在实际应用中这样经常出错，特别是遇到在double和uint8型之间的转换时。所以最好还是养成用前定义的习惯，避免出现不必要的错误。在做完一定量的运算后，一般要用Clear清除内存变量，以防影响后面的程序运行。一般默认路径在matlab安装文件下的work文件夹中，程序和资源都放在此，如果需要引用新的文件，应事先指明路径。应注意语句结尾“;”的使用，特别是图像处理中。应做好程序的注释工作。二.应用matlab解决三维血管重建中的图像处理问题n1.问题的来源及说明（涉及的是图像处理方面）n2.多幅图像的批量读入解决方案n3.重建方法一：最原始的，采用图像数据扫描算法n4.重建方法二：简单介绍寻找各图像中包络球算法，模拟半径和位置进行拟合的算法。n5.说明（这个问题的处理思想各异，本次课程介绍了一种方法，大家可以自己开拓思路想出更好的方法。）问题的来源及说明n题目来源：题目来源：2001年全国大学生数学建模竞赛年全国大学生数学建模竞赛A题。题。n（题目及附件）（题目及附件）n题目说明：本次课程主要侧重于应用题目说明：本次课程主要侧重于应用matlab对图像对图像文件的建模和处理，所以在本题中涉及到的三维数文件的建模和处理，所以在本题中涉及到的三维数学图形里如球的包络，以及拟合中轴线等可能介绍学图形里如球的包络，以及拟合中轴线等可能介绍的较为简略。的较为简略。A题题 血管的三维重建血管的三维重建根据拍照得到的平行切片数字图象，运用计算机重建三维形态。血管可视为一类特殊的管道，该管道的表面是由球心沿着某一曲线（称为中轴线）的球滚动包络而成。例如圆柱就是这样一种管道，其中轴线为直线，由半径固定的球滚动包络形成。现有某管道100张平行切片图象，格式均为BMP，512*512象素。假设：管道中轴线与每张切片有且只有一个交点；球半径固定；切片间距以及图象象素的尺寸均为1。计算管道的中轴线与半径，给出具体的算法，并绘制中轴线在XY、YZ、ZX平面的投影图。多幅图像的批量读入解决方案一幅图片的读取使用matlab，在本题中有100幅图，如何将这100幅图自动读入各对应的矩阵中，以便进一步处理，方案如下：（1）将100幅图放入work文件夹等相应的操作文件夹中。（2）使用for循环结合字符串和数值之间的转换函数，循环各图片文件名。代码如下：forb=0:99m(:,:,b+1)=imread(int2str(b),.bmp);endb为数值，可以循环将b转换成字符串，用于标示图片名图片文件格式由于图片文件是一个二值图像，每一幅图片文件读入到matlab中是一个二维矩阵，且数值为0或者1，代表黑和白。所以上面的代码，最终把每一张图片的数据都记录到一个三维矩阵m中，其中m每一层记录一张图片的数据，一共100层。矩阵m第100层m(:,:,100)矩阵m第50层m(:,:,50)矩阵m第1层m(:,:,1)矩阵m第1幅图（0.bmp）第50幅图（49.bmp）第100幅图（99.bmp）注意：文件读入时间较长，为了以后操作方便，可以将生成的矩阵注意：文件读入时间较长，为了以后操作方便，可以将生成的矩阵m保存为保存为.mat文件，以便以后再次调用。文件，以便以后再次调用。重建方法一：最原始的，采用图像数据扫描算法n有了图片的数据矩阵m后，可以重建图像，由于m的每一层都记录了一副图片，最容易想到的方法是扫描算法，即利用matlab扫描出m的每一层，将离散的100幅图立体再现出来，这样就可以重现血管了。nfork=0:99%k循环m的层数，即每一幅图片nfori=1:512%i,j循环每一层中的数据值nforj=1:512nif(m(i,j,k+1)=0)%如果数据是0，即图片上的黑色nplot3(i-257,j-257,k+1,b-.);nholdon%在相同位置画蓝色，并保持画点nend,end,end,endnholdoff减减257的原因的原因在于在于:依题意，依题意，每幅图原点在每幅图原点在中心，而中心，而i,j是是从从1开始，故开始，故而要减去而要减去257得到三维血管立体效果图：注意：三维扫描速度较慢，可以采用采样描点的方法，例如将步长增为5来进行描点。fori=1:5:512,forj=1:5:512在程序中可以使用rotate3d，或者按figure上的按钮来手工旋转三维图像，得到xoy，xoz和yoz上的投影图。或者编写程序得到各面上的投影图如下：旋转按钮旋转按钮xoy平面平面投影图投影图xoz平面平面投影图投影图yoz平面平面投影图投影图n=100;m1=size(m,1);m2=size(m,2);z=ones(1,1,n);for i=1:m1 for j=1:m2 if isequal(m(i,j,:),z)%注意：注意：if m(i,j,:)=z不行不行 plot(i-257,j-257,b.),hold on;end,end,endtitle(xoy平面投影图平面投影图);hold offn=100;m1=size(m,1);m2=size(m,2);z=ones(m1,1);for i=1:n for j=1:m2 if isequal(m(:,j,i),z)plot(j-257,i,b.),hold on;end,end,endtitle(xoz平面投影图平面投影图);hold offn=100;m1=size(m,1);m2=size(m,2);z=ones(1,m2);for i=1:n for j=1:m1 if isequal(m(j,:,i),z)plot(j-257,i,b.),hold on;end,end,endtitle(yoz平面投影图平面投影图);hold off关于程序的简单说明：（1）为了提高程序的扩展性，对于512和100，应用函数size进行替换，如：fori=1:size(m,1),forj=1:size(m,2);（2）为了提高运算，利用一个变量首先进行记录，以替代每次计算，如：m1=size(m,1),m2=size(m,2),fori=1:m1,forj=2:m2（3）判断两个数组是否相同时，判断条件的不可使用=或者=，因为0,1,1,0和1,1,1,1在判断时结果都为假。而是使用函数isequal判断两数组相同，用isequal判断两数组不同。关于扫描算法的总结：优点：思路简单，操作方便，代码容易编写。缺点：多维扫描速度非常缓慢；得到的图像数据点叠加较密，不易观察；由于是扫描的结果，得不到半径与中轴线。重建方法二：简单介绍寻找最大内切圆和拟合中轴线的算法。目的：目的：要得到管道的半径和绘制出中轴线及其投影。分析：分析：管道可以看成是一球体沿空间一曲线（中轴线）移动所形成。球体被一截平面连续截取所得的圆在该截平面上形成的包络即为截面。由截面形成过程可知，截得的包含球心的圆就是截面的最大内切圆，其圆心坐标为截面与中轴线的交点坐标。通过求最大内切圆的平均半径得出球体半径，即管道半径。按上述同样思路求得所有切片与中轴线的交点，将所有的交点拟合后便可得到中轴线的近似轨迹。说明说明：（1）各类定理需要证明，当年优秀的论文都对相应定理进行了理论的证明，增加了论文的质量。（2）由于程序较长，所以ppt上不再列出。假设：（1）若管道中轴线与每张切片有且只有一个交点；（2）球半径固定；（3）切片间距以及图象象素尺寸为1；（4）血管无断裂无突变，即管道表面光滑且连续；（5）对切片拍照的过程中不存在误差，数据误差仅与切片数字图象分辨率有关；图像轮廓线形成的图像轮廓线形成的图形，反映血管造图形，反映血管造型更为清晰型更为清晰下面的工作就是寻找每一幅图像截面的最大内切圆。下面的工作就是寻找每一幅图像截面的最大内切圆。步骤具体程序是扫描轮廓线算法思想如下（以第算法思想如下（以第66幅图为例）：幅图为例）：利用利用edge找出图像边缘找出图像边缘找出四个顶点，以及最小图像区域找出四个顶点，以及最小图像区域在此区域内扫描截面内部的点，在此区域内扫描截面内部的点，找到这个点到边界的最小距离，找到这个点到边界的最小距离，当取完所有内部点后所有这些当取完所有内部点后所有这些最小距离中的最大者，就为最最小距离中的最大者，就为最大内切圆半径，对应的点就是大内切圆半径，对应的点就是中轴线与截面的交点。中轴线与截面的交点。具体程序具体程序一个存在的问题：一个存在的问题：有相同的最大内切半径，都可以作为待选对象。一个截面图中可能出现多个最大内切圆，其半径相同，但对应的中心点不同。在130图中这样的情况不多，因为截面近似于一个圆形，但是在70100图中出现待选的点情况很多，使用哪一个内切圆都是符合道理的，所以这里选择的方法也比较多。程序中选择的方法是取中间位置的圆心。选择选择c1=sortrows(jilu,3);p,q=size(c1);c2=find(c1(:,3)=c1(p,q);c3=c1(c2,:);c3(round(size(c3,1)/2),:)所得结果：所得结果：内切圆心坐标（x,y）内切圆半径r第1幅图（0.bmp）（-161，0）29.0689第2幅图（1.bmp）（-161，0）29.0689第99幅图（98.bmp）（189，-1）30第100幅图（99.bmp）（189，-6）30得到中轴线坐标后，直接作中轴线图形以及在各平面上的投影效果：得到中轴线坐标后，直接作中轴线图形以及在各平面上的投影效果：三维中轴线三维中轴线中轴线在中轴线在xoy面投影面投影中轴线在中轴线在yoz面投影面投影中轴线在中轴线在xoz面投影面投影具体程序具体程序注意：通过图像可注意：通过图像可以看见，得到的曲以看见，得到的曲线不光滑，因为这线不光滑，因为这是通过是通过plot3和和plot函数得到，只函数得到，只是散点的连接，所是散点的连接，所以下面为了得到更以下面为了得到更准确的模拟图，还准确的模拟图，还要进一步对数据进要进一步对数据进行拟合。行拟合。对得到的中轴线坐标进行拟合的方法，在当年的论文中很多，这里对得到的中轴线坐标进行拟合的方法，在当年的论文中很多，这里只介绍使用多项式进行拟合，只介绍使用多项式进行拟合，matlab中命令是中命令是polyfit。m从零递增时，但当但当m达到某一个达到某一个值值M之后，之后，对对mM，有，有即即m继续继续增加增加时时与与的的值值相差不大。此相差不大。此时时就是所要确定的最小二乘多就是所要确定的最小二乘多项项式次数。式次数。若若m过小则可能不足以表示函数的固有规律，而过小则可能不足以表示函数的固有规律，而m过大时又失去了过大时又失去了光滑性，且产生抖动现象，实践中往往采用下面的方法确定最小二光滑性，且产生抖动现象，实践中往往采用下面的方法确定最小二乘多项式的次数乘多项式的次数m，令：，令：显著减少。显著减少。根据上述理论定出多项式阶数：根据上述理论定出多项式阶数：f1为为8阶多项式，阶多项式，f2为为6阶多项式。阶多项式。YZ平面上：y=f1(z)XZ平面上：x=f2(z)Xy平面上：由上面两式联立得到。具体程序具体程序绘图结果：绘图结果：拟合中轴线在拟合中轴线在yoz面投影面投影拟合中轴线在拟合中轴线在xoz面投影面投影拟合中轴线拟合中轴线具体程序具体程序拟合中轴线在拟合中轴线在xoy面投影面投影将拟合前和拟合后的图放在一起比较：将拟合前和拟合后的图放在一起比较：带有拟合中轴线的轮廓切面图空间效果图带有拟合中轴线的轮廓切面图空间效果图最终管道半径的计算方法：各个截面最大内切圆半径的平均由于前面模拟出来的血管图像都是原始切片图，现在考虑已经得到了血管拟合中轴线以及血管半径，那么可以得到一个较为形象的模拟血管，下面就模拟血管。下面就使用下面就使用matlab结合圆形算法和曲面管道画法得到图形：结合圆形算法和曲面管道画法得到图形：下面是将中轴线也标注在血管中得到图形：下面是将中轴线也标注在血管中得到图形：下面是将中轴线也标注在血管中并将血管半透明化得到图形：下面是将中轴线也标注在血管中并将血管半透明化得到图形：n灵敏度分析的思路与方法（1）由于机器中的像素点是由单个像素构成所以是一种离散的状态，所以距离以及边缘的判别都存在误差，可以增加灵敏度分析，思路是：在取边缘和中轴线坐标时建立扰动机制，分别把像素点活动一个位置，有8种选择，再对结果进行分析。（2）图像中的截面，圆面和距离都是由像素点计算和表示的，存在误差，也可以采用上述的灵敏度分析的方法进行判别。12345678（3）待选点理论上都是中轴线的坐标点，都可以选择，所以选择那个点作为坐标点也会引起误差，灵敏度分析的方法可以采用逐个实验，再分析结果。其他算法的总结题目所使用的工具多种多样，有matlab，mathmatic，fortain语言和c语言，也有用photoshop叠加100幅图的。题目所使用的方法也很多：（1）有用最小二乘法拟合每张图切片的上下曲线，并计算图片与实物的实际比例的。（2）求最大内切球方法：血管中轴线与截面的交点必落在截面边界上的一对点的连线的中点上，过这对点的截面边界的两条切线是相互平行的。也就是最大内切球的原理。（3）除去“弯月”的两端外，图形的宽度保持不变，这一宽度就是直径，选择均匀的中段在左侧定下A，在另一边找B可以得到直径。（4）先求粗略的半径最大值rmax和最小值rmin，再求出切片的边缘点的位置，根据rmin,求出切片内所含的标准圆的圆心O可能存在的位置区域，用二分法不断缩小标准圆半径r0可能的取值范围，从而求出r0。最后用求出切片内所含的标准圆的圆心的位置。（5）将图的边界点分为上弧、下弧。对于上弧固定一点，求出该点到下弧点的距离，当取遍上弧中每一个点，得到的最短距离就是管道直径；处理多个待选点的方法：取一个待选点到其余待选点距离之和最小的点作为截面的中心坐标点。（6）最大覆盖法：找到的最佳圆心位置和半径长度是由此所形成的圆面能最大面积的覆盖截面的图形，具体实现是对区域内部的点进行遍历找到能覆盖面积最大的内切圆。（8）此外还有截面投影法；搜索点法；搜索点法里的步长减半法；边缘相交相切法，最大最小半径法；枚举法；外推法；滚球法；变换法；细化法等大家也可以拓展自己的思路想出更好的方法。（7）投影法：将各截面的图像叠加在xoy平面上，形成血管再xy平面上的投影，其中心线就是血管中轴线在xy平面上的投影。类似的再xoz和yoz投影。并在投影中计算半径值。在拟合的方法选取中，方法也比较多，多是采用多项式拟合，并计算了多项式的阶，阶不同结果也不同；也有的采用Bezier曲线进行拟合。求解管道半径有多种方法，建模竞赛的优秀论文中以“平均法”和“抽样法”居多。平均法是求出每张图像内的最大内切圆半径，再取管道半径为它们的算术平均值。但此方法是一种搜索算法，圆心遍历区域内点，半径由小到大搜索，其计算量相当庞大。抽样法就是利用滚动球半径是常数，取前几片横断面图像内的最大内切圆半径的平均值为管道半径的值。此方法计算较少，但是抽样样本数目选取的合理性较难确定，存在计算误差。一个启发算法：先将100张图片叠合形成一张图片，由于切片垂直Z轴的，此图片是血管在XOY面上的投影图像。因此，它也是滚动球在XOY面上的投影滚动圆，沿中轴线的投影线滚动形成的二维包络图，且滚动球的半径与滚动圆的半径相同，因此只需求出滚动圆半径。具体图像叠合的方法可以用matlab，也可以用Photoshop控件进行叠加。然后在此投影面上使用上下弧搜索的方法得到半径。这样做的好处是：很多方法是对一张图片（而非叠合图片）进行扫描，计算出半径值后在对所有半径再取均值。相对而言，此方法在计算量和精确度上都较为优化。三.图像处理练习（应用matlab解决双目定位问题中的图像处理）n1.问题的来源和说明n2.题目的进一步描述和操作提示n3.问题的讨论n4.练习要求问题的来源和说明本题出自2008年全国大学生数学建模竞赛A题（题目）题目说明：本题不仅涉及到从给定的图像中提取信息，还要涉及到坐标变换、位置成像、建立坐标系等方面的内容。题目的进一步描述和操作提示题目所要解决的问题：（1）建立数学模型和算法以确定靶标上圆的圆心在该相机像平面的像坐标,这里坐标系原点取在该相机的光学中心，x-y平面平行于像平面；（2）对由图2、图3分别给出的靶标及其像，计算靶标上圆的圆心在像平面上的像坐标,该相机的像距（即光学中心到像平面的距离）是1577个像素单位(1毫米约为3.78个像素单位)，相机分辨率为1024768；（3）设计一种方法检验你们的模型，并对方法的精度和稳定性进行讨论；（4）建立用此靶标给出两部固定相机相对位置的数学模型和方法。任务：确定物体某任务：确定物体某些特征点的位置些特征点的位置目标点在目标点在相机一上的坐标相机一上的坐标两部相机两部相机精确的相对位置精确的相对位置目标点在目标点在相机二上的坐标相机二上的坐标相机一相机一相机二相机二方法：方法：两部相两部相机双目定位法机双目定位法几何方法几何方法结果：结果：特征点在一部特征点在一部相机坐标系中的坐标相机坐标系中的坐标特征点的位置特征点的位置关键：系统标定关键：系统标定物平面上若干个圆物平面上若干个圆（称为靶标，圆心为目标点）（称为靶标，圆心为目标点）目标点在相机一目标点在相机一像平面上的像点像平面上的像点目标点在相机二目标点在相机二像平面上的像点像平面上的像点几何方法几何方法失真失真精确地找到精确地找到靶标上这些圆的像的圆心靶标上这些圆的像的圆心一个1008*894的bmp位图，为原始图像用相机拍摄的相图，由于相机的分辨率是1024*768，所以此图的大小也是1024*768像素。两个图像文件都有自己的坐标系，所以如果要将两个文件放在一起处理，就需要建立一个统一的坐标系，也就是两个坐标系的换算。下面提出一个简单的思路：下面提出一个简单的思路：(大家自己拓展思路去思考其他的方法大家自己拓展思路去思考其他的方法)已建立统一的坐标系下进行：将原图和相图分成五块区域，分别进行考虑，下面对其中的一个区域进行分析：这一块matlab程序的实现可以采用边缘函数edge，然后扫描建立相应的区域矩形，参考三维重建算法。a=imread(yuntu.bmp);b=edge(rgb2gray(a);imshow(b)注意：由于图片是rgb格式下的bmp文件，所以使用边缘函数edge时要先将rgb文件转换成灰度图像（rgb2gray），但我们可以看见这样图像边缘不完整，这是由于图像转换所致，可能会影响到下面计算的精确度，如果要得到完整的边缘可以不用edge函数，而直接用程序扫描边缘并记录到矩阵中。一种扫描纪录图像边缘算法：一种扫描纪录图像边缘算法：For循环，图像的这一行上颜色都没有发