【裂缝识别】基于matlab组合BCOSFIRE过滤器墙体裂缝识别【含Matlab源码 321期】

一、简介

1 基于阔值分割的裂缝检测方法
基于阈值分割的裂缝检测方法是一种综合的裂缝检测算法，主要分为两步：首先对输入的图像进行预处理，目的是去除混凝土表面图像上的背景噪声，同时增强裂缝特征；其次对预处理的图像进行阈值分割，提取裂缝。
例如，（1）预处理中利用图像减操作去除微小的变化，如不均匀光照和阴影的干扰；（2）然后利用线性增强处理，增强线性结构，消除块状结构等噪声；（3）最后利用Otsu算法进行阈值分割，其阈值获取依据灰度直方图。
该算法能够有效去除阴影干扰，但是取决于预处理时图像减操作的中值滤波器大小、线性增强的参数大小及阈值选择的准确性，在实际应用中不容易获取。

2 基于形态学的裂缝检测方法
基于形态学的裂缝检测方法是数学形态学和曲率评估检测复杂环境中类似裂缝的模型。通过定义一个矩形或其他形状的结构元素，对混凝土表面图像进行腐蚀、膨胀、开运算和闭运算等操作，使得图像中不是裂缝的边缘信息尽量腐蚀掉，同时对裂缝边缘进行填充，再进行裂缝识别。
基于形态学的裂缝检测方法通常分三步：
（1）通过增强图像背景中黑色像素来提高图像对比度；
（2）采用数学形态学中描述的条件进行裂缝增强，即先使用高斯滤波去噪，然后进行阈值分割；
（3）采用线性结构元素去除被增强的噪声模型，利用较大的结构元素进行闭操作以去除比结构元素小的扭曲的裂缝部分。
该方法对复杂背景的混凝土表面图像进行裂缝检测时能够有效利用裂缝的形状属性、连通性、曲率等，在裂缝检测研究中应用较多。

3 基于区域生长的裂缝检测方法
区域生长是将单个像素或者像素点组成的区域发展成更大区域的过程。
首先选择种子像素点或者种子区域，然后对这些种子像素点进行区域增长是通过将与每个种子点或者区域有相似属性（如强度、灰度级及纹理等）的相邻像素点或者相邻区域合并到种子区域，不断迭代，直到处理完每个像素点。
使用区域生长算法进行裂缝检测时，首先对裂缝图像进行区域划分，寻找种子区域，然后对裂缝区域进行矢量化，根据裂缝区域间的相似属性与种子区域连接，不断迭代，直到所有裂缝区域连接完毕。
该算法能够实现断裂裂缝的连接，但计算量大，噪声区域的存在导致部分噪声区域被误连接，影响裂缝检测的精确率。
4 基于渗流模型的裂缝检测方法
渗流算法是一种可变的局部处理方法。
基于渗流模型的裂缝检测算法的基本思想是根据中心像素点的特征值判断其是否为裂缝上的点，若中心像素点属于裂缝像素，则渗流区域将线性生长，若中心像素点是背景像素，则渗流区域将向各方向生长。
渗流算法充分考虑了邻域像素点间灰度值的连续性，从而保证了裂缝的连续性，此外，该算法同时利用裂缝的灰度特征和形状特征，能够精确检测到不清晰的裂缝，但该算法处理时间长，不能很好地满足实时性要求。
5 基于小波变换的裂缝检测算法
小波变换是一种新的局部变换分析方法，能够提供一个随频率改变的窗口，即其形状、时间窗和频率窗都可以随时频的改变而改变。
采用小波变换进行图像处理时，首先将一幅固定的图像经过变换分解为一系列运动的图像，然后把函数看做是视野的移动。
当一幅图像发生伪运动时，其小波系数会发生震荡，震荡的幅度会给出函数的频率信息，该方法将图像变换到不同方向的子带来进行裂缝提取。
6 基于神经网络的裂缝检测方法
基于神经网络的裂缝检测算法是一种非线性的分类算法。
首先对图像进行方格采样，获取图像的方格影像，并对每个方格影像依次进行直方图均衡处理和二值化处理，然后提取方格特征，构成用于训练和分类的特征向量。
该方法对网络结构和训练方案的选择要求较高，否则会导致直接用于信息提取的结果鲁棒性差。此外，裂缝形状的多样性及宽度的不确定性导致训练过程复杂，运行效率低。

二、部分源代码

