Matlab：通过补丁的“X-Ray”情节线

Question

问题

我正在尝试可视化 3-D 路径，以及它周围的“云”，它代表数据的标准偏差。我希望能够看到一条粗黑线作为路径，其周围有一个均匀的灰色区域，线条没有任何混浊，就好像看到这条线像 X 射线一样穿过云层。

尝试

我使用plot3 创建了一条粗线，并使用patch 创建了一系列围绕该线的每个点的框（图中还有一些额外的东西来表示开始/停止和方向，我也会喜欢他们很容易看到）。我尝试使用patch 中的alpha，但这会导致线条浑浊，因此无论视线中有多少灰色框，灰色的亮度都会发生变化。我希望alpha 为 1，这样每个灰色框的颜色都完全相同，但我希望找到某种方法可以使线条均匀地穿过云层。

最小示例

根据要求，这是一个最小的示例，它产生了下面的图。

% 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

对于 for 循环中非常长的代码段表示歉意，patch 命令有相当多的参数。前三行简单地定义了定义立方体的 8 个顶点的 x、y 和 z 坐标，方法是指定中心点加上或减去立方体的某个半宽 cloud(i)。其余的由各自的 cmets 解释。

感谢您的帮助！

Answer 1

这是我在评论中提到的解决方案的实现（只画了一个 surface 云）。

它没有优化，有一些for 循环可以通过巧妙地使用 bsxfun 或这些辅助函数系列来避免，但它运行正常。在每个点上找到曲线的切线并定位（旋转）每个横截面的数学可能也可以简化，但这不是我的强项，所以如果专家愿意，我会将其留给专家。

基本上，它定义了一个圆（在代码中通常称为“横截面”），其半径与某物（应用程序中的标准偏差，示例中的随机值）成比例。然后每个圆在 3D 中旋转，使其在平移点处与曲线垂直。然后将所有这些圆圈用作单个 surface 图形对象的信封。

当表面多次重叠时（取决于视角），主中心线仍有一点阴影，但主线始终可见。此外，您只需管理一个图形对象。

结果如下所示：

当然，您可以根据自己的喜好改变表面的AlphaValue。我定义了一个与颜色信息数据大小相同的完整矩阵。目前全部设置为0（所以它指向默认颜色图中的绿色），但是这种方式也很容易，如果你想制作另一个参数的颜色功能，只需相应地调整颜色矩阵（以及随之而来的颜色图）。

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

代码如下：

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

---

这只是一个带有补丁的快速缩放视图（代码以较少的点数重新运行，否则它很快就会使整个图变得混乱）：

Answer 2

我的 Matlab 版本很旧，但希望这样的东西也能对你有用：

1. 首先像上面那样绘制多云斑块（还不是实心路径）。

2. 存储轴句柄：

   patch_axis = gca;
   

3. 创建一个新的重叠轴：

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

4. 绘制实心路径线（在这个新轴上）。

5. 将 patch_axis（后面）的旋转和限制链接到 path_axis（前面）的旋转和限制：

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

6. 如果您手动旋转视图，那么它应该在您第一次调整后对齐两个轴并保持它们对齐。但是，如果您使用命令设置旋转和限制，那么您可以在 set 命令中包含两个轴句柄（patch_axis 和 path_axis），或者在之后复制设置：

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

请注意，要调整轴属性（刻度标签等），您需要对可见的 patch_axis 而不是不可见的 path_axis 进行调整。如果您希望它通过手动旋转实时交互，我不确定如何让对齐功能在每次重绘时执行。