Matlab公路裂缝检测

一、简介

1 案例背景

随着国家对公路建设的大力投入，我国的公路通车总里程己经位居世界前列，这样进一步促进了我国经济建设的发展。随着公路的大量投运，公路日常养护和管理已经成为制约公路运营水平提高的瓶颈，特别是路面状态采集、检测维护等工作更是对传统的公路运维模式提出了挑战。路面裂缝是公路日常养护管理中最常见的路面损坏，也是影响公路状态评估和进行必要的公路维修的重要因素!。一般而言，如果路面裂缝能够在被恶化成坑槽之前得到及时修补，则可以大大节约公路的养护成本。传统的公路裂缝检测主要是人工检测，需要配置一定规模的人力、设备等资源来进行定期巡检。但是，面对日益增长的公路建设需求，人工检测具有运营效率低、主观性影响大、危险性较高等不足，已无法满足公路破损快速检测的要求。

随着计算机硬件设备和数字图像处理技术的发展，基于视觉的目标定位及检测技术也取得了不断的进步，由于其具有定位准确、检测快速、自动化操作、易于安装部署等特点，已经广泛应用于工业自动化检测过程，特别是在目标表面质量检测、目标物测量等领域的应用。因此，基于数字图像的路面裂缝检测技术可以提供一个安全、高效、成本低廉的道路状态监控服务，已有多种图像处理方法应用于路面裂缝检测并在一定程度上取得了实际应用。

2 理论基础

路面裂缝检测从视觉上来看是典型的线状目标检测，因此路面裂缝图像的增强与定位属于线状目标检测的研究领域，路面裂缝与一般线状目标相比，具有其自身的特点：目标宽度相对较小、图像对比度较低、具有自然间断、具有分叉和杂点等，并且路面裂缝只是在视觉总体上呈现出线状特征.传统的裂缝自动检测算法，如基于阈值分割、边缘检测、小波变换等算法，往往都假设路面裂缝在整幅图像中具有较高的对比度和较好的连续性，但这种假设在实际的工程项目中往往不成立。由于受拍摄天气、路面损耗、裂缝退化等因素的影响，有一定比例的裂缝相对于路面背景具有极低的对比度，这也会引起传统裂缝检测算法的失效，因此需要在裂缝图像处理前加入一定的预处理步骤。图像预处理一般是应用于图像识别、图像表示等领域的一种前期处理。在图像的采集和传输过程中，往往会因为某些原因导致图像质量降低。

例如，从视觉主观上观察图像中的物体，可能会发觉其轮廓位置过于鲜艳而显得突兀：从被检测目标物的大小和形状来看，图像特征比较模糊、难以定位；从图像对比度的角度来看，可能会受到某些噪声的影响；从图像整体来看，可能会发生某种失真、变形等，因此，待处理图像在视觉直观性和处理可行性等方面可能存在诸多干扰，我们不妨将其统称为图像质量问题。图像预处理正是用于图像质量的改善处理，通过一定的计算步骤进行适当的变换进面突出图像中某些感兴趣的信息，消除或降低干扰信息，如图像对比度增强、图像去噪或边缘提取等处理印，一般情况下，由于裂缝图像的采集需要涉及室外作业，所得图片难免会存在一定的噪声干扰、畸变等各种问题，直接进行裂缝目标的检测和提取往往会遇到困难。因此，本案例首先将裂缝图像进行预处理，改善图像质量，进而提高实验的优化效果。图像预处理的基本方法有图像灰度变换、频域变换、直方图变换、图像去噪、图像锐化、图像色彩变换等。本案例将选择其中的部分方法来进行裂缝图像的预处理操作。

2.1 图像灰度化

自然界中绝大部分的可见光谱均能通过红®、绿(G)、蓝(B)三色光按不同比例和强度进行混合而得到， 我们将其称为RGB色彩模式。该模式以RGB模型为基础， 对图像的每个像素值的RGB分量均分配一个Uint 8类型(0~255) 的强度值。例如， 纯红色的R值为255， G值为0， B值为0； 品红色的R值为255， G值为0， B值为255。RGB图像的红、绿、蓝分量各占8位，因此是24位图像，并且不同亮度的基色混合后，会产生出256x256x 256-16777216种颜色。RGB模型图形化表示如图所示。

