[数字图像处理](五)AHE

先把一个代码放这里，、、、、

%ahe which I implemented

clc;
close all;
clear all;

rgbim = imread('p1.jpg');

subplot(1,2,1);
imshow(rgbim);
ah = ahe(rgbim,4,256);
subplot(1,2,2);
imshow(ah);



    function img_new=ahe(img,grid,limit) 
    
    img_new=img; %img_new is the img after AHE
    [m,n]=size(img);  
    
    %ECR
    
    grid_cols=grid;
    grid_rows=grid;

    grid_width=int32(fix(m/grid_cols));
    grid_height=int32(fix(n/grid_rows));

    
    map=zeros(grid_cols,grid_rows,256);
    
    %for each grid,we create their mapping function
    for i=1:grid_cols
        for j=1:grid_rows
            map(i,j,:)=MakeHistogram(img,1+(i-1)*grid_width,1+(j-1)*grid_height,grid_width,grid_height,limit); 
        end
    end
    
    

%interpolate
%boundary cases I followed the Karel Zuiderveld's implement(C version)
xi = 1; 
for i = 1:grid_cols+1 
    if i == 1 
        subx = grid_width/2; 
        xu = 1; 
        xd = 1; 
    elseif i == grid_cols+1 
        subx = grid_width/2; 
        xu =  grid_cols; 
        xd =  grid_cols; 
    else 
        subx = grid_width; 
        xu = i - 1; 
        xd = i; 
    end 
    yi = 1; 
    for j = 1:grid_rows+1
        if j == 1 
            suby = grid_height/2; 
            yl = 1; 
            yr = 1; 
        elseif j == grid_rows+1
            suby = grid_height/2; 
            yl = grid_rows; 
            yr = grid_rows; 
        else 
            suby = grid_height; 
            yl = j - 1; 
            yr = j; 
        end 
        UL = map(xu,yl,:); 
        UR = map(xu,yr,:); 
        DL = map(xd,yl,:); 
        DR = map(xd,yr,:); 

        subimg = img(xi:xi+subx-1,yi:yi+suby-1); 
        subimg = Interpolate(subimg,UL,UR,DL,DR,subx,suby); 
        img_new(xi:xi+subx-1,yi:yi+suby-1) = subimg; 
        yi = yi + suby; 
    end 
    xi = xi + subx; 
end 

    end
    
    
   

function map=MakeHistogram(img,startx,starty,width,height,limit) 
%make histogram
hist=zeros(1,256); 
        for i=startx:startx+width-1 
            for j=starty:starty+height-1
                value=img(i,j);
                hist(value+1)=hist(value+1)+1;
            end 
        end 
%clip and allienate noise

        if (limit>0)
        excess = 0; 
        %get the excess
        for degree = 1:256 
            excess=excess+max(0,hist(degree)-limit); 
        end 
        ts=limit-excess/256;
        avginc=excess/256;
        %cut them and compensate other gray level
        for degree=1:256
            if (hist(degree)>limit)
                hist(degree)=limit;
            end
            %when compensating,should not larger than limit
            if (hist(degree)>ts)
                excess=excess-(limit-hist(degree));
                hist(degree)=limit;
            else
                hist(degree)=hist(degree)+avginc;
                excess=excess-avginc;
            end
            
        end
        %if there are still some excess,equal divide them to all of the
        %gray levels
        avginc=excess/256;
        for degree=1:256
            hist(degree)=hist(degree)+avginc;
        end
        end
%create mapping function


        

        map=zeros(1,256);
        map(1)=hist(1);
        for i=2:256
            map(i)=map(i-1)+hist(i);
        end
        for i=1:256  
            map(i) = map(i)*256/(width*height);
            if map(i) > 256  
                map(i) = 256;  
            end  
        end

         
end

function new=Interpolate(img,lutUL,lutUR,lutDL,lutDR,width,height) 
%interpolate points according to the following formula
%f(x,y)=f(0,0)(1-x)(1-y) +f(0,1)(1-x)y+ f(1,1)xy+ f(1,0)x(1-y)
new=zeros(width,height);
for i=1:width
    for j=1:height
        value=img(i,j);
        a=double(width);zz
        c=double(i-1);
        b=double(height);
        d=double(j-1);
        
        x=c/a;
        y=d/b;
                
        new(i,j)=fix(lutUL(value+1)*(1-x)*(1-y)+lutDL(value+1)*x*(1-y)...
            +lutUR(value+1)*(1-x)*y+lutDR(value+1)*x*y);
        
    end
end

end