子块部分重叠直方图均衡算法(POSHE算法)MATLAB实现
写在前面
POSHE算法原理解读及c++代码实现请看:https://blog.csdn.net/weixin_40647819/article/details/88416512
这次用MATLAB实现。
代码
1.POSHE.m
%% This is a function of POSHE.
%@ date : 2019/4/27
%@ author : Xiaowu
%@ CSDN :https://blog.csdn.net/weixin_40647819/article/details/88416512
%@ github :https://github.com/2209520576/Image-Processing-Algorithm
%//////////////////////////////////////
%POSHE
%src - 输入图像
%dst - 输出图像
%s - 子块尺寸=原图尺寸/s ,取2的倍数
%k - 子块移动步长=子块尺寸/k, 取2的倍数
%/////////////////////////////////////
%%
function dst=POSHE(src, s ,k)
if ndims(src)==3
src=rgb2gray(src);
end
[nr,nc]=size(src);
%边界扩充
newnr = ceil(nr / s / k)*s*k;
newnc = ceil(nc / s / k)*s*k;
newsrc=padarray(src,[newnr - nr newnc - nc],'replicate','post');
dst=zeros(newnr,newnc,'uint16'); %因为8位无符号整型的范围是0~255,而在累加过程中会超过这个范围,所以采用16位无符号整型
%创建子块,确定子块移动步长
sub_block_r = newnr / s; %子块高度
sub_block_c = newnc / s; %子块宽度
step_r = sub_block_r / k; %垂直移动步长
step_c = sub_block_c / k; %水平移动步长
%对当前子块进行直方图均衡
HE_frequency=zeros(newnr,newnc,'uint16'); %均衡频率计数矩阵
newsrc_draw=newsrc;
figure;imshow(newsrc_draw); title('子块移动网格'); %画子块产生的网格
for sub_block_x= 1 : step_r : (newnr - sub_block_r + 1)
for sub_block_y= 1 : step_c : (newnc - sub_block_c + 1)
sub_block=newsrc(sub_block_x:sub_block_x+sub_block_r-1 ,sub_block_y:sub_block_y+sub_block_c-1);
sub_block_HE=histeq(sub_block); %子块直方图均衡
rectangle('Position',[sub_block_y,sub_block_x,sub_block_c-1,sub_block_r-1],'edgecolor','c'); %画子块产生的网格
%将直方图均衡后子块的像素值映射至输出图像
sub_block_HE_i=1;
for i=sub_block_x:(sub_block_x + sub_block_r-1)
sub_block_HE_j = 1;
for j=sub_block_y:(sub_block_y + sub_block_c-1)
dst(i,j)=dst(i,j) + uint16(sub_block_HE(sub_block_HE_i,sub_block_HE_j));
HE_frequency(i,j)=HE_frequency(i,j) + 1;
sub_block_HE_j=sub_block_HE_j + 1;
end
sub_block_HE_i=sub_block_HE_i + 1;
end
end
end
dst=dst./HE_frequency;
dst=uint8(dst); %数据类型转换
%BERF
step_levels = 2;
large_dir = 15;
dst=BERF(newsrc, dst, step_r, step_c, step_levels, large_dir);
dst=imcrop(dst,[1,1,nc-1,nr-1]);
end
1.BERF.m
%% This is a function of BERF.
%@ date : 2019/4/27
%@ author : Xiaowu
%@ CSDN :https://blog.csdn.net/weixin_40647819/article/details/88416512
%@ github :https://github.com/2209520576/Image-Processing-Algorithm
%///////////////////////////////////////////////////
%BERF滤波器
%src - 输入图像
%dst - 目标图像
%step_r - 子块垂直移动步长
%step_c - 子块水平移动步长
%step_levels - 灰度级缓慢变化的步长
%large_dir - 最大灰度差异
%//////////////////////////////////////////////////
%%
function dst=BERF(src, dst,step_r, step_c, step_levels, large_dir)
[dstrows,dstcols]=size(dst);
%处理子块的行边界(横向边界)
for i=1:step_r:dstrows
for j=1:dstcols
if (i>1 && i<dstrows ) %图像边界判断,先排除第一行和最后一行
%块效应检查
dfacter_newsrc=abs(src(i,j)-src(i+1,j)) + abs(src(i,j)-src(i-1,j)); %原图子块边界与相邻两侧像素的灰度差异
dfacter_dst=abs(dst(i,j)-dst(i+1,j)) + abs(dst(i,j)-dst(i-1,j)); %POSHEed图子块边界与相邻两侧像素的灰度差异
if (abs(dfacter_dst - dfacter_newsrc) > step_levels)
%存在块效应执行BERF
b = 0;
ave_bound=uint8( (uint16(dst(i - 1, j)) + uint16(dst(i + 1, j))) / 2 );
dst(i,j)=ave_bound;
%垂直于子块边界,向上递增方向( increasing rule)
if (dst(i-1,j)> dst(i+1,j))
if (i - 2 - b >= 1) %图像边界判断(下一位置像素)
pixel_present_add = dst(i - 1 - b, j);
pixel_next_add = dst(i - 2 - b, j);
pixel_present_add = ave_bound + step_levels;
while(i - 2 - b >= 1 && pixel_next_add - pixel_present_add >= step_levels && pixel_next_add - pixel_present_add <= large_dir)
pixel_next_add = pixel_present_add + step_levels;
dst(i - 1 - b, j)=pixel_present_add;%%%%%%
dst(i - 2 - b, j)=pixel_next_add;%%%%%%
b =b+1;
if(i - 2 - b >= 1) %图像边界判断(下一位置像素)
pixel_present_add = dst(i - 1 - b, j);
pixel_next_add = dst(i - 2 - b, j);
end
end
end
%垂直于子块边界,向下递减方向( decreasing rule)
b=0;
if (i + 2 + b <=dstrows)
pixel_present_dec = dst(i + 1 + b, j);
pixel_next_dec = dst(i + 2 + b, j);
pixel_present_dec = ave_bound - step_levels;
while (i + 2 + b <=dstrows && pixel_present_dec - pixel_next_dec >= step_levels && pixel_present_dec - pixel_next_dec <= large_dir)
pixel_next_dec = pixel_present_dec - step_levels;
dst(i + 1 + b, j)=pixel_present_dec;%%%%%%
dst(i + 2 + b, j)=pixel_next_dec;%%%%%%
b = b+1;
if (i + 2 + b <=dstrows) %图像边界判断
pixel_present_dec = dst(i + 1 + b, j);
pixel_next_dec = dst(i + 2 + b, j);
end
end
end
end
%垂直于子块边界,向下递增方向( increasing rule)
b=0;
if (dst(i - 1, j) < dst(i + 1, j))
if (i + 2 + b <= dstrows) %图像边界判断
pixel_present_add = dst(i + 1 + b, j);
pixel_next_add = dst(i + 2 + b, j);
pixel_present_add = ave_bound + step_levels;
while (i + 2 + b <= dstrows && pixel_next_add - pixel_present_add >= step_levels && pixel_next_add - pixel_present_add <= large_dir)
pixel_next_add = pixel_present_add + step_levels;
dst(i + 1 + b, j)=pixel_present_add;%%%%%%
dst(i + 2 + b, j)=pixel_next_add;%%%%%%
b =b + 1;
if (i + 2 + b <= dstrows) %图像边界判断
pixel_present_add = dst(i + 1 + b, j);
pixel_next_add = dst(i + 2 + b, j);
end
end
end
%垂直于子块边界,向上递减方向( decreasing rule)
b=0;
if (i - 2 - b >= 1)
pixel_present_dec = dst(i - 1 - b, j);
pixel_next_dec = dst(i - 2 - b, j);
pixel_present_dec = ave_bound - step_levels;
while (i - 2 - b >= 1 && pixel_present_dec - pixel_next_dec >= step_levels && pixel_present_dec - pixel_next_dec <= large_dir)
pixel_next_dec = pixel_present_dec - step_levels;
dst(i - 1 - b, j)=pixel_present_dec;%%%%%%
dst(i - 2 - b, j)=pixel_next_dec;%%%%%%
b =b + 1;
if (i - 2 - b >= 1) %图像边界判断
pixel_present_dec = dst(i - 1 - b, j);
pixel_next_dec = dst(i - 2 - b, j);
end
end
end
end
end
end
end
end
%处理子块的行边界(纵向向边界)
for j=1:step_c:dstcols
for i=1:dstrows
if (j>1 && j<dstcols ) %图像边界判断,先排除第一列和最后一列
%块效应检查
dfacter_newsrc=abs(src(i,j)-src(i,j+1)) + abs(src(i,j)-src(i,j-1)); %原图子块边界与相邻两侧像素的灰度差异
dfacter_dst=abs(dst(i,j)-dst(i,j+1)) + abs(dst(i,j)-dst(i,j-1)); %POSHEed图子块边界与相邻两侧像素的灰度差异
if (abs(dfacter_dst - dfacter_newsrc) > step_levels)
%存在块效应执行BERF
b = 0;
ave_bound=uint8( (uint16(dst(i, j-1)) + uint16(dst(i, j+1))) / 2 );
dst(i,j)=ave_bound;
%垂直于子块边界,向左递增方向( increasing rule)
if (dst(i,j-1)> dst(i,j+1))
if (j - 2 - b >= 1) %图像边界判断(下一位置像素)
pixel_present_add = dst(i,j - 1 - b);
pixel_next_add = dst(i,j - 2 - b);
pixel_present_add = ave_bound + step_levels;
while(j - 2 - b >= 1 && pixel_next_add - pixel_present_add >= step_levels && pixel_next_add - pixel_present_add <= large_dir)
pixel_next_add = pixel_present_add + step_levels;
dst(i,j - 1 - b)=pixel_present_add;%%%%%%
dst(i,j - 2 - b)=pixel_next_add;%%%%%%
b =b+1;
if(j - 2 - b >= 1) %图像边界判断(下一位置像素)
pixel_present_add = dst(i,j - 1 - b);
pixel_next_add = dst(i,j - 2 - b);
end
end
end
%垂直于子块边界,向右递减方向( decreasing rule)
b=0;
if (j + 2 + b <dstcols)
pixel_present_dec = dst(i,j + 1 + b);
pixel_next_dec = dst(i,j + 2 + b);
pixel_present_dec = ave_bound - step_levels;
while (i + 2 + b < dstcols && pixel_present_dec - pixel_next_dec >= step_levels && pixel_present_dec - pixel_next_dec <= large_dir)
pixel_next_dec = pixel_present_dec - step_levels;
dst(i,j + 1 + b)=pixel_present_dec;%%%%%%
dst(i,j + 2 + b)=pixel_next_dec;%%%%%%
b = b+1;
if (j + 2 + b <= dstcols) %图像边界判断
pixel_present_dec = dst(i, j + 1 + b);
pixel_next_dec = dst(i, j + 2 + b);
end
end
end
end
%垂直于子块边界,向右递增方向( increasing rule)
b=0;
if(dst(i , j-1) < dst(i , j+1))
if (j + 2 + b < dstcols) %图像边界判断
pixel_present_add = dst(i, j + 1 + b);
pixel_next_add = dst(i, j + 2 + b);
pixel_present_add = ave_bound + step_levels;
while (j + 2 + b < dstcols && pixel_next_add - pixel_present_add >= step_levels && pixel_next_add - pixel_present_add <= large_dir)
pixel_next_add = pixel_present_add + step_levels;
dst(i,j + 1 + b)=pixel_present_add;%%%%%%
dst(i,j + 2 + b)=pixel_next_add;%%%%%%
b =b + 1;
if (j + 2 + b < dstcols) %图像边界判断
pixel_present_add = dst(i, j + 1 + b);
pixel_next_add = dst(i, j + 2 + b);
end
end
end
%垂直于子块边界,向左递减方向( decreasing rule)
b=0;
if (j - 2 - b >= 1)
pixel_present_dec = dst(i, j - 1 - b);
pixel_next_dec = dst(i, j - 2 - b);
pixel_present_dec = ave_bound - step_levels;
while (j - 2 - b >= 1 && pixel_present_dec - pixel_next_dec >= step_levels && pixel_present_dec - pixel_next_dec <= large_dir)
pixel_next_dec = pixel_present_dec - step_levels;
dst(i,j - 1 - b)=pixel_present_dec;%%%%%%
dst(i,j - 2 - b)=pixel_next_dec;%%%%%%
b =b + 1;
if (j - 2 - b >= 1) %图像边界判断
pixel_present_dec = dst(i, j - 1 - b);
pixel_next_dec = dst(i, j - 2 - b);
end
end
end
end
end
end
end
end
3.POSHE_demo.m
%% This is an example of POSHE.
%@ date : 2019/4/27
%@ author : Xiaowu
%@ CSDN :https://blog.csdn.net/weixin_40647819/article/details/88416512
%@ github :https://github.com/2209520576/Image-Processing-Algorithm
%//////////////////////////////////////////
%POSHE_demo
%src - 输入图像
%dst - 输出图像
%s - 子块尺寸=原图尺寸/s ,一般取2的倍数
%k - 子块移动步长=子块尺寸/k, 一般取2的倍数
%///////////////////////////////////////////
%%
clc;clear all;close all
src=imread('vessel.bmp');
if ndims(src)==3
src=rgb2gray(src);
end
figure;imshow(src);title('src');
%%
s=2; %建议2~4
k=16; %建议12~16
dst=POSHE(src, s ,k);
figure;imshow(dst);title('dst');
结果