Matlab：“ X射线”绘图线通过补丁

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).

代码末尾有一个选项，可将每个横截面显示为补丁对象。它不打算在最终结果中使用，但是如果您想进行自己的修改，它可以帮助您了解整个对象的构造方式。

代码如下：

%% // 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

这只是贴有补丁的快速缩放视图（代码以较少的点数重新运行，否则会迅速使整个图形混乱）：