16.直方图均衡化

数字图像处理(17): 直方图均衡化处理

简介

  直方图均衡化是一种简单有效的图像增强技术，通过改变图像的直方图来改变图像中各像素的灰度，主要用于增强动态范围偏小的图像的对比度。当原始图像的灰度分布较为集中的时候，可能造成图像不够清晰，图像会过曝或者过暗。采用直方图均衡化，可以将原本较为集中的灰度分布变换成较为均匀的形式，增加了图像的灰度之间的动态范围，从而增加了图像的对比度。直方图均衡化的主要思想就是对图像灰度较为集中的地方进行拉伸展宽，对图像中像素较少的灰度值进行归并。

直方图均衡化理论基础

  严格来说，图像的灰度直方图是一个一维的离散函数，可以写成：
$$h(k)=nkk=0,1,2,3,...,L-1\text{\hspace{1em}(1)}$$
  公式中$nk$是图像$f(x,y)$中灰度级为$k$的像素的个数。直方图的每一列的高度对应$nk$。直方图提供了原图中各种灰度值分布的情况，在直方图的基础上，进一步定义归一化的直方图为灰度级出现的相对频率$P_r(k)$，即：
$$P_r(k)=\frac{nk}{N}\text{\hspace{1em}(2)}$$
  式中，$N$表示图像$f(x,y)的像素的总数$，$nk$是图像$f(x,y)$中灰度级为$k$的像素的个数。
  得到直方图的相对频率后就可以按照需求对直方图进行拉伸和映射。
  为了方便讨论，以$r$和$s$分别表示归一化了的原图像灰度和经直方图均衡化后的灰度，因为归一化了，所以r和s的取值范围都在0到1之间。所谓的直方图均衡化，就是根据直方图对像素点的灰度值进行变换，属于点操作范围。换而言之，即：已知r，求相应的s。
  在[0:1]区间内任意一个r，经变换函数$T(r)$都可以产生一个对应的s，且
$$s=T(r)\text{\hspace{1em}(3)}$$
  式中，$T(r)$应当满足以下两个条件：
  （1）在$0\le r\le 1$内，$T(r)$为单调递增函数；
  （1）在$0\le r\le 1$内有 $0\le T(r)\le 1$；

直方图均衡化的实际运算过程

  （1）首先遍历图像，统计出每个灰度值各有多少个像素。
  （2）统计每个灰度值的像素各占一整幅图的多少百分比。以及累计直方图。
  （3）将累计直方图百分比进行区间映射（也就是乘以像素级255）得到像素级的映射关系。
  （4）根据映射关系对像素进行重新分配。

MATLAB实现

clear;
clc;
close all;

a = imread('../img/pre_hist_rq.bmp');

imshow(a);
[row,col] = size(a);

hist_lut = zeros(1,256);
lut_add = zeros(1,256);

result = zeros(row,col);

%Statistical gray value   计算一副图像的灰度值分布数量
for i = 1:row
    for j = 1:col
        hist_lut(a(i,j)+1) = hist_lut(a(i,j)+1)+1;
    end
end

%Cumulative histogram     统计累计直方图
for i=1:256
    if(i==1)
        lut_add(i) = hist_lut(i);
    else 
        lut_add(i) = lut_add(i - 1) + hist_lut(i);
    end
end 

%累计直方图归一化与重映射
for i = 1:row
    for j=1:col
        result(i,j) = lut_add(a(i,j)+1) * 255/(row*col);
    end
end

matlab_hist_eq = uint8(floor(result));

hist_lut1 = zeros(1,256);
%图像灰度重映射
for i = 1:row
    for j = 1:col
        hist_lut1(matlab_hist_eq(i,j)+1) = hist_lut1(matlab_hist_eq(i,j)+1) + 1;
    end 
end

subplot(2,2,1);
imshow(uint8(a));

subplot(2,2,2);
bar(hist_lut);

subplot(2,2,3);
imshow(uint8(result));

subplot(2,2,4);
bar(hist_lut1);

MATLAB实现如何转化为FPGA实现

  MATLAB中首先定义了两个数组，用来存储灰度直方图和累计直方图的数据。

hist_lut = zeros(1,256);   %灰度直方图
lut_add = zeros(1,256);    %灰度累计直方图

  FPGA中同样定义一个二维数组，用来存储灰度直方图。并将其初始化为0。以及一个RAM用来存储累计直方图。

reg  	[31:0] 	hist_lut 	[255:0] 	; 	 	//存放图像的灰度分布情况
integer 	i;

initial begin
	for(i=0;i<256;i=i+1)begin
		hist_lut[i] = 32'd0;
	end
end