function [ ] = PavementCrackDelineation( )% NOTE: It requires a compiled mex-file of the fast implementation % of the max-blurring function.if ~exist('./COSFIRE/dilate')BeforeUsing();end% temporarydataset = 2;CRACK_IVC = 1;CRACK_PV14 = 2;if dataset == CRACK_IVCerror('Dataset not yet available.')return;elseif dataset == CRACK_PV14dname = 'CrackPV14';imagesdir = 'cracks14';        gtdir = 'cracks14_gt';prefix_gt = '';end% Binarization thresholdsthresholds = 0.01:0.01:0.99;nthresholds = numel(thresholds);%% Symmetric filter params and configurationx = 101; y = 101; % centerline1(:, :) = zeros(201);line1(:, x) = 1; %prototype line% Parameters determined in the paper % N.Strisciuglio, G. Azzopardi, N.Petkov, "Detection of curved lines % with B-COSFIRE filters: A case study on crack delineation", CAIP 2017sigma = 3.3;len = 14;sigma0 = 2;alpha = 1;% Symmetric filter paramssymmfilter = cell(1);symm_params = SystemConfig;% COSFIRE paramssymm_params.inputfilter.DoG.sigmalist = sigma;symm_params.COSFIRE.rholist = 0:2:len;symm_params.COSFIRE.sigma0 = sigma0 / 6;symm_params.COSFIRE.alpha = alpha / 6;% Orientationsnumoriens = 12;symm_params.invariance.rotation.psilist = 0:pi/numoriens:pi-pi/numoriens;% Configurationsymmfilter{1} = configureCOSFIRE(line1, round([y x]), symm_params);% Prepare the filter setfilterset(1) = symmfilter;%% APplication of B-COSFIRE for crack delineationfiles = rdir(['./data/' dname '/' imagesdir  '/*.bmp']);nfiles = size(files, 1);% Initialize result matrixnmetrics = 3;RESULTS = zeros(nfiles + 1, nmetrics, nthresholds);for n = 1:nfilesfprintf('Processing image %d of %d. ', n, nfiles);% Read imageimageInput = double(imread(files(n).name)) ./ 255;% Read groud truth[p, name, ext] = fileparts(files(n).name);gt = double(imread(['./data/' dname '/' gtdir '/' prefix_gt name '.bmp'])) ./ 255;imageInput = imcomplement(imageInput);% Pad input image to avoid border effectsNP = 50; imageInput = padarray(imageInput, [NP NP], 'replicate');% Filter responseinhibFactor = 0;tuple = computeTuples(imageInput, filterset);[response, rotations] = applyCOSFIRE_inhib(imageInput, filterset, inhibFactor, tuple);response = response{1};response = response(NP+1:end-NP, NP+1:end-NP);% Cropping out the central part (unpadding)rotations_final = zeros(size(response, 1), size(response, 2), size(rotations, 3));for j = 1:size(rotations, 3)rotations_final(:,:,j) = rotations(NP+1:end-NP, NP+1:end-NP, j);end% Evaluationfprintf(' Result evaluation...\n');for j = 1:nthresholds% Thinning and Histeresis thresholding (using different% thresholds). The threshold of the CAIP17 paper is th=49. Here% we compute the performance for differen thresholds anyway to% build and show the ROC curve.binImg = binarize(rotations_final, thresholds(j));binImg2 = bwmorph(binImg, 'close');binImg2 = bwmorph(binImg2,'skel',Inf);figure(j)imshow(binImg2);% Compute the result metrics for a tolerance of d = 2 (as in% the paper)[cpt2, crt2, F2] = evaluate(binImg2, gt, 2);%[cpt3, crt3, F3] = evaluate(binImg2, gt, 3);RESULTS(n, :, j) = [cpt2, crt2, F2];endend% Average Resultsavg_results = reshape(mean(RESULTS(1:nfiles, :, :)), nmetrics, nthresholds)';[M, idx] = max(avg_results(:,3));fprintf('\nResults of the CAIP17 paper\n');fprintf('Pr: %.3f, Re: %.3f, F: %.3f\n', avg_results(idx, 1), avg_results(idx, 2), avg_results(idx, 3));PrintROCcurve(avg_results);
endfunction [cpt, crt, F] = evaluate(binImg, gt, d)A = zeros(d*2+1, d*2+1); A(d+1,d+1) = 1; B = bwdist(A) <= d;[m, n] = size(binImg);%binImg = bwmorph(binImg,'skel',Inf);gt = padarray(gt, [d d], 0);%binImg = padarray(binImg, [d d], 0);Lr = 0;bad = zeros(size(binImg));for x = 1:mfor y = 1:n%if gt(x, y) == 1 if binImg(x,y) == 1%patch = binImg(x:x+2*d, y:y+2*d); %patch = gt(x+d-d:x+d+d, y-d+d:y+d+d);s = sum(patch(:) .* B(:));if s > 0Lr = Lr + 1;elsebad(x,y) = 1;endendendendLgt = sum(gt(:));Ln = sum(binImg(:));cpt = min(1, Lr / Lgt);crt = min(1, Lr / Ln);F = 2 * cpt * crt / (cpt + crt);
endfunction [binarymap] = binarize(rotoutput1, highthresh)%%%%%%%%%%%%%%%%% BEGIN BINARIZATION %%%%%%%%%%%%%%%%%%% compute thinningorienslist = 0:pi/12:pi-pi/12;[viewResult, oriensMatrix] = calc_viewimage(rotoutput1,1:numel(orienslist), orienslist);thinning = calc_thinning(viewResult, oriensMatrix, 1);%figure; imagesc(thinning);% % % Choose high threshold of hysteresis thresholding% if nargin == 4%     bins = 64;p = 0.05; %Keep the strongest 10% of the pixels in the resulting thinned image%     f = find(thinning > 0);%     counts = imhist(thinning(f),bins);%     highthresh = find(cumsum(counts) > (1-p)*length(f),1,'first') / bins;% end% binarymap = calc_hysteresis(thinning, 1, 0.5*highthresh*max(thinning(:)), highthresh*max(thinning(:)));%figure;imagesc(binarymap);colormap gray; axis image;% show binarized image% figure;% subplot(1,2,1);imagesc(img);axis off;axis image;colormap(gray);% subplot(1,2,2);imagesc(imcomplement(binarymap));axis off;axis image;colormap(gray);%%%%%%%%%%%%%%%%% END BINARIZATION %%%%%%%%%%%%%%%%%%%%
endfunction [] = PrintROCcurve(avg_results)pr = avg_results(:,1);re = avg_results(:,2);figure;linewidth = 3;plot(re, pr, 'linewidth',linewidth,'color',[0.25 0.25 0.25],'markersize',10);set(gca,'YGrid','off');set(gca,'XGrid','off');set(gca,'XTick',0:.1:1)set(gca,'XTickLabel',0:.1:1)axis square;% Plot other methods resultshold on;% 1: Zou14 - data set% 2: CrackTree% 3: FoSAPr = [0.872 0.821 0.845;0.842 0.625 0.733;0.846 0.885 0.897; 0.793 0.753 0.756; 0.949 0.845 0.860;0.671 0.780 0.836; 0.960 0.698 0.716; 0.846 0.696 0.749; 0.767 0.722 0.779; 0.833 0.927 0.811;0.833 0.839 0.792; 0.997 0.847 0.868; 0.499 0.775 0.696; 0.848 0.948 0.925]; Re = [0.965 0.691 0.628 ;0.904 0.605 0.568 ;0.905 0.713 0.612 ;0.903 0.776 0.691 ;0.939 0.600 0.577;0.843 0.649 0.647 ;0.915 0.605 0.552 ;0.929 0.668 0.654 ;0.996 0.669 0.636 ;0.961 0.860 0.805;0.993 0.967 0.937 ;0.823 0.923 0.805 ;0.890 0.706 0.663 ;0.988 0.985 0.880];F = [0.916 0.751 0.721 ;0.872 0.614 0.640 ;0.874 0.790 0.728 ;0.845 0.764 0.722 ;0.944 0.700 0.691;0.747 0.708 0.729 ;0.937 0.648 0.623 ;0.886 0.682 0.698 ;0.867 0.695 0.700 ;0.892 0.892 0.808;0.906 0.898 0.858 ;0.893 0.883 0.835 ;0.639 0.739 0.679 ;0.913 0.966 0.901];avgPr = mean(Pr);avgRe = mean(Re);avgF = mean(F);linewidth = 2;plot(avgRe(1), avgPr(1), 'o', 'linewidth',linewidth,'color',[0.25 0.25 0.25],'markersize',10);plot(avgRe(2), avgPr(2), 's', 'linewidth',linewidth,'color',[0.25 0.25 0.25],'markersize',10);plot(avgRe(3), avgPr(3), 'd', 'linewidth',linewidth,'color',[0.25 0.25 0.25],'markersize',10);xlabel('Recall');ylabel('Precision');legend({'COSFIRE', 'Zou14', 'CrackTree', 'FoSA'}, 'Location', 'southwest');title('ROC curve');fprintf('\nResults of Zou et al.\n');fprintf('Pr: %.3f, Re: %.3f, F: %.3f\n', avgPr(1), avgRe(1), avgF(1));fprintf('\nResults of CrackTree.\n');fprintf('Pr: %.3f, Re: %.3f, F: %.3f\n', avgPr(2), avgRe(2), avgF(2));fprintf('\nResults of FoSA.\n');fprintf('Pr: %.3f, Re: %.3f, F: %.3f\n', avgPr(3), avgRe(3), avgF(3));end

三、运行结果

四、matlab版本及参考文献

1 matlab版本
2014a

2 参考文献
[1] 蔡利梅.MATLAB图像处理——理论、算法与实例分析[M].清华大学出版社，2020.
[2]杨丹,赵海滨,龙哲.MATLAB图像处理实例详解[M].清华大学出版社，2013.
[3]周品.MATLAB图像处理与图形用户界面设计[M].清华大学出版社，2013.
[4]刘成龙.精通MATLAB图像处理[M].清华大学出版社，2015.