Matlab: "X-Ray" storyline via patch

Here is the implementation of the solution I mentioned in the comment (just a drawing of one cloud surface

).

It is not optimized, there are a few loops for

that could have been avoided by clever use of bsxfun or these helper function families, but it works fine as it is. The math for defining the tangent curve at each point and orienting (rotating) each cross section could also be simplified, but that's not my strong suit, so I'll leave that to the experts if they like it.

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

.

As long as the surface overlaps itself for a long time (depending on the viewing angle), there is still a little shading of the main center line, but the main line is always visible. Also, you only have one graphical object to manage.

The result looks like this:

Of course, you can change the AlphaValue

surfaces to your liking. I defined the full matrix the same size as the data for the color information. Everything is currently set to 0

(which is why it indicates green in the default color palette), but just as easy if you want to make a color function of another parameter, just adjust the color matrix accordingly (and the color palette that comes with it).

At the end of the code, there is an option to display each cross section as a patch object. It is not intended to be used in the end result, but it can help you understand how it all works if you want to make your own changes.

Here is 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 preview with patches (repeating code with fewer dots, otherwise it will quickly clutter up the whole drawing):