How to display multiple ROI (masks) in different colors on one image?

I have problem displaying multiple ROI(masks) on one image. I am able to display one mask on the image, and then imfuse the images at the end. But this means that the colors of each mask becomes increasingly transparent, and i would like to avoid that. Please support with ideas and upgrades to the current code.

clc;  % Clear command window.
clear;  % Delete all variables.
close all;  % Close all figure windows except those created by imtool.
imtool close all;  % Close all figure windows created by imtool.
workspace;  % Make sure the workspace panel is showing.
fontSize = 16;

Image = imread('Report_Picture.png');
[rows, columns, numberOfColorChannels] = size(Image);
imshow(Image, []);
axis on;
title('Report Image', 'FontSize', fontSize);
% set(gcf, 'Position', get(0,'Screensize')); % Maximize figure.
message = sprintf('1 freehand');
uiwait(msgbox(message));
roi = images.roi.AssistedFreehand;
draw(roi);
message = sprintf('2 freehand');
uiwait(msgbox(message));
roi2 = images.roi.AssistedFreehand;
draw(roi2);
message = sprintf('3 freehand');
uiwait(msgbox(message));
roi3 = images.roi.AssistedFreehand;
draw(roi3);
message = sprintf('4 freehand');
uiwait(msgbox(message));
roi4 = images.roi.AssistedFreehand;
draw(roi4);
mask1 = createMask(roi);
mask2 = createMask(roi2);
mask3 = createMask(roi3);
mask4 = createMask(roi4);

% % Display the freehand mask.
% imshow(mask3);
% axis on;
% title('Binary mask of the region', 'FontSize', fontSize);

% If it's grayscale, convert to color
if numberOfColorChannels < 3
  rgbImage = cat(3, Image, Image, Image);
else
  % It's really an RGB image already.
  rgbImage = Image;
end
% Extract the individual red, green, and blue color channels.
redChannel = rgbImage(:, :, 1);
greenChannel = rgbImage(:, :, 2);
blueChannel = rgbImage(:, :, 3);
redChannel2 = rgbImage(:, :, 1);
greenChannel2 = rgbImage(:, :, 2);
blueChannel2 = rgbImage(:, :, 3);
redChannel3 = rgbImage(:, :, 1);
greenChannel3 = rgbImage(:, :, 2);
blueChannel3 = rgbImage(:, :, 3);
redChannel4 = rgbImage(:, :, 1);
greenChannel4 = rgbImage(:, :, 2);
blueChannel4 = rgbImage(:, :, 3);
% Specify the color we want to make this area.
desiredColor1 = [235, 18, 22]; % Red
desiredColor2 = [0, 255, 0]; % Green
desiredColor3 = [251, 255, 0]; % Yellow
desiredColor4 = [251, 255, 0]; % Yellow
% Make the red channel that color
redChannel(mask1) = desiredColor1(1);
greenChannel(mask1) = desiredColor1(2);
blueChannel(mask1) = desiredColor1(3);
redChannel2(mask2) = desiredColor2(1);
greenChannel2(mask2) = desiredColor2(2);
blueChannel2(mask2) = desiredColor2(3);
redChannel3(mask3) = desiredColor3(1);
greenChannel3(mask3) = desiredColor3(2);
blueChannel3(mask3) = desiredColor3(3);
redChannel4(mask4) = desiredColor4(1);
greenChannel4(mask4) = desiredColor4(2);
blueChannel4(mask4) = desiredColor4(3);
% Recombine separate color channels into a single, true color RGB image.
rgbImage1 = cat(3, redChannel, greenChannel, blueChannel);
imwrite(rgbImage1,'C:\Users\phillip.wulff\Ergo Files\Report\Pictures\Picture1.png');
rgbImage2 = cat(3, redChannel2, greenChannel2, blueChannel2);
imwrite(rgbImage2,'C:\Users\phillip.wulff\Ergo Files\Report\Pictures\Picture2.png');
rgbImage3 = cat(3, redChannel3, greenChannel3, blueChannel3);
imwrite(rgbImage3,'C:\Users\phillip.wulff\Ergo Files\Report\Pictures\Picture3.png');
rgbImage4 = cat(3, redChannel4, greenChannel4, blueChannel4);
imwrite(rgbImage4,'C:\Users\phillip.wulff\Ergo Files\Report\Pictures\Picture4.png');
% Display the image.

