% CLAHE
%
%  same as AHE, except clips (by cliplevel) 

function [ out ] = CLAHE( I, n, cliplevel, windowsize )

r = floor((windowsize-1)/2);
minI = min(min(I));
maxI = max(max(I));
step = (maxI - minI) / (n - 1);



H = zeros(n,1);

[height width] = size(I);
out = zeros(height, width);
area = 0;

for j=1:height
    lowj = max(1, j - r);
    highj = min(height, j + r);
    
    
    for i=(-r+1):(width+r+1)
        subi = i - r - 1;
        addi = i + r;
        
       
        if ( subi >= 1 )
           
            for jj = lowj:highj
                idx = floor(I(jj, subi) / step) + 1;
                H(idx) = H(idx) - 1;
            end
            
            area = area - (highj - lowj + 1);
        end
        
        
        if ( addi <= width )
            for jj = lowj:highj
                idx = floor(I(jj, addi) / step) + 1;
                H(idx) = H(idx) + 1;
            end
            area = area + (highj - lowj + 1);
        end
        
        if ( i >= 1 && i <= width )
           
            idx = floor(I(j, i) / step) + 1;
            val = 0;
            
           
            cropped = 0;
            SH = H;
			
			%clipping routine
			cliplevel = area * cliplevel;
            
            for k=1:n
                if ( SH(k) > cliplevel )
                    cropped = cropped + (SH(k) - cliplevel);
                    SH(k) = cliplevel;
                end
                
                
                
            end

            
            spread = cropped / n;
            for k=1:(idx-1)
                val = val + SH(k) + spread;
            end
         
            val = (val / area) * (n - 1);
            
            out(j, i) = val;
        end
    end
    
end


out = out / n;

end

