histogram of RGB combinations in matlab

Question

I have an RGB image (of type uint8) and I want to perform the following process:

1. Quantize the image to 16 levels (for each layer)
2. calculate the histogram of each RGB combination: first I need a 16*16*16 matrix to hold those values - i.e if we denote the histogram matrix hist then hist(2,9,3) should hold the amount of pixels in the histogram with intensity levels of r=2, g=9, b=3 etc. the next stage is to reshape it into a 4096*1 vector (but that's the easy part)

for now I have the following implementation:

function hist = compRGBHist(I)
    I = I./16 % quantize to 16 color levels
    hist = zeros(16,16,16);
    [m, n, ~] = size(I);
    for i=1:m
        for j = 1:n
            rgb = impixel(I,j,i);
            r = rgb(1);
            g = rgb(2);
            b = rgb(3);
            A(r,g,b) = A(r,g,b) + 1;
        end
    end
    hist = reshape(hist, 4096, 1);
end

this implementation works, but it is VERY slow - I need to repeat the process a 100 times (as part of implementing a particle filter), and even though performing it on quite small images (actually image portions) of size ~80*40 it takes very long time. I'm looking for a more efficient way to do this. thanks!

Answer 1

You can turn the 4-bit values of RGB into a 12-bit value using bit shifts:

I = uint16(bitshift(I,-4));   % make sure the integer is large enough
twelvebit = bitshift(I(:,:,1),8) + bitshift(I(:,:,2)),4) + I(:,:,3);

This will result in values that look like this:

0 0 0 0   R R R R   G G G G   B B B B
-------   -------   -------   -------
 unused     red      green      blue

Now that your values are in the range [0..4095] you can calculate the histogram over 4096 bins.

Answer 2

If you divide a uint8 (0..255) by 2^4 you end up with a range of (0..16) instead of (0..15). You should first subtract 2^3 I think, to make sure you end up with exactly 16 bins.

I = randi([0 255],512,512,3,'uint8'); %example data
I = (I-2^3)./2^4;
I = uint16(I);
I(:,:,2)=I(:,:,2)*2^4;
I(:,:,3)=I(:,:,3)*2^8;
I = sum(I,3,'native');
h=hist(I(:),[0:4095]);
h=reshape(h,[16 16 16]);

Answer 3

I like accumarray here :

function hist = compRGBHist(I)
    I = bitshift(I, -4); % quantize to 16 color levels
    R = I(:,:,1) + 1; % +1 to be in [1..16]
    G = I(:,:,2) + 1;
    B = I(:,:,3) + 1;
    A = accumarray([R(:), G(:), B(:)], 1, [16,16,16]);
    hist = reshape(A, 4096, 1);
end

PS : Beaker's, bitshift solution is probably the best.