假设F(i， j) 为RGB模型中的某像素， 若其3种基色的亮度值相等， 则会产生灰度颜色，将该R=G=B的值称为灰度值(或者称为强度值、亮度值)。因此，灰度图像就是包含多个量化灰度级的图像。假设该灰度级用Uint 8类型数值表示， 则图像的灰度级就是256(即2°=256)。本案例所选择的灰度图像灰度级均为256，其像素灰度值为0~255的某个值，当亮度值都是255时产生纯白色，当亮度值都是0时产生纯黑色，并且亮度从0到255呈现逐渐增加的趋势。RGB图像包含了由红、绿、蓝三种分量组成的大量的色彩信息， 灰度图像只有亮度信息而没有色彩信息。针对路面裂缝图像的检测要求，一般需要去除不必要的色彩信息， 将所采集到的RGB图像转换为灰度图像。RGB图像的灰度化方法有以下几种。

（1）分量值

选取像素F(i.j)的R、G、B分量中的某个值作为该像素的灰度值，即

式中，Fg(i，j)为转换后的灰度图像在(i.j)处的灰度值。

（2）最大值

选取像素F(i，j)的R、G、B分量中的最大值作为该像素的灰度值，即

（3）平均值

选取像素F(i，j)的R、G、B分量的亮度均值作为该像素的灰度值，即

（4）加权平均值

选取像素F(i，j)的R、G、B分量的亮度加权均值作为该像素的灰度值，权值的选取一般是根据分量的重要性等指标，将3个分量以加权平均的方式进行计算得到灰度值。人眼在视觉主观上一般对绿色分量敏感度较高， 对蓝色分量敏感度较低， 因此对RGB三分量进行加权平均能得到较合理的灰度图像，常用的计算公式如下：

F.(i.j)=0.299R(i，j)+0.587G(i，j)+0.114*B(i，j)

采用加权平均计算灰度图像的方式对裂缝图像进行灰度化，所得结果如图所示。

2.2 图像滤波

裂缝图像在采集或传输的过程中往往会受到成像设备与传输介质等因素的干扰而产生噪声，因此待处理的裂缝图像可能会存在边缘模糊、黑白杂点等问题，这在一定程度上会对裂缝目标的检测和识别产生影响，干扰实验结果的判断，因此需要对裂缝图像进行滤波去噪。本节将从均值滤波和中值滤波两方面来进行图像去噪的处理。均值滤波也称为邻域平均滤波，该方法假设待处理图像是由许多灰度值为常量的小区域组成的，并且相邻区域间存在较高的空间相关性，而噪声则显得相对独立。因此，通过将单个像素及其指定邻域内的所有像素按某种规则计算平均灰度值，再作为新图像中的对应像素值，可达到滤波去噪的目的，这一过程被称为均值滤波。邻域平均法属于非加权邻域平均范畴，是最常用的均值滤波操作。

图像边缘一般集中了图像的细节和高频信息，如果通过邻域平均法进行去噪，则往往会引起图像边缘的模糊，这也会对裂缝目标的检测带来不利影响。中值滤波是常用的非线性滤波方法，其主要思想是对像素邻域向量化取中值来进行滤波，具有运算简单、高效，能有效去除脉冲噪声的特点，在去噪的同时它也可以有效地保护图像的边缘细节信息。因此，本案例将采用中值滤波的方法来对裂缝图像进行去噪处理，其处理步骤如下。

(1)定位

在图像中移动模板，将模板中心与图像中的某个像素重合。

(2)计算

选择模板对应于图像的各像素灰度值，进行向量化，并将其进行排序。

(3)赋值

选择序列的中间值，作为输出赋予模板中心对应的像素。

如图2所示，根据中值滤波器形状和维数的不同，其模板有线形、十字形、方形、菱形等，不同形状的窗口也会产生不同的滤波效果。在对裂缝图像进行中值滤波处理时，其关键在于选择合适的模板形状和模板大小。