A = imfuse(rgbImage1,rgbImage2,'blend','Scaling','none');
A = imfuse(A,rgbImage3,'blend','Scaling','none');
A = imfuse(A,rgbImage4,'blend','Scaling','none');
imshow(A);
title('Image with color inside the mask region', 'FontSize', fontSize);

Answer 1

[DGM]

I know this is probably too late to help, but consider this:

fontSize = 16;
Image = imread('Report_Picture.png');
[rows, columns, numberOfColorChannels] = size(Image);

imshow(Image, []);
axis on;
title('Report Image', 'FontSize', fontSize);
disp('Make freehand selection 1 of 4')
%message = sprintf('1 freehand');
%uiwait(msgbox(message)); % wow that's annoying
roi = images.roi.AssistedFreehand;
draw(roi);
disp('Make freehand selection 2 of 4')
roi2 = images.roi.AssistedFreehand;
draw(roi2);
disp('Make freehand selection 3 of 4')
roi3 = images.roi.AssistedFreehand;
draw(roi3);
disp('Make freehand selection 4 of 4')
roi4 = images.roi.AssistedFreehand;
draw(roi4);

mask1 = createMask(roi);
mask2 = createMask(roi2);
mask3 = createMask(roi3);
mask4 = createMask(roi4);

% If it's grayscale, convert to color
if numberOfColorChannels < 3
    rgbImage = cat(3, Image, Image, Image);
else
    % It's really an RGB image already.
    rgbImage = Image;
end

% Specify the color we want to make this area.
desiredColor1 = [235, 18, 22]; % Red
desiredColor2 = [0, 255, 0]; % Green
desiredColor3 = [251, 255, 0]; % Yellow
desiredColor4 = [251, 255, 0]; % Yellow

% generate overlays from masks
olay1 = 255-mask1.*permute(255-desiredColor1,[1 3 2]);
olay2 = 255-mask2.*permute(255-desiredColor2,[1 3 2]);
olay3 = 255-mask3.*permute(255-desiredColor3,[1 3 2]);
olay4 = 255-mask4.*permute(255-desiredColor4,[1 3 2]);
alloverlays = (olay1 + olay2 + olay3 + olay4)/4;

% combine overlay with image
A = uint8((double(rgbImage) + alloverlays)/2);

% Display the image.
imshow(A);
title('Image with color inside the mask region', 'FontSize', fontSize);

There are a lot of ways you can combine images, but using something like imfuse() can become cumbersome and limiting. The above is the source image averaged with the average of overlays. In this way all overlays are transparent, they have the same weight and the source is subdued.

If you don't want the background subdued, just do a darken-only blend:

A = min(uint8(alloverlays),rgbImage);

If you actually want the masked regions to not be transparent, then it'll have to be decided whether they should be transparent with respect to each other, otherwise overlapping regions will obscure each other (might not matter if they never overlap). Let's assume that they can blend with each other, but not the source image:

allmasks = repmat(mask1|mask2|mask3|mask4,[1 1 3]);
A = rgbImage;
A(allmasks) = uint8(alloverlays(allmasks));

If they should stack atop each other, then alloverlays would be calculated differently.

[Gwendolyn Williams]

I realize I'm pretty late to the party here, but I could use some help - you stated that if we didn't want the background subdued, we could " just do a darken-only blend". When I changed the line defining A to do this according to your code, it gives me the following error:

"Integers can only be combined with integers of the same class, or scalar doubles."

