How can I apply the jet-colormap to an image matrix with non-integer entries in an imwrite command?

Hello everybody :-),

I am currently working on my master's thesis and in desperate need of some help.

I have 800 grey scale images of a bone which I imread (works fine) .

for z=1:800
        filename = [num2str(z) '.tif'];
        s=([folder_path filename]);
        slice = imread(s);
        dim_x=size(slice,1);
        dim_y=size(slice,2);
        data2D=zeros(dim_x,dim_y,'uint8');
        data2D(:,:)=uint8(slice);

The grey values < GV_thr are 0 and the grey values GV_thr:256 define the osseous matter. Then, I want to replace values GV_thr:256 of the image matrix (='data2D') with values of the vector 'nu' which depend on the grey values, nu(GV); min(nu)=0.4544 and max(nu)=0.4657. (10200 is a factor to get rid of the decimal places and gain a certain range in respect of the integer-based colormap.) (works fine too, but very slowly)

        nu_rounded=round(nu*10200) % factor to get certain range for the colormap

        for i=GV_thr:256
            data2D(data2D(:,:) == i) = nu_rounded(i)
        end 

        maxVal = nu12_rounded(GV_thr); % (=4750)
        minVal = nu12_rounded(256); % (=4635)
        range  = maxVal-minVal; % (=115)

Now, I want to save the image with the same dimensions and everything, but in color. That can be achieved with 'imwrite' defining a colormap. However, if I apply 'cmap = [zeros(minVal,3); jet(range)]' I get an error, because maxVal > 256. So, I employed a trick and re-replaced the entries of the 'nu'-Vector in the image matrix (='data2D') with new grey values. (see below) This new grey value distribution is based on the intensity of 'nu' this time, which is the whole point of this.

        Vector_nu12=[minVal:maxVal]
        Vector_GV=[GV_thr:256]

        for ii=1:range
            data2D(data2D(:,:) == Vector_nu12(ii)) = Vector_GV(ii)
        end

        minVal1 = min(Vector_GV); % (=141)
        maxVal1 = max(Vector_GV); % (=256)
        range1  = maxVal1-minVal1;   % (=115)

The zeros of the image shall remain black and the bone tissue shall be demonstrated by the jet colormap.

        cmap = [zeros(minVal1,3); jet(range1)];

       % Imwrite new images
        imwrite(data2D, colormap(cmap), [num2str(fpath),num2str(z),'.tiff'])
        disp(['saved ',num2str(z),'/',num2str(num_images)])
end

However, the surrounding area of the bone is black as I wanted it to be, but the bone itself is just red and not in the jet color range. I don't understand what can be the reason for this as I already used this method applying only 1 replacement- step and it worked just fine.

I hope someone can spare some time to look through the code and can give me some hints. :-)

Thank you very much in advance!

[Image Analyst]

These are 16 bit images I presume? Can you attach one? Do you want each scaled independently according to its own max value or do you want a global max for all the images? And you want black up until the threshold, blue at the threshold and then jet up until the max value, either the max value of that particular image or of the global max, which ever you're using?

[Mohammad Abouali]

If all the bone materials are 256, then there would be only one color assigned to it. As image analyst said, perhaps it would be the best to upload a sample data2D in .mat format so we can see how that can be done.

[Gwen]

Thank you very much for your answers!!!

Actually, the images are already 8 bit grey scales, so converting them again as I did above is useless...

The threshold grey value of all 800 images is 140 and the values below are 0 and black (see image below). Yes, I want the new color images to be scaled according to the replaced values. The background has to stay black (=0), the minumum or new threshold value(=0.4544) I want to be blue and the maximum value (=0.4657) to be red; the range between the last 2 shall be defined by the jet colormap. So, I want a global max for all images and a global min.

The multiplication with 10200, the rounding and the second replacement is only an attempt to achieve this.

[Image Analyst]

If they are 8 bit, then how are you getting a max of 4750 and a min of 4635?

[Gwen]

I only get the max of 4750 and the min of 4635 through the substitution of the grey values in the image matrix (=data2D) by the values of the vector 'nu' which has the value 0.4657 at the GV_thr=140 and 0.4544 at GV=256:

nu(140)=0.4657

nu(256)=0.4544.

Then I multiplicated (factor 10200) and rounded the values of the nu vector and gained 4750 as max and 4635 as min. So I gained a new image matrix (=data2D) with new intensities and with enties ranging from 4635 to 4750 (bone) and with 0 (background). This was only an attempt to get integers for the jet colormap.

If there is a way to directly apply the jet colormap to the new threshold or minumum value (=0.4544) and maximum value (=0.4657) it would be even better.

