Drawing Multiple surf Plots with Different Color Map and Different Transparencies

Hello MATLAB Community! I've been working on CT slice visualising. I've extracted CT slices from DICOM file and plotted RT contours on top of it using the hold on. I used surf function to generete slice view and plot function to plot contours on top. Right now the view is something like this,

This is the main viewing angle, from the side plot looking like this,

I used surf function to plot the slice because the geometry manupulation of the CT geometry is easier and its working well so far. What I like to do now is to overlay dose map on top of this. Without CT, surf dose looking like this,

From the side,

The dose array has the same structure as the CT slices. I want to overlay this on top of the CT slice and when I looking from the 90 angle (like in the main viewing angle from the first picture) I want to see CT slice, the dose heatmap and contours as well. I tried changing alpha values with hold on but it did not work. How can I achieve such thing?

[Rahul]

can you share the code to generate the above plots?

[DGM]

Do you actually need to view the plots off-axis, or is a 2D representation sufficient?

[Arda Kayaalp]

2D Representation is sufficient but I found easier to manipulate with coordinates with the surf command, but I'm open to any suggestions.

[Arda Kayaalp]

This is the snippet, I added some comments,

% Import the CT Array as Array
% CT Array has dimension of 200x200x303
% Meaning that on Axial View each image
% has 200x200 pixels and 303 slices in total
Array = dicomData.CT.Array;

% Lets go through each slice in the array
for slice = 1:303
    structureZCoordinate = zCT(slice); % Found the desired slice's z-coordinate
    [X, Y] = meshgrid(xCT, yCT); % Create a mesh array from x and y coordinates
    Z = squeeze(Array(:,:, slice)); % Extract the desired slice from the array
    surf(X, Y, Z, 'LineStyle', 'none'); % Generate a surf plot from them
    cmax = ceil(max(Array(:))); % Find Max
    caxis([0,cmax]); % Set limits
    colormap(gray(1024)); % Use colormap gray but in the future Hounsfield Unit will be used
    %
    %% Options
    if showScale
        colorbar('location','eastoutside');
    end
    %                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                               
    if showLabels
        ylabel( {'\bf Y: Posterior \Leftarrow Anterior (mm)'});
        xlabel({'\bf X: Right \Rightarrow Left (mm)'});
    end
    %
    axis( [min(xCT) max(xCT) min(yCT) max(yCT)]) % Set axis limits
    view( 0,-90) % Set viewing angle
    %
    if showTitle
        title( {['\bf Z-Axis : ' num2str(structureZCoordinate),' (mm)'], ['\bf Slice: ' num2str(slice)]});
    end
    hold on %% hold on to print contours on top
    %% Print Contours
    
    if ~isempty(selectedStructures) % Check whether slice contour is empty or not
        for i = 1:length(selectedStructures) % Go throug all contours defined in the current slice
            structure = dicomData.Structures.(selectedStructures{i});
            B = bwboundaries(bwperim(squeeze(structure.indicatorMask(:, :,slice))),'noholes');% Create a boundary
            for j=1:size(B,1) % Go through boundary coordinates
                poly = B{j};
                structureXCoordinates = (poly(:,2)-1)*(xCT(2)-xCT(1))+xCT(1); % Extract the boundary coordinates 
                structureYCoordinates = (poly(:,1)-1)*(yCT(2)-yCT(1))+yCT(1);
                % Plot the boundary coordinates
                plot(structureXCoordinates,structureYCoordinates, 'LineWidth', 2, 'Color', contourColour(i,:));
            end
        end
    end
    % hold off
    % drawnow
    
    % => Until this point code is working as expected now I want to draw
    % dose array on top of all this, with CT and structure contours visible
    
    %% Print Dose Overlay
    doseArray = dicomData.Dose.Array; % Extract the dose array
    % Extract the dose coordinates
    xDose = dicomData.Dose.xDose; 
    yDose = dicomData.Dose.yDose;
    zDose = dicomData.Dose.zDose;
    % Extracted the dose map of the current slice
    structureZCoordinate = zDose(slice);
    [X, Y] = meshgrid(xDose, yDose);
    Z = doseArray(:,:,slice);
    s = surf(X, Y, Z,'LineStyle','none');
    % Same as before
    cmax = ceil(max(doseArray(:)));
    caxis([0,cmax]);
    % I want to use jet colormap, with highest dose represented with red
    % and lowest/none dose areas showed with blue but this also messes up
    % the view
    %   colormap(jet(1024));
    % I tried the blend in the view but not working
    alpha(s,0.1)
    hold off
    drawnow
    