2.3 图像增强

路面裂缝图像的采集一般在室外进行，容易受到大气、光照、机械振动等因素的影响，采集到的裂缝图像可能存在整体偏暗或偏亮等问题，进而产生对比度较低的图像。此类图像的特点是灰度分布范围较小，集中在少量的灰度区间内，这也给后续的裂缝检测和识别带来了不利影响，因此需要对此类图像进行增强处理来提高对比度。直方图作为图像灰度级分布的统计表，能在一定程度上反映图像的对比度详情。图像的灰度直方图表示该图像所属灰度类型中不同灰度级像素出现的相对频率，并且直方图的横坐标表示灰度，纵坐标表示灰度出现的次数或概率。直方图均衡化利用灰度直方图进行图像对比度的调整，以达到增强图像视觉效果的目标。直方图均衡化的基本思想是通过某种变换，将原始图像的灰度直方图从集中于某个较小的灰度区间变成在更大灰度区间内均匀分布的形式，得到灰度级差式分布，从而达到增强图像整体对比度的目标。裂缝图像区域通常属于颜色较暗的灰度区间，背景区域则属于相对较亮的灰度区间。但在采集裂缝图像的过程中，往往会由于天气干扰、曝光不足等原因而造成图像整体偏暗，使裂缝区域与背景区域亮度特征相近而不易辨别，如图所示。从原始裂缝图像的灰度直方图可以看出，其灰度值分布主要集中在0100的低层次灰度区间。因此，为了提高裂缝与背景的对比度，需要将原图像的灰度值范围进行扩大，形成较为明显的灰度级差，进而增加裂缝图像的对比度。经过灰度直方图均衡化处理，裂缝图像的灰度范围扩到了0255，得到对比度增强后的裂缝图像，更加突出了裂缝与背景的差异程度。

2.4 图像二值化

灰度图像二值化是指通过约定一个灰度阈值来分割目标与背景，在阀值之内的像素于目标即记为1，其他则属于背景即记为0。在裂缝目标检测与识别的过程中，可以采用裂缝边缘、面积等特征来进行判别，也可以采用裂缝目标与周围背景的灰度差异值作为个判别依据，这就要求引入阈值来进行图像二值化处理。假设一幅灰度裂缝图像用f(x，y)表示，其中，(x，y)表示图像中像素的位置坐标，T为阈值，则阈值分割后的二值图像b(x，y)满足：

裂缝目标或背景区域的像素灰度通常是高度相关的，但裂缝目标与背景区域之间的灰度值则通常存在较大差异，一般包含明显的边缘等特征。因此，为了从更大程度上分割裂缝目标与背景，则需要对其进行灰度阈值分割选取合适的阈值。阈值计算方法根据其计算过程可以分为两种：全局阈值和基本自适应阈值，如下所述。

(1)全局阈值是最常见的阈值计算方法，它一般以图像的直方图或灰度空间分布为基础来确定一个阈值，进而实现灰度图像的二值化。特别是当图像的灰度直方图分布呈双峰时，全局阈值法可以明显地将目标和背景分量，得到较为理想的图像分割效果。但裂缝图像一般具有光照不均匀、噪声干扰等特点，其灰度直方图往往不会呈双峰分布，因此全局阈值分割方法效果较差。

(2)基本自适应阀值是一种比较基础的图像自适应分割方法，它一般以图像像素自身及其邻域灰度变化的特征为基础进行阈值分割，进而实现灰度图像的二值化。该方法充分考虑了每个像素邻域的特征，所以一般能更好地突出目标和背景的边界。

裂缝图像的背景在多数情况下比较固定，如路面、桥面、墙体等。但由于图像采集一般在室外进行，会受到拍摄条件、路面杂物等因素的影响，所以图像容易出现退化或噪声干扰。本案例通过分析裂缝图像目标和背景的特点，采用自定义和迭代法优化相结合的方法。

3 程序实现

