Matlab: "X-Ray" plot line through patch

Question

Problem

I'm attempting to visualize a 3-D path, along with a "cloud" around it which represents the standard deviation of data. I would like to be able to see a thick black line as the path with an EVENLY gray area around it, without any cloudiness to the lines, as if seeing the line through the cloud like an x-ray.

Attempt

I used plot3 to create a thick line and patch to create a series of boxes centered around each point of the line (there are some extra things on the plot to represent start/stop and direction, I would also like them to be seen easily). I tried playing with the alpha of the patch, but that creates a cloudiness to the line, such that the brightness of the gray color changes with however many gray boxes are in the line of sight. I would like the alpha to be 1, so that each gray box is exactly the same color, but I was hoping to find some way to make line seen through the cloud evenly.

Minimal Example

As requested, here is a minimal example, which produces the plot below.

% Create a path as an example (a circle in the x-y plane, with sinusoidal deviations in the z-axis)
t = 0:1/100:2*pi;
x = sin(t);y = cos(t);
z = cos(t).*sin(5*t);
figure;
plot3(x,y,z,'k','linewidth',7);

% Draw patches
cloud = .1*rand(size(t)); % The size of each box (make them random, "like" real data)
grayIntensity = .9; % Color of patch
faceAlpha = .15; % Alpha of patch

for i = 1:length(x)
patch([x(i) - cloud(i); x(i) + cloud(i); x(i) - cloud(i); x(i) + cloud(i); x(i) - cloud(i); x(i) + cloud(i); x(i) - cloud(i); x(i) + cloud(i)],... % X values
[y(i) - cloud(i); y(i) - cloud(i); y(i) + cloud(i); y(i) + cloud(i); y(i) - cloud(i); y(i) - cloud(i); y(i) + cloud(i); y(i) + cloud(i)],... % Y values
[z(i) + cloud(i); z(i) + cloud(i); z(i) + cloud(i); z(i) + cloud(i); z(i) - cloud(i); z(i) - cloud(i); z(i) - cloud(i); z(i) - cloud(i)],... % Z values
grayIntensity*ones(1,3),... % Color of patch
'faces', [1 2 4 3;5 6 8 7;1 2 6 5; 8 7 3 4;1 5 7 3;2 6 8 4],... % Connect vertices to form faces (a box)
'edgealpha',0,... % Make edges invisible (to get continuous cloud effect)
'facealpha',faceAlpha); % Set alpha of faces
end

Apologies for the VERY long stretch of code in the for loop, there are quite a large number of arguments to the patch command. The first three lines are simply defining the x, y, and z coordinates of the 8 vertices which define a cube, by specifying the center point plus or minus some half-width of the cube, cloud(i). The rest should be explained by their respective comments.

Thank you for any help!

Answer 1

Here is an implementation of the solution I was mentioning in the comment (just drawing one single surface cloud).

It is not optimized, there are a few for loop which could be avoided with clever use of bsxfun or these family of helper function, but it runs ok as it is. The math for finding the tangent of the curve at each point and orienting (rotating) each cross section could probably be simplified as well, but this is not my strong suit so I leave that to the experts if they feel like it.

Basically, it define a circle (often referred as "cross section" in the code) with a radius proportional to something (the standard deviation in your application, a random value in the example). Each circle is then rotated in 3D so it is normal to the curve at the point where it is translated. All these circles are then used as an envelope for a single surface graphic object.

There is still a bit of shadowing of the main centre line when the surface overlap itself many time (depending on the view angle), but the main line is always visible. Plus you only have one graphic object to manage.

The result looks like so:

Of course you can vary the AlphaValue of the surface to your liking. I defined a full matrix the same size as the data for the colour information. At the moment it is all set to 0 (so it points to the green colour in the default colormap), but this way it is also easy if you want to make the colour function of another parameter, just adjust the colour matrix accordingly (and the colormap that goes with it).

There is an option at the end of the code to display each cross section as a patch object. It is not intended to be used in the final result, but it can help you understand how the whole thing is constructed if you want to make your own modifications.

Here goes the code:

%% // Create a path as an example (a circle in the x-y plane, with sinusoidal deviations in the z-axis)
nPts = 180 ;
t = linspace(0,359,nPts)*pi/180;
x = sin(t); y = cos(t);
z = cos(t).*sin(2*t);

figure;
h.line = plot3(x,y,z,'k','linewidth',2,'Marker','none');
hold on
xlabel('X')
ylabel('Y')
zlabel('Z')

%% // Define options
%// cloud = .1*rand(size(t)) ; % The size of each box (make them random, "like" real data)
%// I used another randomization process, make that function of your stdev
r.min = 0.1 ; r.max = 0.2 ;
radius = r.min + (r.max-r.min).* rand(size(t)) ;

%// define surface and patch display options (FaceAlpha etc ...), for later
surfoptions  = {'FaceAlpha',0.2 , 'EdgeColor','none' , 'EdgeAlpha',0.1 , 'DiffuseStrength',1 , 'AmbientStrength',1 } ;
patchoptions = {'FaceAlpha',0.2 , 'EdgeColor','k'    , 'EdgeAlpha',0.2 , 'DiffuseStrength',1 , 'AmbientStrength',1 } ;
patchcol     = [1 0 0] ;  % Color of patch

