Matlab Notes_图像的边缘检测

function Td=EdgeDetec
tic
eI=imread('JupiterEarthComp.jpg');
eIg=rgb2gray(eI);                             % 读取rgb图像并转化为灰度图
sigma=1;                                      % 设定高斯标准差 sigma=1
filterExtent = ceil(4*sigma);                 % 根据高斯标准差计算滤波器长度
x = -filterExtent:filterExtent;               % 生成一维高斯核 gK=gaussKernel
c = 1/(sqrt(2*pi)*sigma);
gK = c * exp(-(x.^2)/(2*sigma^2));   
gK = gK/sum(gK);                              % 标准化
                                              % 对图像进行高斯滤波平滑 Smooth
BS1=imfilter(eIg,gK,'conv','replicate');      % 沿着X轴卷积
BS2=imfilter(BS1,gK','conv','replicate');     % 沿着Y轴卷积
DgK = gradient(gK);                           % 梯度 DgK=derivGaussKernel 
negVals = DgK < 0;
posVals = DgK > 0;
DgK(posVals) = DgK(posVals)/sum(DgK(posVals));
DgK(negVals) = DgK(negVals)/abs(sum(DgK(negVals)));
Bdx = double(imfilter(BS2, DgK, 'conv','replicate'));
Bdy = double(imfilter(BS2, DgK', 'conv','replicate'));
mag=hypot(Bdx,Bdy); 
magmax = max(mag(:));
if magmax>0
    magGrad = mag / magmax;                   % 梯度标准化
end
[m,n]=size(eIg);
counts=imhist(magGrad, 64);
highThresh = find(cumsum(counts) > 0.7*m*n, 1 ,'first' ) / 64;
lowThresh = 0.4*highThresh;
Bm = thinAndThreshold(Bdx, Bdy, magGrad, lowThresh, highThresh);
figure('Name','  Canny图像边缘检测')      
subplot(121);   imshow(eIg);     xlabel('原Gray图像');
subplot(122);  imshow(Bm);  xlabel('Canny边缘检测');
fprintf('\n高阈值TL: %f \n', highThresh);
fprintf('\n低阈值TH: %f \n', lowThresh);
imwrite(Bm,'JupiterEarthComp-Bm.png');
Td=toc;

 function H = thinAndThreshold(dx, dy, magGrad, lowThresh, highThresh)
%-% [Modified-01] E = cannyFindLocalMaxima(dx,dy,magGrad,lowThresh);  % 非极大值抑制
E = nonMaximalSupp(dx,dy,magGrad,lowThresh,highThresh);  % 非极大值抑制
if ~isempty(E)
    [rstrong,cstrong] = find(magGrad>highThresh & E);
    if ~isempty(rstrong) 
        H = bwselect(E, cstrong, rstrong, 8);            % 选定强边缘8连通目标
        H = bwmorph(H, 'thin', 1);                       % 边缘细化
    else
        H = false(size(E));
    end
else
    H = false(size(E));
end
end

 function idxLocalMax = cannyFLM(direction,ix,iy,mag)
%-% [Modified-02] function idxLocalMax = cannyFindLocalMaxima(direction,ix,iy,mag)
[m,n] = size(mag);
switch direction
 case 1
  idx = find((iy<=0 & ix>-iy)  | (iy>=0 & ix<-iy));
 case 2
  idx = find((ix>0 & -iy>=ix)  | (ix<0 & -iy<=ix));
 case 3
  idx = find((ix<=0 & ix>iy) | (ix>=0 & ix<iy));
 case 4
  idx = find((iy<0 & ix<=iy) | (iy>0 & ix>=iy));
end

if ~isempty(idx)
  v = mod(idx,m);
  extIdx = (v==1 | v==0 | idx<=m | (idx>(n-1)*m));
  idx(extIdx) = [];
end
 
ixv = ix(idx);  
iyv = iy(idx);   
gradmag = mag(idx);


switch direction
 case 1
  d = abs(iyv./ixv);
  gradmag1 = mag(idx+m).*(1-d) + mag(idx+m-1).*d; 
  gradmag2 = mag(idx-m).*(1-d) + mag(idx-m+1).*d; 
 case 2
  d = abs(ixv./iyv);
  gradmag1 = mag(idx-1).*(1-d) + mag(idx+m-1).*d; 
  gradmag2 = mag(idx+1).*(1-d) + mag(idx-m+1).*d; 
 case 3
  d = abs(ixv./iyv);
  gradmag1 = mag(idx-1).*(1-d) + mag(idx-m-1).*d; 
  gradmag2 = mag(idx+1).*(1-d) + mag(idx+m+1).*d; 
 case 4
  d = abs(iyv./ixv);
  gradmag1 = mag(idx-m).*(1-d) + mag(idx-m-1).*d; 
  gradmag2 = mag(idx+m).*(1-d) + mag(idx+m+1).*d; 
end
idxLocalMax = idx(gradmag>=gradmag1 & gradmag>=gradmag2); 

 function H = nonMaximalSupp(dx, dy, magGrad, lowThresh, highThresh) 
magmax = max(magGrad(:));
    if magmax > 0
        magGrad = magGrad / magmax;
    end
[row, col] = size(magGrad);
E = false(row, col);
idxStrong = [];
for dir = 1:4
    %-% [Modified-03] idxLocalMax = cannyFindLocalMaxima(dir,dx,dy,magGrad);
    idxLocalMax = cannyFLM(dir,dx,dy,magGrad);
    idxWeak = idxLocalMax(magGrad(idxLocalMax) > lowThresh);
    E(idxWeak)=1;
    idxStrong = [idxStrong; idxWeak(magGrad(idxWeak) > highThresh)]; %#ok<AGROW>
end
[m,n] = size(E);
if ~isempty(idxStrong)         % result is all zeros if idxStrong is empty
    rstrong = rem(idxStrong-1, m)+1;
    cstrong = floor((idxStrong-1)/m)+1;
    H = bwselect(E, cstrong, rstrong, 8);
else
    H = zeros(m, n);
end