So, I tried changing my base image to a uint8 as well, and got out an image that only showed a very small portion of my ROIs (the first ROI didn't show up at all - see attached picture).

I also tried changing the alloverylays to uint16 before taking the min, but that resulted in a completely black image. I'm not sure why, so I realize my understanding of exactly what's happening in regards to the coloring and use of overlays is insufficient here but I can't seem to figure it out on my own.

There was one difference in my sitation than that outlined above. My ROIs were loaded in from NIFTI files as uint8, which I converted into logicals using the logical() function rather than using createMask() as you had done with the freehand ROIs. I assume this was a valid workaround because I am able to get an image out thats completely subdued but shows the three ROIs correctly with the base image (image below).

My base image was a black and white uint16, which I converted to RGB using the if statement given above, and I only had 3 ROIs so I changed the 4 divisor to a 3 in the alloverlays definition.

Other than that, everything I'm doing is exactly the same as outlined above, yet using the darken-olnly blend fails to show my ROIs transparently on top of my base image.

To be safe, my code:

pons = niftiread('V03_pons.nii.gz'); %read ROIs as uint18
cereb = niftiread('V03_cereb.nii.gz');
spinal = niftiread('V03_SC.nii.gz');

pons_rot = imrotate(pons,90); %sets nifti ROI to the right orientation
cereb_rot = imrotate(cereb,90);
spinal_rot = imrotate(spinal,90);

mag_image = dicomread('V03_mag-0001.dcm'); %read in dicom base image

[rows, columns, numberOfColorChannels] = size(mag_image);

if numberOfColorChannels < 3
    rgbImg = cat(3, mag_image, mag_image, mag_image);
else
    rgbImg = mag_image;
end

desiredColor1 = [235, 18, 22]; % Red
desiredColor2 = [0, 255, 0]; % Green
desiredColor3 = [251, 255, 0]; % Yellow

mask1 = logical(pons_rot);
mask2 = logical(cereb_rot);
mask3 = logical(spinal_rot);

olay1 = 255-mask1.*permute(255-desiredColor1,[1 3 2]);
olay2 = 255-mask2.*permute(255-desiredColor2,[1 3 2]);
olay3 = 255-mask3.*permute(255-desiredColor3,[1 3 2]);

alloverlays = (olay1 + olay2 + olay3)/3;

A = uint8((double(rgbImg) + alloverlays)/2);
% A = min(uint8(alloverlays),uint8(double(rgbImg)));

imshow(A)

If I need to make a new thread for this question, I'm happy to.

Thanks in advance!

[DGM]

A couple quick things:

A lot of the image blending/compositing code I post on the forums is easily broken by changes in input class. As you've found out, there are limitations to mixed-class arithmetic. To avoid wrapping every example in a bunch of distracting input validation/casting, I simply omit that part. I'm used to using MIMT tools that have simple and compact ways of dealing with this, but since that's a custom toolbox, I usually avoid making my answers dependent on it.

The darken-only blend works well with the lineart because the core features of the image objects are all dark lines. Your image is both visually inverted (comparitively speaking) and continuous. It has a continuum of gray levels with no discrete foreground or background. The inversion means that darken-only blending really is the wrong approach to use. A lighten-only blend may be more effective, but the fact that the image is relatively continuous may mean that such a blend would produce unsatisfactory results due to the creation of new edge features at the edge of the blended area.

If you give me a bit, I will come up with an expanded example.

[DGM]

For this example, I'm simply going to use MIMT tools. As mentioned, this simplifies everything dramatically, but that also means I can't run this in the forum editor.

I took the first of your images and simply desaturated it. In this example, inpict is MxNx1 uint8, and the overlays are (together) MxNx3x4 uint8. This will still work even if inpict is floating point or if it's MxNx3.

The masks are MxNx1 logical, but they could also be MxNx1 unit-scale (0-1) floating point.

MIMT imblend() can blend images of equal geometry (height and width) using any one of many blend modes and compositing modes, as one would blend layers in Photoshop. Its primary parameters are opacity (scalar 0-1), the blend mode name, and an optional blend parameter. The user can also optionally specify a compositing mode and parameter. The documentation lists and describes the modes and options.

For this section, I'm going to stick to applying the overlays primarily using image blending. Subsequent examples will go into more compositing-oriented approaches.

inpict = imread('xray_blend.jpg');
%inpict = im2double(inpict); % imblend is class-agnostic
%inpict = repmat(inpict,[1 1 3]); % and it can combine mixed-depth images

sz = imsize(inpict,2); % explicitly get size of first two dimensions

% i'm assuming all the masks have the same geometry as inpict
nmasks = 4;
maskstack = zeros([sz 1 nmasks]);
maskstack(:,:,:,1) = imread('xrmask1.png');
maskstack(:,:,:,2) = imread('xrmask2.png');
maskstack(:,:,:,3) = imread('xrmask3.png');
maskstack(:,:,:,4) = imread('xrmask4.png');

% Specify the color we want to make each area.
desiredColortable = [235, 18, 22;  % red
    0, 255, 0; % green
    251, 255, 0; % yellow
    200, 65, 255]; % purple

% generate overlays from masks 
% this is a 4D image containing all four RGB overlays
olaystack = uint8(repmat(maskstack,[1 1 3 1]).*permute(desiredColortable,[4 3 2 1])); % black BG
%olaystack = uint8(255-repmat(maskstack,[1 1 3 1]).*permute(255-desiredColortable,[4 3 2 1])); % white BG

outpict = inpict;
for f = nmasks:-1:1
    outpict = imblend(olaystack(:,:,:,f),outpict,1,'lightenrgb',1);
end

imshow(outpict)

This is the result of a lighten-only blend. This mode is primarily just a per-channel max() operation, so it's easy to see that certain regions rapidly get saturated at white. Even in unsaturated regions (the yellow, green) the nonlinear behavior of a max() operation means that much of the result is flat. Local detail has been lost.

One way to reduce the effect of a blending operation is simply to reduce the opacity of the final composition:

outpict = imblend(olaystack(:,:,:,f),outpict,0.5,'lightenrgb',1);

This makes BG detail more visible again, though the result tends to look a bit washed out and flat still, as the output is still a function of that heavily-saturated intermediate image.

This can be done with other modes. Consider 'screen'.

outpict = imblend(olaystack(:,:,:,f),outpict,0.5,'screen',1);

Most other modes have an optional blend parameter which can be used to adjust the behavior without resorting to tweaking the output composition via opacity. This is 'linear dodge' (addition) at full opacity and half-slope.

outpict = imblend(olaystack(:,:,:,f),outpict,1,'lineardodge',0.5);

Similarly with 'color dodge'.

outpict = imblend(olaystack(:,:,:,f),outpict,1,'colordodge',0.5);

Note that 'colordodge' cannot raise black areas, whereas 'lineardodge' and others can. Any modes with a neutral response at BG = 0 will behave this way.

Since these overlays have a black background, any blend mode with a neutral response locus (NRL) at FG = 0 should work. That includes all dodge modes, hybrid glow modes, relational lighten modes, screen, addition, subtraction, difference, negation, etc.

You could alternatively (though I doubt it's desired) try blending toward darker colors instead. If you switch out the appropriate line, you can create the overlays with a white background instead of black.

olaystack = uint8(255-repmat(maskstack,[1 1 3 1]).*permute(255-desiredColortable,[4 3 2 1])); % white BG

This would allow you to use modes with a NRL at FG = 1 (burn modes, hybrid shadows, darkenrgb, multiply, etc).

outpict = imblend(olaystack(:,:,:,f),outpict,0.5,'multiply',1);

outpict = imblend(olaystack(:,:,:,f),outpict,1,'linearburn',0.5);

As you might guess, these two groups of modes are complements of each other.

Once installed from the FEX link above, the primary documentation for imblend() can be found with

help imblend

The webdocs have some particular examples, though mostly of unrelated component modes:

http://mimtdocs.rf.gd/manual/html/imblend.html

The FEX package includes a PDF with tabular and graphical descriptions of most modes, including information about the neutral response characteristics.

[DGM]

The above describes the use of imblend() primarily as one might do the task in Photoshop or GIMP; i.e. the primary focus is in blending, leveraging the NRL behavior of the modes to control what parts of the image are unaffected. One can also approach this task using compositing or any combination of both.

In this example, the overlay stack contains four frames, each is an RGBA image where the color channels are solid. Only the alpha content describes the individual ROIs.

inpict = imread('xray_blend.jpg');

sz = imsize(inpict,2);

% i'm assuming all the masks have the same geometry as inpict
nmasks = 4;
maskstack = zeros([sz 1 nmasks]);
maskstack(:,:,:,1) = imread('xrmask1.png');
maskstack(:,:,:,2) = imread('xrmask2.png');
maskstack(:,:,:,3) = imread('xrmask3.png');
maskstack(:,:,:,4) = imread('xrmask4.png');

% Specify the color we want to make each area.
desiredColortable = [235, 18, 22;  % red
    0, 255, 0; % green
    251, 255, 0; % yellow
    200, 65, 255]; % purple

% generate overlays from masks 
% this is a 4D image containing all four RGBA overlays
olaystack = uint8(ones([sz 3 nmasks]).*permute(desiredColortable,[4 3 2 1])); % solid color fields
olaystack = cat(3,olaystack,imcast(maskstack,'uint8')); % with alpha

outpict = inpict;
for f = nmasks:-1:1
    outpict = imblend(olaystack(:,:,:,f),outpict,0.25,'normal',1,'gimp',1);
end
outpict = outpict(:,:,1:3); % strip alpha

imshow(outpict)

By default, imblend() replicates GIMP compositing (approximately SRC-OVER or SRC-ATOP depending on the blend mode). You could alternatively try something else

outpict = imblend(olaystack(:,:,:,f),outpict,0.25,'normal',1,'translucent',1);

As mentioned, this can be combined with blending. As the area affected by the operations is defined by the overlay alpha instead of its color, you could use any blend mode regardless of its NRL.

[DGM]

If alpha compositing is all that's needed, imblend() ends up being a bit overkill. Within MIMT, there are simpler tools. Included in the same package is replacepixels(). This takes a background image, a foreground image or color tuple, and a mask. The mask does not need to be binarized, and it can be RGB.

inpict = imread('xray_blend.jpg');

sz = imsize(inpict,2);

% i'm assuming all the masks have the same geometry as inpict
nmasks = 4;
maskstack = zeros([sz 1 nmasks]);
maskstack(:,:,:,1) = imread('xrmask1.png');
maskstack(:,:,:,2) = imread('xrmask2.png');
maskstack(:,:,:,3) = imread('xrmask3.png');
maskstack(:,:,:,4) = imread('xrmask4.png');

% Specify the color we want to make each area.
desiredColortable = [235, 18, 22;  % red
    0, 255, 0; % green
    251, 255, 0; % yellow
    200, 65, 255]; % purple

opacity = 0.4; % the opacity for all overlays

outpict = inpict;
for f = nmasks:-1:1
    outpict = replacepixels(uint8(desiredColortable(f,:)),outpict,opacity*maskstack(:,:,:,f));
end

imshow(outpict)

The code is marginally simplified, but the drastic simplification of parameter options may make this approach more appealing.

Note that in this example, the input classes don't match. Outpict is uint8, whereas the individual masks are double (they start out as logical, but multiplying by opacity implicitly casts them as double). The only reason that the color tuple being fed to replacepixels() is cast uint8 isn't to make it compatible with the class of outpict, but rather because the values in the color tuple are appropriate for uint8 (0-255). As with many tools, the data range should correspond to the class. Similarly, the overlays in the prior example were cast uint8 not to match the background image, but rather to match its own data range.

[Gwendolyn Williams]

Thank you for all of this! Truly, the multiple examples and detailed explanations are extremely helpful. I tend to get confused with understanding what items need to be of what types and why, and the nuanced differences between functionality of different data types. I was able to get each of your examples to work on my end, and MIMT is going to be very helpful for me going forwards with this specific project. Again, thank you so much for this detailed and informative response!!

[DGM]

FWIW, that package also includes imcast(), which is basically a replacement for im2uint8(), im2double(), etc. Its advantage is that the class name is a parameter. This allows you to simplify workflows to avoid having to deal with potential dependency on the class of the input image. Trying to do the same with the IPT tools requires a bunch of conditionals and is kind of a hassle.

% say you read in some image that might be of any class
inpict = imread('cameraman.tif');

% regardless of class of inpict, cast & scale to double
% but keep track of its original class
[workingpict inclass] = imcast(inpict,'double');

% do a bunch of hypothetical operations with other floating-point arrays
% that might break if workingpict were integer-class
npict = rand(size(workingpict));
workingpict = workingpict.*npict;

% at the end of the script or function
% cast & scale the result so that it matches the original image
outpict = imcast(workingpict,inclass);

Of course, always working in floating point is expensive in terms of memory, but it is at least consistent and generally more convenient to write.