end

here is the data set that I've been using,

https://drive.google.com/file/d/1VOhEfUr38IknfqaH-o6pzQ28K8w6yg69/view?usp=sharing

It's around 77 Mb, hence I uploaded to the drive.

Answer 1

[DGM]

Well, my network connection is super terrible, so I'm not going to be trying to download a bunch of stuff.

As given, the problem requires two (independently) colormapped objects in the same axes. That's not normally an option. There are possible ways around the issue:

• Merging all the colormaps into one long colormap and then trying to restrict the rendering of each object to a certain part of the composite colormap
• Converting one of the objects to an RGB image and then displaying it with image()/imagesc()/imshow(). Such an object is no longer colormapped, so that frees things up.
• Overlaying two axes objects atop each other, each with its own colormap.

Since you prefer the use something like surf() for handling the background image, I'm going to assume using image() or the like isn't an option. Tools like surf() and pcolor() can utilize full explicit XY data, whereas image() will always assume that XY coordinates are evenly spaced. If your XY data is evenly spaced, option 2 might be possible, though it would require changes to your existing work.

Following the third option, this might be one example.

% get an image with associated coordinates
imz = flipud(imread('cameraman.tif'));
[imx imy] = meshgrid(linspace(-200,200,256));

% get some placeholder dose map data
[X Y] = meshgrid(linspace(-100,100,100));
Z = exp(-(hypot(X,Y)/80).^2); % some Z data (a gaussian)

% get some placeholder boundary lines
boundaryxy = [3 -74;23 -22;86 -21;30 23;85 65;4 58;-9 101;
    -35 44;11 37;-5 13;-84 21;-92 -32;-11 -17;-81 -59;3 -74];

% create two axes
hax1 = axes;
hax2 = axes;

% pcolor and surf have similar syntax and behavior
% draw the background image in the lower axes
hp1 = pcolor(hax1,imx,imy,imz);
hp1.LineStyle = 'none';

% and draw the overlay in the upper axes, set its alpha accordingly
hp2 = pcolor(hax2,X,Y,Z);
hp2.LineStyle = 'none';
hp2.FaceAlpha = 0.5;

% draw the boundary curves in the upper axes, atop the dose map
hold on;
plot(boundaryxy(:,1),boundaryxy(:,2),'g')

% set colormaps
colormap(hax1,gray(256));
colormap(hax2,jet(256));

% wrangle axes setup
linkaxes([hax1 hax2]);
axis(hax1,'equal','tight'); 
set(hax2,'visible','off','xtick',[],'ytick',[]) % hide top axes decorations
set(hax2,'position',hax1.Position,'xlim',hax1.XLim,'ylim',hax1.YLim, ...
    'plotboxaspectratio',hax1.PlotBoxAspectRatio); % try to align both axes

Of course, trying to subsequently edit the axes positions or add things like legends or colorbars may provoke further alignment issues that will need to be resolved. In the thread linked above, there is some discussion of adding colorbars.

If it's desired to truncate the periphery of the dose map overlay, it should be possible to do that via the AlphaData property of hp2. You'd just need to generate a 2D map that describes which parts should be rendered.

[Arda Kayaalp]

Thank you for that detailed answer. Because my data isn't equally spaced, I found surf is the best way to go. The third option and the example you've provided helped a lot, with little change to the current code. Thanks a lot!