Don't do it pixel by pixel. That'll take forever. At most, do it channel by channel or frame by frame. There are a lot of ways to do something like this. I'm going to choose to do it in a fairly basic and direct way.
% this is just my test image; not included
inpict=imread('sources/table.jpg');
graypict=im2double(rgb2gray(inpict));
% we want how many gray levels
graylevels=5;
% what's our colormap?
cmap=[0 0 1; 0 1 0; 1 1 0; 1 0 0; 0 0 0];
% number of graylevels should correspond to the length of cmap
% get map of nans
nanmap=isnan(graypict);
% allocate output array
outpict=zeros([size(graypict,1) size(graypict,2) 3]);
% permute the colormap
cmap=permute(cmap,[1 3 2]);
for gl=1:graylevels
% make a logical mask for this range of gray values
mask=(graypict>=((gl-1)/graylevels)) & (graypict<=(gl/graylevels));
% step through the image channels and assign color values to the mask region
% this can be done with bsxfun, but let's not get confusing.
for c=1:3
thischannel=outpict(:,:,c);
thischannel(mask)=cmap(gl,1,c);
outpict(:,:,c)=thischannel;
end
end
% inherit all those wonderful NaNs
for c=1:3
thischannel=outpict(:,:,c);
thischannel(nanmap)=NaN;
outpict(:,:,c)=thischannel;
end
imshow(outpict)
It's not going to win any awards for efficiency or elegance, but it's readable. It takes about 9 seconds to process a 4500x5500 image on a ~10 year old computer.
The per-channel loops could be avoided with bsxfun and/or array replication. The MIMT (on the file exchange) has a multilevel mask generation tool (mlmask()) and a tool for replacing image regions with a color tuple based on a single-channel mask (replacepixels()). There are probably other ways.