[Guillaume]

Not answering your question at all, but avoid using the jet color map. See for example this page or this one

[Gwen]

Thank you! :-) I already read some of these pages and tried the range of one color, but my prof wants to have the jet color range and it would be consistent with my previous figures.

[Gwen]

However, I now know that the first replacement didn't work:

Instead of replacing the values 140 with 4750... and finally 255 with 4635, Matlab took the highest integer 255 on every non-zero place. So the jet colormap is successfully applied in the end and the bone is red, because there is only 1 integer left (255).

Summarizing that, my main problems here are that I can't replace the 8 bit grey values of the image matrix with values above 255 and that the jet colormap can't be applied to non-integers. Is there a solution to either of the two problems or another way around?

Answer 1

[Image Analyst]

Try this:

minValue = 0.4544
maxValue = 0.4657
% Make image with values from minValue to maxValue.
grayImage = (maxValue - minValue) * rand(480, 640) + minValue;
% Make some 0's in it.
zeroLocations = randperm(numel(grayImage), 10000);
grayImage(zeroLocations) = 0;
% Now we have our starting image.  
% It has some 0's with most values between .4544 and .4657.

% Map 0's to minValue so mat2gray will normalize over the range we want.
zeroLocations = grayImage == 0;
grayImage(zeroLocations) = minValue;
% Normalize the image to [0,255].
normalizedImage = uint8(255 * mat2gray(grayImage)); % Map to [0,255]
% Put zeros back in.
normalizedImage(zeroLocations) = 0;
% Convert to color image
rgbImage = ind2rgb(normalizedImage, jet(256));
imshow(rgbImage);

The first half is just to get an image with values in the "nu" range plus some scattered zeros. The second half does what you want I think.

[Gwen]

Thank you so much for your efforts!!! I will take a closer look at it in the evening.

The crucial point is however that I have to imread the image of the bone, get the image matrix with the GVs and then I want to replace the i. grey value with the i. element of the vector nu. This is important to settle the places and intensities of the entries in the new image matrix.

Answer 2

[Image Analyst]

Why not just use intlut() and ind2rgb()?

[Gwen]

The problem is that I don't have integer values in my image matrix any more but instead the values of nu = 0.4544 - 0.4657 (and 0). The reason why I replaced the GV with the nu- values at all is because the nu- values have a different distribution than the grey values: they are not linearly spaced like the grey values.

Answer 3

[Guillaume]

Instead of replacing the values 140 with 4750... and finally 255 with 4635, Matlab took the highest integer 255 on every non-zero place

That would be because you've defined data2D as uint8. The maximum value of a uint8 is 255. If you convert data2D to uint16 (or double) beforehand, you won't have that problem

works fine too, but very slowly

Use matrix indexing instead, avoids the loop. With your first example:

nu_rounded=round(nu*10200)
rescaled = nu_rounded(data2D); %as simple as that!

%to restore the values below gv_thr to their original values:
belowthreshold = data2D < GV_thr
rescaled(belowthreshold) = data2D(belowthreshold);

%or if you want to set all the values below the threshold to 0:
nu_rounded=round(nu*10200)
nu_rounded(1:GV_thr-1) = 0;
rescaled = nu_rounded(data2D)

[Gwen]

thank you, it would be really great to avoid the loop there. The image I imread is already 8 bit (my whole work is based on 8 bit images) but when I write 'uint16' instead of 'uint8' in the code and run it, I also get the error message: 'Subscript indices must either be real positive integers or logicals.'

[Guillaume]

Reading the values as uint8 is fine, but since you do this in a loop:

       data2D(data2D(:,:) == i) = nu_rounded(i)

The values of nu_rounded are coerced to 8 bit. So prior to your loop you would need to convert data2D to uint16 or double.

If you use matrix indexing as in

data2D = nu_rounded(data2D);

it's not a problem anymore as data2D will be whatever type nu_rounded is. The problem as you found out is that values in data2D start at 0, which is not a valid subscript index. You would have run into the same problem with your original code if your GV_thr was set to 0 by the way.

One possible workaround:

nu_rounded = [0 round(nu*10200)]; %shift values by one index
nu_rounded(1:GV_thr) = 0; %to set everything below GV_thr to 0
%or
nu_rounded(1:GV_thr) = 0:GV_thr to keep everything below GV_thr the same

data2D = nu_rounded(uint16(data2D)+1); %widen data2D so 1 can be added without overflow. Add 1 to avoid index 0 (hence the shift on the first line)

[Gwen]

Thanks! I will try that. :-)