自适应滤波器

自适应局部降噪滤波器

均值是计算平均值的区域上的平均灰度，方差是该区域上的图像对比度

$g(x,y)$噪声图像在$(x,y)$处的值
${\sigma }_{\eta }^2$ 为噪声的方差，为常数，需要通过估计得到
${\bar{z}}_{S_{xy}}$ 为局部平均灰度
${\sigma }_{S_{xy}}^2$ 为局部方差

$$\hat{f}(x,y)=g(x,y)-\frac{{\sigma }_{\eta }^2}{{\sigma }_{S_{xy}}^2}[g(x,y)-{\bar{z}}_{S_{xy}}]\text{\hspace{1em}(5.32)}$$

当${\sigma }_{\eta }^2>{\sigma }_{\eta }^2$ 时，将比率设置为1，这样做会使得滤波器是非线性的，但可阻止因缺少图像噪声方差的知识而产生无意义的结果（即负灰度级）。

若${\sigma }_{\eta }^2$估计值太低时，算法会因校正量小于就有的值而返回与原图像接近的图像。估计值太高会使得方差的比率在1.0处被水平，与正常情况相比，算法会更频繁地从图像中减去平均值。若允许为负值，且最后重新标定图像，则如前所述，结果 将损失图像的动态范围。

整个图像的相关值：
${\mu }_n={\sum }_{i=0}^{L-1}(r_i-m{)}^{n}p(r_i)$ 为灰度值r相对于其均值同的第n阶矩
$m={\sum }_{i=0}^{L-1}{r}_{i}p(r_i)$ 均值
${\sigma }^2={\mu }_2={\sum }_{i=0}^{L-1}(r_i-m{)}^{2}p(r_i)$ 方差

邻域内的相关值：
$m_{S_{xy}}={\sum }_{i=0}^{L-1}{r}_i{P}_{S_{xy}}(r_i)$ 均值
${\sigma }_{S_{xy}}^2={\mu }_2={\sum }_{i=0}^{L-1}(r_i-m_{S_{xy}}{)}^2{P}_{S_{xy}}(r_i)$ 方差

def calculate_sigma_eta(image):""":param image: input image:return: sigma eta of image"""hist, bins = np.histogram(image.flatten(), bins=256, range=[0, 256], density=True)r = bins[:-1]m = np.sum(r * hist)sigma = np.sum((r - m)**2 * hist)return sigmadef calculate_sigma_sxy(block):"""param block: input blockreturn: sigma sxy of block"""#==========公式法==========hist, bins = np.histogram(block.flatten(), bins=256, range=[0, 256], density=True)r = bins[:-1]m = (r * hist).sum()sigma = ((r - m)**2 * hist).sum()#===========等价于========
#     m = np.mean(block)
#     sigma = np.var(block)return sigma, m

def adaptive_local_denoise(image, kernel, sigma_eta=1):"""adaptive local denoising math: $$\hat{f}(x,y) = g(x,y) - \frac{\sigma_{\eta} ^ 2 }{\sigma_{S_{xy}}^2} \big[ g(x,y) - \bar z_{S_{xy}}\big]$$param: image:  input image for denoisingparam: kernel: input kernel, actually only use kernel shape, just want to keep the format as mean filterreturn: image after adaptive local denoising """epsilon = 1e-8height, width = image.shape[:2]m, n = kernel.shape[:2]padding_h = int((m -1)/2)padding_w = int((n -1)/2)# 这样的填充方式，可以奇数核或者偶数核都能正确填充image_pad = np.pad(image, ((padding_h, m - 1 - padding_h), \(padding_w, n - 1 - padding_w)), mode="edge")img_result = np.zeros(image.shape)for i in range(height):for j in range(width):block = image_pad[i:i + m, j:j + n]gxy = image[i, j]z_sxy = np.mean(block)sigma_sxy = np.var(block)rate = sigma_eta / (sigma_sxy + epsilon)if rate >= 1:rate = 1img_result[i, j] = gxy - rate * (gxy - z_sxy)return img_result

# 自适应局部降噪滤波器处理高斯噪声
img_ori = cv2.imread('DIP_Figures/DIP3E_Original_Images_CH05/Fig0513(a)(ckt_gaussian_var_1000_mean_0).tif', 0) #直接读为灰度图像kernel = np.ones([7, 7])
img_arithmentic_mean = arithmentic_mean(img_ori, kernel=kernel)
img_geometric_mean = geometric_mean(img_ori, kernel=kernel)
img_adaptive_local = adaptive_local_denoise(img_ori, kernel=kernel, sigma_eta=1000)plt.figure(figsize=(10, 10))plt.subplot(221), plt.imshow(img_ori, 'gray'), plt.title('With Gaussian noise'), plt.xticks([]),plt.yticks([])
plt.subplot(222), plt.imshow(img_arithmentic_mean, 'gray'), plt.title('Arithmentic mean'), plt.xticks([]),plt.yticks([])
plt.subplot(223), plt.imshow(img_geometric_mean, 'gray'), plt.title('Geomentric mean'), plt.xticks([]),plt.yticks([])
plt.subplot(224), plt.imshow(img_adaptive_local, 'gray'), plt.title('Adaptive local denoise'), plt.xticks([]),plt.yticks([])
plt.tight_layout()
plt.show()

自适应中值滤波器

$z_{\text{min}}是S_{xy}$邻域中的最小灰度值；$z_{\text{max}}是S_{xy}$邻域中的最大灰度值；$z_{\text{med}}是S_{xy}$邻域中的中值， $z_{xy}$是坐标$(x,y)$处的灰度值；$S_{\text{max}}$是$S_{xy}$允许的最大尺寸，是大于1的奇正整数。

自适应中值滤波算法在点$(x,y)$处使用两个处理层次，分别表示为层次$A$和层次$B$:

层次$A$:
若$z_{\text{min}}<z_{\text{med}}<z_{\text{max}}$ 则转到层次B
否则，增$S_{xy}$的尺寸，
若$S_{xy}\le S_{\text{max}}$, 则重复层次A
否则，输出$z_{\text{med}}$

层次$B$：
若$z_{\text{min}}<z_{xy}<z_{\text{max}}$，则输出$z_{xy}$
否则，输出 $z_{\text{med}}$

这算法有3个主要的目的：去除椒盐（冲激）噪声，平滑其他非常冲激噪声，减少失真（如目标边界的过度细化）。该算法统计上认为$z_{\text{min}}$和$z_{\text{max}}$是区域$S_{xy}$的“类冲激”噪声分量，即使它们不是图像中的最小像素和最大像素值。

能很好的处理椒盐噪声，但对高斯噪声不敏感

def adaptive_median_denoise(image, sxy=3, smax=7):"""adaptive median denoising param: image: input image for denoisingparam: sxy  : minimum kernel sizeparam: smax : maximum kernel sizereturn: image after adaptive median denoising """epsilon = 1e-8height, width = image.shape[:2]m, n = smax, smaxpadding_h = int((m -1)/2)padding_w = int((n -1)/2)# 这样的填充方式，可以奇数核或者偶数核都能正确填充image_pad = np.pad(image, ((padding_h, m - 1 - padding_h), \(padding_w, n - 1 - padding_w)), mode="edge")img_new = np.zeros(image.shape)for i in range(padding_h, height + padding_h):for j in range(padding_w, width + padding_w):sxy = 3     #每一轮都重置  k = int(sxy/2)block = image_pad[i-k:i+k+1, j-k:j+k+1]zxy = image[i - padding_h][j - padding_w]zmin = np.min(block)zmed = np.median(block)zmax = np.max(block)if zmin < zmed < zmax:if zmin < zxy < zmax:img_new[i - padding_h, j - padding_w] = zxyelse:img_new[i - padding_h, j - padding_w] = zmedelse:while True:sxy = sxy + 2k = int(sxy / 2)if zmin < zmed < zmax or sxy > smax:breakblock = image_pad[i-k:i+k+1, j-k:j+k+1]zmed = np.median(block)zmin = np.min(block)zmax = np.max(block)if zmin < zmed < zmax or sxy > smax:if zmin < zxy < zmax:img_new[i - padding_h, j - padding_w] = zxyelse:img_new[i - padding_h, j - padding_w] = zmedreturn img_new

# 自适中值滤波器处理椒盐噪声
img_ori = cv2.imread('DIP_Figures/DIP3E_Original_Images_CH05/Fig0514(a)(ckt_saltpep_prob_pt25).tif', 0) #直接读为灰度图像kernel = np.ones([7, 7])img_arithmentic_mean = median_filter(img_ori, kernel=kernel)
img_adaptive_median = adaptive_median_denoise(img_ori)plt.figure(figsize=(15, 10))
plt.subplot(231), plt.imshow(img_ori, 'gray'), plt.title('Salt pepper noise'), plt.xticks([]),plt.yticks([])
plt.subplot(232), plt.imshow(img_arithmentic_mean, 'gray'), plt.title('Arithmentic mean'), plt.xticks([]),plt.yticks([])
plt.subplot(233), plt.imshow(img_adaptive_median, 'gray'), plt.title('Adaptive median'), plt.xticks([]),plt.yticks([])
plt.tight_layout()
plt.show()

# 自适中值滤波器处理高斯噪声
img_ori = cv2.imread('DIP_Figures/DIP3E_Original_Images_CH05/Fig0513(a)(ckt_gaussian_var_1000_mean_0).tif', 0) #直接读为灰度图像kernel = np.ones([7, 7])
img_arithmentic_mean = median_filter(img_ori, kernel=kernel)
img_adaptive_median = adaptive_median_denoise(img_ori)plt.figure(figsize=(15, 10))plt.subplot(231), plt.imshow(img_ori, 'gray'), plt.title('With Gaussian noise'), plt.xticks([]),plt.yticks([])
plt.subplot(232), plt.imshow(img_arithmentic_mean, 'gray'), plt.title('Arithmentic mean'), plt.xticks([]),plt.yticks([])
plt.subplot(233), plt.imshow(img_adaptive_median, 'gray'), plt.title('Adaptive median'), plt.xticks([]),plt.yticks([])
plt.tight_layout()
plt.show()