写在前面

 POSHE算法原理解读及c++代码实现请看:https://blog.csdn.net/weixin_40647819/article/details/88416512

 这次用MATLAB实现。

代码

git 下载地址:https://github.com/2209520576/Image-Processing-Algorithm/tree/master/Grayscale%20transformation/POSHE_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');

结果