自適應濾波器

自適應局部降噪濾波器

均值是計算平均值的區域上的平均灰度，方差是該區域上的圖像對比度

$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()