我执行了一个简单的 MATLAB 实现，并得到了不错的结果。
MATLAB 代码：
I = imread('02_giorgos_1_f_M_30_830.tif');
I = im2double(uint8(I));
J = imdiffusefilt(I);
K = imtophat(J, ones(3));
figure;imshow(imadjust(K, stretchlim(K)));
结果：
我不知道你是否了解 MATLAB，但我使用了imdiffusefilt 的默认参数（相当于你代码中的anisodiff）。
默认的 MATLAB 参数等价于：
- 输入图像在 [0, 1] 范围内，而不是 [0, 255] 范围内。
-
niter=5（注意：您只使用了 1 次迭代，这还不够）。
kappa=0.1gamma=0.125- MATLAB 默认为 8 个邻居连接（不是
anisodiff 中使用的 4 个邻居）。
8 个邻居连接：
- 为了获得与 MATLAB 中相同的结果，我实现了一个 8 个邻居连通性各向异性扩散（基于 MATLAB 源代码）。
注意：使用 4 个邻居连接它可以工作，但结果不如使用 8 个邻居好。
显示输出图像：
- 为了正确显示输出图像，我使用了
imadjust(K, stretchlim(K))。
该命令拉伸输入图像的范围，使百分位 1 变为 0，百分位 99 变为 1（线性拉伸）。
还有一件事：
-
我没有使用固定阈值200，而是使用百分位 95 阈值：
t = np.percentile(first_tophat_img, 95)
ret, thresh1 = cv2.threshold(first_tophat_img, t, 255,
cv2.THRESH_BINARY)
这是代码（使用cv2.imshow 进行测试）：
import cv2
import numpy as np
import matplotlib.image as mpimg
import warnings
face_img = mpimg.imread('02_giorgos_1_f_M_30_830.tif')
def anisodiff8neighbors(img, niter=5, kappa=0.1, gamma=0.125):
""" See https://www.mathworks.com/help/images/ref/imdiffusefilt.html
Anisotropic diffusion filtering with 8 neighbors
Range of img is assumed to be [0, 1] (not [0, 255]).
"""
if img.ndim == 3:
m = "Only grayscale images allowed, converting to 2D matrix"
warnings.warn(m)
img = img.mean(2)
img = img.astype('float32')
imgout = img.copy()
for ii in range(niter):
# MATLAB source code is commented
#paddedImg = padarray(I, [1 1], 'replicate');
padded_img = np.pad(imgout, (1, 1), 'edge')
#diffImgNorth = paddedImg(1:end-1,2:end-1) - paddedImg(2:end,2:end-1);
#diffImgEast = paddedImg(2:end-1,2:end) - paddedImg(2:end-1,1:end-1);
#diffImgNorthWest = paddedImg(1:end-2,1:end-2) - I;
#diffImgNorthEast = paddedImg(1:end-2,3:end) - I;
#diffImgSouthWest = paddedImg(3:end,1:end-2) - I;
#diffImgSouthEast = paddedImg(3:end,3:end) - I;
diff_img_north = padded_img[0:-1, 1:-1] - padded_img[1:, 1:-1]
diff_img_east = padded_img[1:-1, 1:] - padded_img[1:-1, 0:-1]
diff_img_north_west = padded_img[0:-2, 0:-2] - imgout
diff_img_north_east = padded_img[0:-2, 2:] - imgout
diff_img_south_west = padded_img[2:, 0:-2] - imgout
diff_img_south_east = padded_img[2:, 2:] - imgout
#case 'exponential'
#conductCoeffNorth = exp(-(abs(diffImgNorth)/gradientThreshold).^2);
#conductCoeffEast = exp(-(abs(diffImgEast)/gradientThreshold).^2);
#conductCoeffNorthWest = exp(-(abs(diffImgNorthWest)/gradientThreshold).^2);
#conductCoeffNorthEast = exp(-(abs(diffImgNorthEast)/gradientThreshold).^2);
#conductCoeffSouthWest = exp(-(abs(diffImgSouthWest)/gradientThreshold).^2);
#conductCoeffSouthEast = exp(-(abs(diffImgSouthEast)/gradientThreshold).^2);
conduct_coeff_north = np.exp(-(np.abs(diff_img_north)/kappa)**2.0)
conduct_coeff_east = np.exp(-(np.abs(diff_img_east)/kappa)**2.0)
conduct_coeff_north_west = np.exp(-(np.abs(diff_img_north_west)/kappa)**2.0)
conduct_coeff_north_east = np.exp(-(np.abs(diff_img_north_east)/kappa)**2.0)
conduct_coeff_south_west = np.exp(-(np.abs(diff_img_south_west)/kappa)**2.0)
conduct_coeff_south_east = np.exp(-(np.abs(diff_img_south_east)/kappa)**2.0)
#fluxNorth = conductCoeffNorth .* diffImgNorth;
#fluxEast = conductCoeffEast .* diffImgEast;
#fluxNorthWest = conductCoeffNorthWest .* diffImgNorthWest;
#fluxNorthEast = conductCoeffNorthEast .* diffImgNorthEast;
#fluxSouthWest = conductCoeffSouthWest .* diffImgSouthWest;
#fluxSouthEast = conductCoeffSouthEast .* diffImgSouthEast;
flux_north = conduct_coeff_north * diff_img_north
flux_east = conduct_coeff_east * diff_img_east
flux_north_west = conduct_coeff_north_west * diff_img_north_west
flux_north_east = conduct_coeff_north_east * diff_img_north_east
flux_south_west = conduct_coeff_south_west * diff_img_south_west
flux_south_east = conduct_coeff_south_east * diff_img_south_east
#% Discrete PDE solution
#I = I + diffusionRate * (fluxNorth(1:end-1,:) - fluxNorth(2:end,:) + ...
# fluxEast(:,2:end) - fluxEast(:,1:end-1) + (1/(dd^2)).* fluxNorthWest + ...
# (1/(dd^2)).* fluxNorthEast + (1/(dd^2)).* fluxSouthWest + (1/(dd^2)).* fluxSouthEast);
imgout = imgout + gamma * (flux_north[0:-1,:] - flux_north[1:,:] +
flux_east[:,1:] - flux_east[:,0:-1] + 0.5*flux_north_west +
0.5*flux_north_east + 0.5*flux_south_west + 0.5*flux_south_east)
return imgout
#new_img = anisodiff(face_img, niter=1, kappa=20, gamma=0.1, step=(1., 1.), option=1)
face_img = face_img.astype(float) / 255;
#new_img = anisodiff(face_img, niter=5, kappa=0.1, gamma=0.125, step=(1., 1.), option=1)
new_img = anisodiff8neighbors(face_img, niter=5, kappa=0.1, gamma=0.125)
filterSize =(3, 3)
kernel = cv2.getStructuringElement(cv2.MORPH_RECT,
filterSize)
input_image = new_img
first_tophat_img = cv2.morphologyEx(input_image,
cv2.MORPH_TOPHAT,
kernel)
# Use percentile 95 (of image) as threshold instead of fixed threshold 200
t = np.percentile(first_tophat_img, 95)
ret, thresh1 = cv2.threshold(first_tophat_img, t, 255, cv2.THRESH_BINARY)
cv2.imshow('thresh1', thresh1)
filterSize =(3, 3)
kernel = cv2.getStructuringElement(cv2.MORPH_RECT,
filterSize)
second_tophat_img = cv2.morphologyEx(input_image,
cv2.MORPH_BLACKHAT,
kernel)
#ret, thresh1 = cv2.threshold(second_tophat_img, 200, 255, cv2.THRESH_BINARY)
# Use percentile 95 (of image) as threshold instead of fixed threshold 200
t = np.percentile(second_tophat_img, 95)
ret, thresh2 = cv2.threshold(second_tophat_img, t, 255, cv2.THRESH_BINARY)
cv2.imshow('thresh2', thresh2)
lo, hi = np.percentile(first_tophat_img, (1, 99))
first_tophat_img_stretched = (first_tophat_img.astype(float) - lo) / (hi-lo) # Apply linear "stretch" - lo goes to 0, and hi goes to 1
cv2.imshow('first_tophat_img_stretched', first_tophat_img_stretched)
cv2.waitKey()
cv2.destroyAllWindows()
结果：
thresh1:
thresh2:
first_tophat_img_stretched: