function [b, pz, n] = histspec(imr) % histogram specification % r --> sr ---> sz ---> z % imr = imread('kids.tif', 'tiff'); % imz = imread('cameraman.tif', 'tiff'); % imr = imread('barbara.jpg', 'jpeg'); % imz = imread('building.jpg', 'jpeg'); [countsr, n] = imhist(imr); n = 0:255; countsz = 0.54 - 0.46*cos(n*2*pi/255); totsampr = prod(size(imr)); pr = countsr/totsampr; pz = countsz/sum(countsz); sr = cumsum(pr); % modify pz (and indirectly sz) to obtain distinc absisca for interp1 function k = pz == 0; pz(k) = 0.001; sz = cumsum(pz); sz = sz / sz(end); m = n/255; t = interp1(sz, m, sr, 'linear'); a = double(imr); [row, col] = size(a); b = zeros(row,col); c = ones(row, col); for n = 0:255 k = (a == n); b(k) = round(255*t(n+1)*c(k)); end b = uint8(b); n = 0:255; % figure; % subplot(211); plot(n, pz); % subplot(212); imhist(b);