//存储灰度直方图的ram
hist_ram u1_hist_ram (
	.clka		(clk			),  	
	.wea		(wea			),   	
	.addra		(addra			), 	
	.dina		(lut_add		),  	
	.douta		(				), 	

	.clkb		(clk			),  	
	.web		(1'b0			),   	
	.addrb		(pre_data		), 	
	.dinb		(32'd0			),  	
	.doutb		(hist_sum		)  	
);

  然后MATLAB开始统计一整幅图像的灰度直方图，利用图像的灰度值作为数组的索引，由于MATLAB下标是1到256，而灰度值是0到255，所以需要将图像的灰度值$a(i,j)+1$才能对应上索引。

%Statistical gray value   计算一副图像的灰度值分布数量
for i = 1:row
    for j = 1:col
        hist_lut(a(i,j)+1) = hist_lut(a(i,j)+1)+1;
    end
end

  这段MATLAB代码在FPGA中也很好实现。

//当数据到来时相应的灰度值数量+1  当累计直方图计算完毕后依次清除数据
always @(posedge clk)
	if(clc_de)
		hist_lut[clc_cnt] <= 32'd0;
	else if(pre_de==1)
		hist_lut[pre_data] <= hist_lut[pre_data] + 1;

  除了这段代码之外，由于我们还需要知道一副图像上面时候完结的，方便我们拿到完整的图像直方图后进行下一步的运算。

always @(posedge clk)
    if(rst_n == 0)
        hcnt    <=  0;
    else if(hcnt == COL - 1 && pre_de == 1'b1)
        hcnt    <=  0;
    else if(pre_de == 1'b1)
        hcnt    <=  hcnt+1;
    else 
        hcnt    <=  hcnt;

always @(posedge clk)
    if(rst_n == 0)
        vcnt    <=  0;
    else if(hcnt == COL - 1 && vcnt == ROW - 1 && pre_de == 1'b1)
        vcnt    <=  0;
    else if(hcnt == COL - 1)
        vcnt    <=  vcnt + 1;
    else 
        vcnt    <=  vcnt    ;

  MATLAB拿到灰度直方图后开始进行累计直方图的运算。

%Cumulative histogram     计算一副图像灰度值占比   
for i=1:256
    if(i==1)
        lut_add(i) = hist_lut(i);
    else 
        lut_add(i) = lut_add(i - 1) + hist_lut(i);
    end
end

  FPGA中也会进行累计直方图的运算，但是这个运算会消耗256个时钟周期，在一帧图像结束后，FPGA会拉高一个累计直方图的运算使能，并进行累计直方图的运算计数。依次将结果写入RAM。这里分五个always语句块来实现。
  （1）第一个语句块用来控制累计直方图的计算使能信号，在一幅图的灰度直方图计算完毕后拉高，累计直方图计算完毕后拉低

//当一帧数据结束后，开始统计直方图使能
always @(posedge clk)begin
	if(!rst_n)begin
		inter_de 	<= 	0 ;
	end
	else if(inter_cnt == 8'd255)begin
		inter_de 	<= 	0;
	end
	else if(hcnt == COL -1 && vcnt == ROW - 1)begin
		inter_de 	<= 	1'b1;
	end
	else begin
		inter_de 	<= 	inter_de ;
	end
end

  （2）第二个always语句块用来对累计直方图的计算进行计数，由于灰度值有0到255共256个，因此这个计数为0到255。

//依次统计直方图的灰度使能
always @(posedge clk)begin
	if(!rst_n)begin
		inter_cnt 	<= 0;
	end
	else if(inter_cnt == 8'd255)begin
		inter_cnt 	<= 0;
	end
	else if(inter_de == 1)begin
		inter_cnt 	<= inter_cnt + 1'b1;
	end
	else begin
		inter_cnt 	<= inter_cnt;
	end
end

  （3）第三个always语句块就根据前面两个信号做累计直方图。

//依次统计直方图灰度
always @(posedge clk)
	if(rst_n == 0)
		lut_add 	<= 	0;
	else if(inter_de)
		lut_add 	<= 	lut_add 	+ 	hist_lut[inter_cnt];
	else 
		lut_add 	<= 	'd0;

  （4）第四个always语句块就拉高ram的写入使能

//ram写使能拉高
always @(posedge clk )
	if(rst_n == 0)
		wea 	<= 	1'b0;
	else 
		wea 	<= 	inter_de;

  （5）第五个always语句块就控制写入地址

//ram写地址
always @(posedge clk)
	if(rst_n == 0)
		addra 	<= 'd0 ;
	else if(addra == 8'd255)
		addra 	<= 	'd0;
	else if(wea == 1'd1)
		addra 	<= 	addra + 1'b1;
	else 
		addra 	<= 	addra ;

  当直方图统计完成后就需要清除灰度直方图的数据来准备统计下一副图像的灰度直方图了。

//统计完成后清除灰度直方图数据
always @(posedge clk)
	if(rst_n == 0)
		clc_de 	<= 	1'b0;
	else 
		clc_de 	<= 	inter_de ;

always @(posedge clk)
	if(rst_n == 0)
		clc_cnt 	<= 0;
	else 
		clc_cnt 	<= 	inter_cnt ;

  下面就是根据累计直方图来计算直方图均衡化后的像素值了。先看MATLAB公式。

%直方图均衡化 将灰度值占比的0~1的256个值映射成为0~255
%先求出（0~255）像素的占比  然后将占比乘以255得到原始与均衡后的映射结果
for i = 1:row
    for j=1:col
        result(i,j) = lut_add(a(i,j)+1)*255/(row*col);
    end
end

  在FPGA中，至此，累计直方图以及计算完毕，接下来就要根据累计直方图来计算新的灰度图了。这里我们来看一下这个双端口RAM。RAM的B端口是一个只读端口，读的地址是根据输入的像素来的，由于RAM内存储的是累计直方图，于是RAM的输出是根据输入像素的灰度值所在的累计直方图的数据。这里先看我们定义的一个参数coe，他对应MATLAB中的放大了$2^{25}$次方的一个乘数，做了一个定浮点数的转化。然后将灰度值对应的累计直方图，也就是RAM的B端口根据原图像灰度值索引的累计直方图做一个乘法。计算出直方图均衡化的一个对应灰度值。注意，这里由于我们将这个乘数在定义时扩大了$2^{25}$次方，于是在计算出结果后我们也需要将数据进行位截取来再做输出。

reg  [0:1]	pre_de_r 	;

reg  	[63:0] 	mult 	;

localparam 	[63:0]  	coe 	= 	(2**25)*255/(COL*ROW) ; //累计直方图归一化以及重映射需要乘以的数据   乘数是累计直方图的灰度值

always @(posedge clk)
	if(rst_n == 0)
		pre_de_r 	<= 	0;
	else 
		pre_de_r 	<= 	{pre_de_r[0],pre_de};

//图像的数据进行重映射
always @(posedge clk)
	if(rst_n == 0)
		mult 	<= 	0;
	else if(pre_de_r[0] == 1'b1)
		mult 	<= 	hist_sum * coe;
	else 
		mult 	<= 	'd0;

always @(posedge clk)
	if(rst_n == 0)
		post_de 	<= 	1'b0;
	else 
		post_de 	<= 	pre_de_r[1];


//图像数据的位截取
always @(posedge clk)
	if(rst_n == 0)
		post_data 	<= 	1'b0;
	else if(pre_de_r[1] == 1'b1)
		post_data 	<= 	mult[25+7:25];
	else 
		post_data 	<= 	'd0;

  至此，MATLAB转FPGA算法的转化已经全部完成。后面贴上完整的FPGA算法模块。

FPGA实现

module hist_eq#(
	parameter 	ROW 	= 	1920 		,
	parameter 	COL 	= 	1080 
)(
	input 	wire  			clk  		,

	input 	wire  			pre_vs 		,
	input 	wire  			pre_de 		,
	input 	wire  	[7:0] 	pre_data	,

	output 	wire   			post_vs 	,
	output 	wire  			post_de 	,
	output 	wire  	[7:0] 	post_data 
);

reg  	[31:0] 	hist_lut 	[255:0] 	; 	 	//存放图像的灰度分布情况

integer 	i;

initial begin
	for(i=0;i<256;i++)begin
		hist_lut[i] = 32'd0;
	end
end

reg  			pre_vs_r 	;

wire  			vs_pose 	;
wire  			ext_vs_pose ;

reg  			rst_n 		;

reg  	[15:0] 	hcnt 		;
reg  	[15:0] 	vcnt 		;

reg  			inter_de 	;
reg  	[7:0] 	inter_cnt 	;

reg  			clc_de 		;
reg  	[7:0] 	clc_cnt 	;

wire  	[31:0] 	hist_sum 	;

always @(posedge clk)
	pre_vs_r 	<= 	pre_vs 	;

assign 	vs_pose 	<= 	~pre_vs_r && pre_vs ;

data_sync_ext   u1_data_sync_ext(
    .clka           (clk        ),
    .rst_n          (1'b1       ),
    .pulse_a        (vs_pose    ),
    .ext_pulse_a    (ext_vs_pose) 
);

/***********************行列计数*************************************/
always @(posedge clk)
    rst_n <=    ~ext_vs_pose ;

always @(posedge clk)
    if(rst_n == 0)
        hcnt    <=  0;
    else if(hcnt == COL - 1 && pre_de == 1'b1)
        hcnt    <=  0;
    else if(pre_de == 1'b1)
        hcnt    <=  hcnt+1;
    else 
        hcnt    <=  hcnt;

always @(posedge clk)
    if(rst_n == 0)
        vcnt    <=  0;
    else if(hcnt == COL - 1 && vcnt == ROW - 1 && pre_de == 1'b1)
        vcnt    <=  0;
    else if(hcnt == COL - 1)
        vcnt    <=  vcnt + 1;
    else 
        vcnt    <=  vcnt    ; 

//当数据到来时相应的灰度值数量+1  当数据结束时依次清除数据
always @(posedge clk)
	if(clc_de)
		hist_lut[clc_cnt] <= 32'd0;
	else if(pre_de==1)
		hist_lut[pre_data] <= hist_lut[pre_data] + 1;


//当一帧数据结束后，开始统计直方图使能
always @(posedge clk)begin
	if(!rst_n)begin
		inter_de 	<= 	0 ;
	end
	else if(inter_cnt == 8'd255)begin
		inter_de 	<= 	0;
	end
	else if(hcnt == COL -1 && vcnt == ROW - 1)begin
		inter_de 	<= 	1'b1;
	end
	else begin
		inter_de 	<= 	inter_de ;
	end
end

//依次统计直方图的灰度
always @(posedge clk)begin
	if(!rst_n)begin
		inter_cnt 	<= 0;
	end
	else if(inter_cnt == 8'd255)begin
		inter_cnt 	<= 0;
	end
	else if(inter_de == 1)begin
		inter_cnt 	<= inter_cnt + 1'b1;
	end
	else begin
		inter_cnt 	<= inter_cnt;
	end
end

//依次统计直方图灰度
always @(posedge clk)
	if(rst_n == 0)
		lut_add 	<= 	0;
	else if(inter_de)
		lut_add 	<= 	lut_add 	+ 	hist_lut[inter_cnt];
	else 
		lut_add 	<= 	'd0;

//ram写使能拉高
always @(posedge clk )
	if(rst_n == 0)
		wea 	<= 	1'b0;
	else 
		wea 	<= 	inter_de;

//ram写地址
always @(posedge clk)
	if(rst_n == 0)
		addra 	<= 'd0 ;
	else if(addra == 8'd255)
		addra 	<= 	'd0;
	else if(wea == 1'd1)
		addra 	<= 	addra + 1'b1;
	else 
		addra 	<= 	addra ;

//统计完成后清除灰度直方图数据
always @(posedge clk)
	if(rst_n == 0)
		clc_de 	<= 	1'b0;
	else 
		clc_de 	<= 	inter_de ;

always @(posedge clk)
	if(rst_n == 0)
		clc_cnt 	<= 0;
	else 
		clc_cnt 	<= 	inter_cnt ;

//存储灰度直方图的ram
hist_ram u1_hist_ram (
	.clka		(clk			),  	
	.wea		(wea			),   	
	.addra		(addra			), 	
	.dina		(lut_add		),  	
	.douta		(				), 	

	.clkb		(clk			),  	
	.web		(1'b0			),   	
	.addrb		(pre_data		), 	
	.dinb		(32'd0			),  	
	.doutb		(hist_sum		)  	
);

reg  [0:1]	pre_de_r 	;

reg  	[63:0] 	mult 	;

localparam 	[63:0]  	coe 	= 	(2**25)*255/(COL*ROW) ; //累计直方图归一化以及重映射需要乘以的数据   乘数是累计直方图的灰度值

always @(posedge clk)
	if(rst_n == 0)
		pre_de_r 	<= 	0;
	else 
		pre_de_r 	<= 	{pre_de_r[0],pre_de};

//图像的数据进行重映射
always @(posedge clk)
	if(rst_n == 0)
		mult 	<= 	0;
	else if(pre_de_r[0] == 1'b1)
		mult 	<= 	hist_sum * coe;
	else 
		mult 	<= 	'd0;

always @(posedge clk)
	if(rst_n == 0)
		post_de 	<= 	1'b0;
	else 
		post_de 	<= 	pre_de_r[1];


//图像数据的位截取
always @(posedge clk)
	if(rst_n == 0)
		post_data 	<= 	1'b0;
	else if(pre_de_r[1] == 1'b1)
		post_data 	<= 	mult[25+7:25];
	else 
		post_data 	<= 	'd0;