Given the way you're trying to do 1D lookup for the final remapping transformation, your input must be integer-valued and nonzero. If your inputs are unit-scale float, direct lookup by array indexing won't work. You'd need to actually do the interpolation.
That said, so long as your input is strictly uint8, your example works. It can be cleaned up significantly:
% setup
inpict = imread('cameraman.tif');
% transform the image
counts = histcounts(inpict,'binmethod','integers','binlimits',[0 255]);
myCDF = cumsum(counts)/numel(inpict);
outpict = 255*myCDF(double(inpict)+1);
outpict = uint8(outpict);
% test the results
opref = histeq(inpict,256); % the reference
immse(outpict,opref)
imshow([inpict opref outpict])
In this example, we're assuming:
- input is single-channel
- input is uint8 class
- number of input and output bins is 256
- target output distribution is uniform
Note that the example and reference outputs are visually similar, but not identical. They differ by only a few LSB, why? The fact is that because histeq() is designed to be more flexible than this integer-only code, it does the transformation calculations and output interpolation differently.
