郑州做装饰的网站,平面设计很难学吗,网站建设基本代码,做网站要不要学ps直方图均衡化是一种在图像处理技术#xff0c;通过调整图像的直方图来增强图像的对比度。 本博客不利用opencv库#xff0c;仅利用numpy、matplotlib来实现直方图均衡、自适应直方图均衡、对比度受限自适应直方图均衡
直方图均衡
包括四个流程
计算图像RGB三通道的归一化直… 直方图均衡化是一种在图像处理技术通过调整图像的直方图来增强图像的对比度。 本博客不利用opencv库仅利用numpy、matplotlib来实现直方图均衡、自适应直方图均衡、对比度受限自适应直方图均衡
直方图均衡
包括四个流程
计算图像RGB三通道的归一化直方图计算变换函数 k k k为输入像素值输出像素值为 s T ( k ) sT(k) sT(k) s T ( k ) ( L − 1 ) ∑ j 0 k p r ( r j ) ( L − 1 ) ∑ j 0 k n j M N L − 1 M N ∑ j 0 k n j sT(k)(L-1)\sum_{j0}^kp_r(r_j)(L-1)\sum_{j0}^k\frac{n_j}{MN}\frac{L-1}{MN}\sum_{j0}^kn_j sT(k)(L−1)∑j0kpr(rj)(L−1)∑j0kMNnjMNL−1∑j0knj 其中L256表示256个像素值 p r ( r j ) p_r(r_j) pr(rj)表示像素值为 j j j的像素点概率MN表示像素点个数 n j n_j nj表示像素值为j的像素点个数 k ∈ [ 0 , L − 1 ] k\in[0, L-1] k∈[0,L−1]对于每个通道的每个像素点通过变换函数实现变换三个通道的所有像素点变换后合并起来组合成新的RGB图像 直方图均衡利用全图的直方图统计信息提升全图的对比度对局部图的对比度提升较小所以需要自适应直方图均衡。
自适应直方图均衡Adaptive Histogram Equalization, AHE
包括4个步骤
设置滑动窗口大小W (window size)以及步长S (stride)需要保证WS滑动窗口从图像左上角开始向右、向下滑动。滑动过程中窗口可能会滑出图像区域。为了避免该情况如果窗口超出图像区域把图像拉回来贴到最近的边缘。如下图蓝色窗口滑出图像区域将其拉回。 对滑动窗口内的像素值统计直方图进行变换。包括 在窗口W*W区域内统计直方图在窗口W*W区域内计算变换函数如果窗口位于图像边缘窗口在图像四个角、四条边上变换函数作用于窗口内所有像素点如果窗口位于图像内变换函数仅作用于窗口中心的S*S区域可以理解为利用窗口的统计信息优化窗口内局部区域的对比度 三个通道的所有像素点变换后合并起来组合成新的RGB图像 自适应直方图均衡对滑动窗口内的图像进行对比度增强提升图像局部的对比度。一方面局部对比度提升过大容易失真另一方面不同局部区域的直方图统计信息差距明显导致增强后的图像有明显的块状结构。
对比度受限自适应直方图均衡Contrast Limited Adaptive histogram equalization, CLAHE
具备6个步骤
设置滑动窗口W直方图阈值T。滑动窗口从图像左上角开始向右、向下滑动。滑动过程中窗口可能会滑出图像区域。为了避免该情况采用和直方图均衡不同的策略对原图做边缘的零填充。对滑动窗口内的像素值进行变换。变换逻辑为 在窗口W*W区域内统计直方图通过对直方图裁剪限制对比度。将直方图中超过直方图阈值T的部分均分到直方图的每个块中。计算每个窗口的变换函数通过对目标窗口的像素值进行插值缓解自适应直方图均衡带来的分块问题 目标窗口D内目标像素点位于4个子窗口a、b、c、d每个子窗口有4个上、下、左、右的领域窗口。例如子窗口a的领域窗口是A、B、C、D子窗口b的领域窗口是B、M、D、E。需要讨论一些边界情况如果窗口D在图像的左上角那么子窗口a的领域窗口只有D如果窗口D在紧挨图像上边缘那么子窗口b的领域窗口只有D和E…通过目标像素点和领域窗口的插值实现变换。例如对于在子窗口a的像素点先计算在A、B、C、D四个窗口变换函数下的变换值通过距离插值得到最终变换值。三个通道的所有像素点变换后合并起来组合成新的RGB图像
实验对比结果
对于自适应直方图均衡参数设置为W64, S20 对于对比度受限自适应直方图均衡参数设置为W40, T10 ori原图 1直方图均衡 2自适应直方图均衡 3对比度受限自适应直方图均衡
样例1 RGB三通道直方图比较
样例2 RGB三通道直方图比较
三种算法比较
直方图均衡自适应直方图均衡对比度受限自适应直方图均衡增强效果对比度较原图增强图像较连续局部区域增强不明显部分局部区域比较暗对比度较原图增强图像不连续有块状效应局部区域增强较明显对比度较原图增强图像连续无块状效应局部区域对比度也有提升直方图比较直方图变得均衡但首或尾可能值特别大直方图变得均衡但首或尾值依旧比其他位置大直方图均衡化程度没有前两者高运行时间快慢中等
整体代码
import matplotlib.pyplot as plot
import PIL
import numpy as np
from pylab import subplots_adjust
depth 3def show_histogram(img):展示三个通道直方图:param img: 图像numpy:return: void但是显示图像plot.subplot(3, 1, 1)plot.hist(img[:, :, 0].flatten(), 256)plot.subplot(3, 1, 2)plot.hist(img[:, :, 1].flatten(), 256)plot.subplot(3, 1, 3)plot.hist(img[:, :, 2].flatten(), 256)plot.show()def get_histogram(img, normed True) :计算img的直方图由于img有三个通道需要在三个通道上做直方图根据norm参数设置为0~1的概率值还是整数值:param img: 三维矩阵长*宽*通道数:param normd: 是否做归一化:return: 在rgb通道上的直方图his np.zeros((3, 256))total_num img.shape[0] * img.shape[1]for dep in range(depth) :for hei in range(img.shape[0]) :for wid in range(img.shape[1]) :his[dep][img[hei][wid][dep]] 1if normed True :his / total_numreturn hisdef histogramEqualization(img, save_pthno) :计算直方图均衡化之后的图像计算保存 直方图均衡化之后的图像:param img: 输入图像(numpy):param save_pth: 保存路径:return: 直方图均衡之后的图像(numpy)### 计算直方图his get_histogram(img, normed True)### 计算累加和his_culsum np.zeros((3, 256))for dep in range(depth) :sum_tmp 0.0for index in range(256) :sum_tmp his[dep][index]his_culsum[dep][index] sum_tmphis_cdf 255.0 * his_culsum### Q1 整数化处理(测试三种取整方式)# his_cdf np.around(his_cdf, 0)his_cdf np.floor(his_cdf)# img_hisequa np.zeros(img.shape)### Q2 保存图片前后的直方图不同意### Q3 numpy类型图片颜色显示不正确### 设置映射new_photo PIL.Image.new(RGB, (img.shape[1], img.shape[0]))img_hisequa np.array(new_photo)for hei in range(img.shape[0]) :for wid in range(img.shape[1]) :r img[hei][wid][0]g img[hei][wid][1]b img[hei][wid][2]img_hisequa[hei][wid][0] ,img_hisequa[hei][wid][1], img_hisequa[hei][wid][2] int(his_cdf[0][r]), int(his_cdf[1][g]), int(his_cdf[2][b])if save_pth ! no :img PIL.Image.fromarray(img_hisequa.astype(uint8)).convert(RGB)img.save(save_pth, quality 95)return img_hisequadef AHE(img, W, S, save_pth no) :对图像做自适应直方图均衡化对图像做自适应直方图均衡化保存到save_pth(当save_pthno时不处理):param img: 输入numpy图像:param W: 窗口边长参数W:param S: 影响区域大小边长/stride:param save_pth: 保存图像路径:return:assert W S### 按照S得到在水平方向和竖直方向滑动的次数如果不能取整加一img_height img.shape[0]img_width img.shape[1]hei_num img_height / Sif hei_num int(hei_num) :hei_num int(hei_num) 1else:hei_num int(hei_num)wid_num img_width / Sif wid_num int(wid_num) :wid_num int(wid_num) 1else :wid_num int(wid_num)print(hei_num, wid_num)# 处理边缘采用窗口为W影响区域为W处理最后一个区域采用同样的W大小可以忽视重叠部分new_photo PIL.Image.new(RGB, (img_width, img_height))img_hisequa np.array(new_photo)for hei_i in range(1, hei_num 1) :sta_hei S * (hei_i - 1)end_hei S * (hei_i)if end_hei img_height :sta_hei img_height - Wend_hei img_heightimg_hisequa[sta_hei:end_hei, 0: W, :] histogramEqualization(img[sta_hei:end_hei, 0: W, :])img_hisequa[sta_hei:end_hei, img_width - W: img_width, :] histogramEqualization(img[sta_hei:end_hei, img_width - W: img_width, :])for wid_i in range(1, wid_num 1) :sta_wid S * (wid_i - 1)end_wid S * (wid_i)if end_wid img_width :sta_wid img_width - Wend_wid img_widthimg_hisequa[0: W, sta_wid: end_wid, :] histogramEqualization(img[0: W, sta_wid: end_wid, :])img_hisequa[img_height - W: img_height, sta_wid: end_wid, :] histogramEqualization(img[img_height - W: img_height, sta_wid: end_wid, :])#其次处理内部### 类似的计算内部区域滑动次数有小数加一pad_num int((W - S) / 2)assert (W - S) % 2 0img_height_inside img_height - 2 * Wimg_width_inside img_width - 2 * Whei_num_inside img_height_inside / Sif hei_num_inside int(hei_num_inside) :hei_num_inside int(hei_num_inside) 1else :hei_num_inside int(hei_num_inside)wid_num_inside img_width_inside / Sif wid_num_inside int(wid_num_inside) :wid_num_inside int(wid_num_inside) 1else :wid_num_inside int(wid_num_inside)### 将窗口W按照S滑动最后一个区域也采用W窗口影响区域为S尽管会重叠但是效果相同base_left Wbase_right img_width - Wbase_up Wbase_down img_height - Wfor hei_i in range(1, hei_num_inside 1) :for wid_i in range(1, wid_num_inside 1) :up base_up S * (hei_i - 1)down base_up S * hei_iif down base_down :up base_down - Sdown base_downleft base_left S * (wid_i - 1)right base_left S * wid_iif right base_right :left base_right - Sright base_rightimg_tmp histogramEqualization(img[up - pad_num: down pad_num, left - pad_num: right pad_num, :])img_hisequa[up: down, left: right, :] img_tmp[pad_num: img_tmp.shape[0] - pad_num, pad_num: img_tmp.shape[1] - pad_num,:]if save_pth ! no :img_save PIL.Image.fromarray(img_hisequa.astype(uint8)).convert(RGB)img_save.save(save_pth, quality 95)return img_hisequadef clip(his, threshold 1.0) :裁剪直方图阈值threshold threshold_value * len(his) / total_num:param his: 直方图numpy一维:param threshold: 阈值 threshold_value * len(his) / total_numthreshold_value表示数量:return:裁剪之后的直方图total_num np.sum(his)# print(total_num)threshold_value total_num / len(his) * thresholdvalue_average 0for i in range(len(his)) :if his[i] threshold_value :value_average his[i] - threshold_valuevalue_average / len(his)his_clip np.zeros(his.shape)for i in range(len(his)) :if his[i] threshold_value :his_clip[i] threshold_value value_averageelse :his_clip[i] his[i] value_averagereturn his_clipdef get_his_clip_histogramEqualization(img, threshold 1.0) :对裁剪后的直方图做均衡在单通道上实现其中调用了clip函数。:param img: 输入图片单通道图片numpy:param threshold:阈值threshold threshold_value * len(his) / total_num:return:返回直方图的累计分布函数即映射函数his np.zeros(256)total_num img.shape[0] * img.shape[1]for hei in range(img.shape[0]):for wid in range(img.shape[1]):his[img[hei][wid]] 1his / total_numhis_clip clip(his, threshold threshold)his_culsum np.zeros(256)sum_tmp 0for index in range(256):sum_tmp his_clip[index]his_culsum[index] sum_tmphis_cdf 255.0 * his_culsumhis_cdf np.floor(his_cdf)return his_cdfdef CLAHE_single_channel(img, W, T):对单通道图像做CLAHE处理即直方图裁剪插值其中调用了get_his_clip_histogramEqualization来得到裁剪后的直方图:param img: input image:param W: window size W*W:param T: histogram threshold threshold_value * len(his) / sum(hist):return: 对单通道图像做了CLAHE处理之后的单通道图像numpy### clipimg_height img.shape[0]img_width img.shape[1]W_hei_num (img_height // W)W_wid_num (img_width // W)img np.pad(img, ((0,W - (img_height - W_hei_num * W)),(0,W - (img_width - W_wid_num * W))),constant,constant_values (0,0))W_hei_num img.shape[0] // WW_wid_num img.shape[1] // Wres_img np.zeros(img.shape)### 对单通道做直方图均衡#### 对整除的部分做处理his_map {}for w_hei_i in range(W_hei_num):for w_wid_i in range(W_wid_num):up w_hei_i * Wdown (w_hei_i 1) * Wleft w_wid_i * Wright (w_wid_i 1) * Whis_tmp get_his_clip_histogramEqualization(img[up: down, left: right], T)his_map[(w_hei_i, w_wid_i)] his_tmp### 插值for i_h in range(W_hei_num * W):for i_w in range(W_wid_num * W):### 计算左上角的块在W_hei_numW_wid_num块中的横纵坐标## 从0开始向下取整left_up_h int((i_h - W / 2) / W)left_up_w int((i_w - W / 2) / W)### 计算距离dis_x i_h - left_up_h * W - W / 2dis_x / Wdis_y i_w - left_up_w * W - W / 2dis_y / W### 左上角if left_up_h 0 and left_up_w 0:res_img[i_h][i_w] his_map[(left_up_h 1, left_up_w 1)][img[i_h, i_w]]### 右上角elif left_up_h 0 and left_up_w W_hei_num - 1:res_img[i_h][i_w] his_map[(left_up_h 1, left_up_w)][img[i_h, i_w]]### 右下角elif left_up_h W_hei_num - 1 and left_up_w W_hei_num - 1:res_img[i_h][i_w] his_map[(left_up_h, left_up_w)][img[i_h, i_w]]### 左下角elif left_up_h W_hei_num - 1 and left_up_w 0:res_img[i_h][i_w] his_map[(left_up_h, left_up_w 1)][img[i_h, i_w]]### 四个边界## 上边elif left_up_h 0 :left his_map[(0, left_up_w)][img[i_h, i_w]]right his_map[(0, left_up_w 1)][img[i_h, i_w]]res_img[i_h][i_w] (1 - dis_y) * left dis_y * right## 下边elif left_up_h W_hei_num - 1 :left his_map[(W_hei_num - 1, left_up_w)][img[i_h, i_w]]right his_map[(W_hei_num - 1, left_up_w 1)][img[i_h, i_w]]res_img[i_h][i_w] (1 - dis_y) * left dis_y * right## 左边elif left_up_w 0 :up his_map[(left_up_h, 0)][img[i_h, i_w]]down his_map[(left_up_h 1, 0)][img[i_h, i_w]]res_img[i_h][i_w] (1 - dis_x) * up dis_x * down## 右边elif left_up_w W_wid_num - 1 :up his_map[(left_up_h, W_wid_num - 1)][img[i_h, i_w]]down his_map[(left_up_h 1, W_wid_num - 1)][img[i_h, i_w]]res_img[i_h][i_w] (1 - dis_x) * up dis_x * downelse :left_up his_map[(left_up_h, left_up_w)][img[i_h][i_w]]left_down his_map[(left_up_h 1, left_up_w)][img[i_h][i_w]]right_up his_map[(left_up_h, left_up_w 1)][img[i_h][i_w]]right_down his_map[(left_up_h 1, left_up_w 1)][img[i_h][i_w]]res_img[i_h][i_w] (1 - dis_y) * ((1 - dis_x) * left_up dis_x * left_down) dis_y * ((1 - dis_x) * right_up dis_x * right_down)return res_img[:img_height, :img_width]### 处理右边下边的边界部分def CLAHE(img, W, T, save_pth no) ::param img: input image:param W: window size W*W:param T: histogram threshold threshold_value * len(his) / sum(hist):return:# save_pth astro_result3.jpgres_img img.copy()res_img[: ,:, 0] CLAHE_single_channel(img[:,:,0],W, T)res_img[:, :, 1] CLAHE_single_channel(img[:,:,1],W, T)res_img[:, :, 2] CLAHE_single_channel(img[:,:,2],W, T)if save_pth ! no :img_save PIL.Image.fromarray(res_img.astype(uint8)).convert(RGB)img_save.save(save_pth, quality95)return res_img### 以下三个函数为测试直方图均衡化自适应直方图均衡化限制对比度自适应直方图均衡化简单调用
def callhistogramEqualization() :img1_path ./rock.jpgimg1 PIL.Image.open(img1_path, r)img1_numpy np.array(img1)histogramEqualization(img1_numpy, rock_result1.jpg)show_histogram(img1_numpy)img1res_path rock_result1.jpgimg1res PIL.Image.open(img1res_path, r)img1res_numpy np.array(img1res)show_histogram(img1res_numpy)def callAHE() :img1_path ./rock.jpgimg1 PIL.Image.open(img1_path, r)img1_numpy np.array(img1)AHE(img1_numpy, W48, S16, save_pthrock_result2.jpg)show_histogram(img1_numpy)img1res_path rock_result2.jpgimg1res PIL.Image.open(img1res_path, r)img1res_numpy np.array(img1res)show_histogram(img1res_numpy)def callCLAHE() :img1_path ./rock.jpgimg1 PIL.Image.open(img1_path, r)img1_numpy np.array(img1)CLAHE(img1_numpy, W40, T 10, save_pthrock_result3.jpg)# CLAHE(img1_numpy, W40, A0, T10, save_pthastro_result3.jpg)show_histogram(img1_numpy)img1res_path rock_result3.jpgimg1res PIL.Image.open(img1res_path, r)img1res_numpy np.array(img1res)show_histogram(img1res_numpy)def getResult() :img1Pth ./astronaut.jpgimg1 PIL.Image.open(img1Pth, r)img1_numpy np.array(img1)img1_numpy_he histogramEqualization(img1_numpy)img1_numpy_ahe AHE(img1_numpy, W64, S20)img1_numpy_clahe CLAHE(img1_numpy, W40, T 10)fg plot.figure()# plot.subplots_adjust(wspace0, hspace0)ax1 plot.subplot(221)ax1.set_title(origin)ax1.axis(off)ax1.imshow(img1_numpy)ax2 plot.subplot(222)ax2.set_title(1)ax2.axis(off)ax2.imshow(img1_numpy_he)ax3 plot.subplot(223)ax3.set_title(2)ax3.axis(off)ax3.imshow(img1_numpy_ahe)ax4 plot.subplot(224)ax4.set_title(3)ax4.axis(off)ax4.imshow(img1_numpy_clahe)plot.tight_layout()plot.savefig(./result_1.jpg)plot.show()img2Pth ./rock.jpgimg2 PIL.Image.open(img2Pth, r)img2_numpy np.array(img2)img2_numpy_he histogramEqualization(img2_numpy)img2_numpy_ahe AHE(img2_numpy, W64, S20)img2_numpy_clahe CLAHE(img2_numpy, W40, T 10)plot.figure()ax1 plot.subplot(221)ax1.set_title(origin)ax1.axis(off)ax1.imshow(img2_numpy)ax2 plot.subplot(222)ax2.set_title(1)ax2.axis(off)ax2.imshow(img2_numpy_he)ax3 plot.subplot(223)ax3.set_title(2)ax3.axis(off)ax3.imshow(img2_numpy_ahe)ax4 plot.subplot(224)ax4.set_title(3)ax4.axis(off)ax4.imshow(img2_numpy_clahe)plot.tight_layout()plot.savefig(./result_2.jpg)plot.show()
if __name__ __main__ :# callhistogramEqualization()# callAHE()# callCLAHE()# img1res_path astor_official_result2.jpg# img1res PIL.Image.open(img1res_path, r)# img1res_numpy np.array(img1res)# show_histogram(img1res_numpy)getResult()
