close?all;
clear;
clc;
[filenamepathname]?=?uigetfile（‘*.*‘‘選擇圖片‘‘E:\裂紋\可分類裂紋\31.png‘）;
if?isequal（filename0）||isequal（pathname0）
????return;
else
????filefullpath?=?[pathnamefilename];
end
%%?傅里葉變換低通濾波，傅里葉反變換。
d?=?30;?
bw?=?imread（filefullpath）;
bw?=?rgb2gray（bw）;
figure（1）;
subplot（221）;imshow（bw）;

F?=?fftshift（fft2（double（bw）））;

F1?=?abs（F）;
T?=?log（F1+1）;
subplot（222）;
imshow（T[]）

[mn]?=?size（bw）;
m_mid?=?fix（m/2）;
n_mid?=?fix（n/2）;
F_lpf?=?zeros（mn）;
h?=?zeros（mn）;
for?i?=?1:m
????for?j?=?1:n
????????distance?=?sqrt（（i-m_mid）^2+（j-n_mid）^2）;
????????if?distance?>=?d
????????????h（ij）?=?0;
????????else
????????????h（ij）?=?1;
????????end
????????F_lpf（ij）?=?h（ij）*F（ij）;
????end
end
bw?=?ifftshift（F_lpf）;
bw?=?uint8（real（ifft2（bw）））;
subplot（223）;imshow（bw）;title（‘低通濾波圖像‘）;

%?%%?直方圖均衡
%?bw?=?histeq（bw128）;
%?subplot（224）;imshow（bw）;title（‘直方圖均衡化后的圖像‘）;
%%?各向異性擴散
diff_im=bw;???
num_iter?=5;
delta_t?=?1/4;
kappa?=4;
dx?=?1;
dy?=?1;
dd?=?sqrt（2）;
%這是4個方向
hN?=?[0?1?0;?0?-1?0;?0?0?0];
hS?=?[0?0?0;?0?-1?0;?0?1?0];
hE?=?[0?0?0;?0?-1?1;?0?0?0];
hW?=?[0?0?0;?1?-1?0;?0?0?0];
%非正方向
%?hNE?=?[0?0?1;?0?-1?0;?0?0?0];
%?hSE?=?[0?0?0;?0?-1?0;?0?0?1];
%?hSW?=?[0?0?0;?0?-1?0;?1?0?0];
%?hNW?=?[1?0?0;?0?-1?0;?0?0?0];
f=diff_im/255;???%將圖像歸一化
%?Anisotropic?diffusion.
for?t?=?1:num_iter
????????nablaN?=?imfilter（diff_imhN‘conv‘）;
????????nablaS?=?imfilter（diff_imhS‘conv‘）;???
????????nablaW?=?imfilter（diff_imhW‘conv‘）;
????????nablaE?=?imfilter（diff_imhE‘conv‘）;???
%?????????nablaNE?=?imfilter（diff_imhNE‘conv‘）;
%?????????nablaSE?=?imfilter（diff_imhSE‘conv‘）;???
%?????????nablaSW?=?imfilter（diff_imhSW‘conv‘）;
%?????????nablaNW?=?imfilter（diff_imhNW‘conv‘）;?
????????????cN?=?1./（1?+?f.*（kappa./nablaN）.^2）;
????????????cS?=?1./（1?+?f.*（kappa./nablaS）.^2）;
????????????cW?=?1./（1?+?f.*（kappa./nablaW）.^2）;
????????????cE?=?1./（1?+?f.*（kappa./nablaE）.^2）;
%?????????????cNE?=?1./（1?+?f.*（kappa./nablaN）.^2）;
%?????????????cSE?=?1./（1?+f.*（kappa./nablaN）.^2）;
%?????????????cSW?=?1./（1?+f.*（kappa./nablaN）.^2）;
%?????????????cNW?=?1./（1?+?f.*（kappa./nablaN）.^2）;
%?????????diff_im?=?diff_im?+?...
%???????????????????delta_t*（...
%???????????????????（1/（dy^2））*cN.*nablaN?+?（1/（dy^2））*cS.*nablaS?+?...
%???????????????????（1/（dx^2））*cW.*nablaW?+?（1/（dx^2））*cE.*nablaE?+?...
%???????????????????（1/（dd^2））*cNE.*nablaNE?+?（1/（dd^2））*cSE.*nablaSE?+?...
%???????????????????（1/（dd^2））*cSW.*nablaSW?+?（1/（dd^2））*cNW.*nablaNW?）;
??????diff_im?=?diff_im?+?...
??????????????????delta_t*（...
??????????????????cN.*nablaN?+?cS.*nablaS?+?...
??????????????????cW.*nablaW?+?cE.*nablaE?）;
????????fprintf（‘\rIteration?%d\n‘t）;
end
figure（2）imshow（diff_im[]）;
I=abs（diff_im-bw）;
figure（3）imshow（I[]）;

mean_I=mean（mean（I））;
std_I=std2（I）;
C=1.5;
threshold=mean_I+C*std_I;

%?threshold=5;
A=size（I）;
for?i=1:A（1）
????for?j=1:A（2）
????????if?I（ij）