%% // get the gradient at each point of the curve
Gx = diff([x,x(1)]).' ;                       %'//damn StackOverflow prettifier 
Gy = diff([y,y(1)]).' ;                       %'//damn StackOverflow prettifier 
Gz = diff([z,z(1)]).' ;                       %'//damn StackOverflow prettifier 
%// get the middle gradient between 2 segments (optional, just for better rendering if low number of points)
G = [ (Gx+circshift(Gx,1))./2 (Gy+circshift(Gy,1))./2 (Gz+circshift(Gz,1))./2] ;

%% // get the angles (azimuth, elevation) of each plane normal to the curve

ux = [1 0 0] ;
uy = [0 1 0] ;
uz = [0 0 1] ;

for k = nPts:-1:1 %// running the loop in reverse does automatic preallocation
    a = G(k,:) ./ norm(G(k,:)) ;
    angx(k) =  atan2( norm(cross(a,ux)) , dot(a,ux))  ;
    angy(k) =  atan2( norm(cross(a,uy)) , dot(a,uy))  ;
    angz(k) =  atan2( norm(cross(a,uz)) , dot(a,uz))  ;

    [az(k),el(k)] = cart2sph( a(1) , a(2) , a(3) ) ;
end
%// adjustment to be normal to cross section plane the way the rotation are defined later
az = az + pi/2 ; 
el = pi/2 - el ;

%% // define basic disc
discResolution = 20 ;
tt = linspace( 0 , 2*pi , discResolution ) ;
xd = cos(tt) ;
yd = sin(tt) ;
zd = zeros( size(xd) ) ;

%% // Generate coordinates for each cross section

ccylX = zeros( nPts , discResolution ) ;
ccylY = zeros( nPts , discResolution ) ;
ccylZ = zeros( nPts , discResolution ) ;
ccylC = zeros( nPts , discResolution ) ;

for ip = 1:nPts
    %// cross section coordinates, with radius function of [rand] in this
    %// example. Make it function of the stdev in your application.
    csTemp = [ ( radius(ip) .* xd )  ; ... %// X coordinates
               ( radius(ip) .* yd )  ; ... %// Y coordinates
                               zd    ] ;   %// Z coordinates

    %// rotate the cross section (around X axis, around origin)
    elev = el(ip) ;
    Rmat = [ 1     0           0     ; ...
             0 cos(elev)  -sin(elev) ; ...
             0 sin(elev)   cos(elev) ] ;
    csTemp = Rmat * csTemp ;

    %// do the same again to orient the azimuth (around Z axis)
    azi = az(ip) ;
    Rmat = [ cos(azi)  -sin(azi) 0 ; ...
             sin(azi)   cos(azi) 0 ; ...
               0            0    1 ] ;
    csTemp = Rmat * csTemp ;

    %// translate each cross section where it should be and store in global coordinate vector
    ccylX(ip,:) = csTemp(1,:) + x(ip) ;
    ccylY(ip,:) = csTemp(2,:) + y(ip) ;
    ccylZ(ip,:) = csTemp(3,:) + z(ip) ;
end

%% // Display the full cylinder
hd.cyl = surf( ccylX , ccylY , ccylZ , ccylC ) ;

%// use that if the graphic object already exist but you just want to update your data:
%// set( hd.cyl , 'XData',ccylX , 'YData',ccylY ,'ZData',ccylZ ) 

set( hd.cyl , surfoptions{:} )


%% // this is just to avoid displaying the patches in the next block
%// comment the "return" instruction or just execute next block if you want
%// to see the building cross sections as patches
return 

%% // display patches
hp = zeros(nPts,1) ;
for ip = 1:nPts
   hp(ip) = patch( ccylX(ip,:) , ccylY(ip,:) , ccylZ(ip,:) , patchcol ) ;
   set( hp(ip) , patchoptions{:} )
end

---

And this is just a quick zoomed view with the patches on (code rerun with a lower number of points, otherwise it quickly clutter the whole figure):

Answer 2

My Matlab version is old, but something like this will hopefully work for you too:

1. Draw the cloudy patches first as you did above (not the solid path yet).

2. Store the axis handle:

   patch_axis = gca;
   

3. Create a new, overlapping axis:

   path_axis = axes('position',get(patch_axis,'position'),'visible','off');
   

4. Draw the solid path line (in this new axis).

5. Link rotation and limits of the patch_axis (behind) to that of the path_axis (in front):

   set(rotate3d(path_axis),'ActionPostCallback', ...
   @(src,evt)set(patch_axis,{'view','xlim','ylim','zlim'}, ...
   get(path_axis,{'view','xlim','ylim','zlim'})));
   

6. If you are rotating your view manually then it should line up the two axes after your first adjustment and keep them lined up. But if you're setting the rotation and limits with a command, then you can either just include both axis handles (patch_axis & path_axis) in your set command or copy the settings after:

   set(patch_axis,{'view','xlim','ylim','zlim'}, ...
   get(path_axis,{'view','xlim','ylim','zlim'})
   

Note that to adjust axis properties (tick labels etc.), you'll need to do it to the visible patch_axis not the invisible path_axis. And if you want it to be interactive in real time with manual rotation, I'm not sure how to get the alignment function to execute with every redraw.