Matlab下3D网格标量函数最慢增长方向的路径值提取方法

实现思路 #

• 从已定位的最小值点出发，采用三维网格贪心步进法：每一步在当前点的26个三维相邻网格点中，筛选满足「F值大于当前点F值」的候选点，选择其中F值最小的点作为下一个路径点，保证每一步的F增量最小，即沿最慢增长方向行进
• 记录每一步的坐标、F值，累计计算与起点的欧氏距离作为r坐标，直到路径点触及三维网格边界停止

完整Matlab实现代码 #

% ---------------- 前置逻辑（与你现有代码一致） ----------------
x = linspace(-10,80,100);
y = linspace(-20,5,100);
z = linspace(-10,10,100);
[X,Y,Z] = meshgrid(x,y,z);

F = some_scalar_function(X, Y, Z);

% 定位最小值点坐标
[minval,ind] = min(F(:));
[ii,jj,kk] = ind2sub(size(F),ind);
xmin = x(jj);
ymin = y(ii);
zmin = z(kk);

% ---------------- 新增最慢增长路径提取逻辑 ----------------
% 初始化路径存储变量
path_ijk = [ii,jj,kk]; % 路径点对应的网格下标
path_xyz = [xmin, ymin, zmin]; % 路径点对应的实际空间坐标
path_F = minval; % 路径点对应的F值
path_r = 0; % 路径点到起点的欧氏距离，即你需要的r坐标

% 生成三维网格26邻域的偏移量（覆盖所有相邻网格点）
[dx, dy, dz] = meshgrid(-1:1, -1:1, -1:1);
neighbor_offset = [dx(:), dy(:), dz(:)];
neighbor_offset(neighbor_offset(:,1)==0 & neighbor_offset(:,2)==0 & neighbor_offset(:,3)==0, :) = []; % 去掉当前点自身

max_step = 1000; % 配置最大步长防止死循环
for step = 1:max_step
    current_ijk = path_ijk(end,:);
    current_F = path_F(end);
    
    % 生成所有邻域点的网格下标
    neighbor_ijk = current_ijk + neighbor_offset;
    % 筛选落在网格有效范围内的邻域点
    valid = neighbor_ijk(:,1)>=1 & neighbor_ijk(:,1)<=size(F,1) ...
        & neighbor_ijk(:,2)>=1 & neighbor_ijk(:,2)<=size(F,2) ...
        & neighbor_ijk(:,3)>=1 & neighbor_ijk(:,3)<=size(F,3);
    neighbor_ijk = neighbor_ijk(valid,:);
    if isempty(neighbor_ijk)
        break % 走到网格边界，停止行进
    end
    
    % 提取所有有效邻域点的F值
    neighbor_F = arrayfun(@(i) F(neighbor_ijk(i,1), neighbor_ijk(i,2), neighbor_ijk(i,3)), 1:size(neighbor_ijk,1));
    % 筛选F值大于当前点的候选点（保证向外增长，不回落）
    candidate_idx = neighbor_F > current_F;
    if ~any(candidate_idx)
        break % 无有效增长方向，停止行进
    end
    candidate_ijk = neighbor_ijk(candidate_idx,:);
    candidate_F = neighbor_F(candidate_idx);
    
    % 选择F值最小的候选点作为下一个路径点，保证增长最慢
    [~, min_idx] = min(candidate_F);
    next_ijk = candidate_ijk(min_idx,:);
    next_xyz = [x(next_ijk(2)), y(next_ijk(1)), z(next_ijk(3))];
    next_F = candidate_F(min_idx);
    next_r = path_r(end) + norm(next_xyz - path_xyz(end,:));
    
    % 写入路径存储变量
    path_ijk = [path_ijk; next_ijk];
    path_xyz = [path_xyz; next_xyz];
    path_F = [path_F; next_F];
    path_r = [path_r; next_r];
end

% ---------------- 结果输出及可视化 ----------------
% 你需要的r坐标序列和对应F值序列
r = path_r;
F_along_r = path_F;

% 可视化等势面和提取的路径
figure;
isosurface(X,Y,Z,F,minval+100);
hold on;
plot3(path_xyz(:,1), path_xyz(:,2), path_xyz(:,3), 'r-', 'LineWidth',2);
xlabel('x'); ylabel('y'); zlabel('z');
legend('F等势面','最慢增长路径');

% 自定义标量函数（与你现有定义一致）
function F = some_scalar_function(X, Y, Z)
F = (X-6).^2 + 10*(Y+2).^2 + 10*Z.^2 + 5*X.*Y;
end

说明 #

• 最终输出的r为沿路径的累计距离坐标，F_along_r为对应位置的F值，可直接通过plot(r, F_along_r)查看F随r的变化曲线
• 如果需要更高精度的路径，可将贪心步进替换为基于梯度的连续追踪：每一步通过imgradient3获取当前点的梯度方向，取与梯度垂直的方向中F增长最小的方向作为行进方向，配合interp3插值获取非网格点的F值即可

内容的提问来源于stack exchange，提问作者teeeeee