根据裂缝图像的特点，在对其进行目标检测和识别之前，需要进行图像预处理，主要包括：直方图均衡化增强、中值滤波去噪、对比度增强、二值化处理、二值图像滤波等步骤。其中，在二值化过程中对阈值的确定选择自定义阈值法与迭代自适应法相结合的方式来计算；二值图像滤波主要是连通区域的面积滤波，通过去除小面积的杂点噪声来进行滤波去噪。裂缝图像经过预处理可以得到突出裂缝目标的二值图像，然后可以根据形态学区域特征来获取裂缝目标并进行检测识别。对于裂缝的形状识别可以通过计算图像中裂缝目标的外接矩形的长宽比来确定。

二、部分源代码

    clc;
    close all;
    clear all;
    %% 读取图片
    I=imread('1.jpg');
    %% 判断图像格式及调整图像尺寸、直方图均衡化
    if(ndims(I)==3)%ndims是求数组维数的函数
        img=rgb2gray(I);
        img=double(img);%转换为双精度
    else
        img=double(I);%转换为双精度
    end
    if(max(size(img))>1024)
        scale=1/8;  %图像调整为原图尺寸的1/8
        I1 = imresize(img, scale);%缩放处理
    else
        I1=img;
    end    
    [m,n]=size(I1);
    f_size=max(m,n);
    %% 去噪、锐化
    [Slog f2]=frequence_get(I1,f_size);
    %% 陷波滤波 去除周期性噪声 
    % % 用陷波滤波模板修改原图傅里叶频谱，滤除高频周期噪声
    % % 构造陷波滤波器模板
    Mask=Muban(Slog,f_size);
    vex=find(Mask==0);
    I_spec=I1;
    I_spec=im2uint8(mat2gray(I_spec));%矩阵归一化，把图像数据类型转换为无符号的8位整形
    I_spec(f_size,f_size)=0;
    Spectrum=fft2(uint8(I_spec));%2维离散傅里叶变换
    Spec=abs(Spectrum);   %幅度谱
    phi=angle(Spectrum);  %相位谱
    Spec=fftshift(Spec);
    Spec1=Spec;
    Spec1(vex)=0;
    Spec2=ifftshift(Spec1);
    CJ=Spec2.*exp(i*phi);  %重建傅里叶频谱
    Re_fft2=ifft2(CJ);     %反傅里叶变换，还原图像
    Restore=uint8(real(Re_fft2));
    %% 对还原图像维纳滤波
% %    维纳滤波     
    K=wiener2(Restore(1:m,1:n),[3 3]);%维纳滤波（wiener filtering) 一种基于最小均方误差准则、对平稳过程的最优估计器。这种滤波器的输出与期望输出之间的均方误差为最小，因此，它是一个最佳滤波系统。它可用于提取被平稳噪声所污染的信号。 
    [result,threshold1]=edge(K,'sobel',0.07);% 若原图没有进行均衡化，阈值用0.07，sobel算子边缘检测裂缝       
    [BW,thresh1]=bwfilter(K,result);%二值化滤波器去干扰
    figure,
    subplot(151);imshow(I1,[]);title('原始图像'); 
    subplot(152);imshow(Restore(1:m,1:n),[]);title('反傅里叶变换图像'); 
    subplot(153);imshow(K,[]);title('维纳滤波降噪');
    subplot(154);imshow(result,[]);title('裂纹检测'); 
    subplot(155);imshow(BW);title('二值化滤波器去干扰');
% %     形态学操作
    Se0=strel('line',3,0);%创建一个相对于邻域中心对称的线性结构元素
    Se90=strel('line',3,90);
    bw_dialte=imdilate(BW,[Se0 Se90]);%膨胀
    bw_fill=imfill(bw_dialte,'holes');%填充  
    bw1=bwmorph(bw_fill,'thin',Inf);%细化  
    bw2=bwmorph(bw1,'spur',3);  %去毛刺     
    bw3=bwareaopen(bw2,5);%移除小目标  
    figure,
    subplot(151);imshow(bw_dialte);title('膨胀');
    subplot(152);imshow(bw_fill);title('填充');
    subplot(153);imshow(bw1);title('细化');
    subplot(154);imshow(bw2);title('去毛刺');
    subplot(155);imshow(bw3);title('移除小目标');

三、运行结果