matlab - 在 3D 网格上找到对应于最慢增加方向的标量函数的值

Question

我在 3D 空间的网格上定义了一些标量函数 F(x,y,z)，并且在数组中的某处有最小的 F。下面给出了生成这样一个函数并定位最小值坐标的示例代码：

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

% Find the minimum of the function on the grid
[minval,ind] = min(F(:));
[ii,jj,kk] = ind2sub(size(F),ind);
xmin = x(jj);
ymin = y(ii);
zmin = z(kk);

figure;isosurface(X,Y,Z,F,minval+100)

% Some sample scalar function (assume it is given on the grid, and the analytic form not known)
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 值的向量，该方向对应于函数 F增长最慢的方向，即从最小值开始并向外“行走”。我也想获得坐标 r 的对应值。我试图在下图中解释我的意思：

沿着 r 以外的任何方向走一条路径都会导致 F 的急剧增加，因此不是正确的路线。任何人都可以展示如何在 Matlab 中完成此操作吗？谢谢！

---

编辑

在 rahnema1 和 Ander Biguri 发表评论后，我运行了命令

[Gmag,Gazimuth,Gelevation] = imgradient3(F);

看看通过 z=0 的平面，函数 F(x,y,z=0) 本身如下所示：

并且输出imgradient3()看起来像这样（同样，从得到的完整 3D 阵列中只有一个平面）：

如何从这些中获得与最慢增长路径相对应的线切割作为 r 的函数？（仍然记住它们是 3D 阵列，方向不一定限于 z=0 平面）。

标签： matlaboptimizationvectorminimum

Answer 1

考虑一个球体，它的中心是 的最小点的位置F。对于球体表面上的每个点，计算到中心点的最短路径。用于imerode查找球体的表面并用于shortestpathtree计算最短路径树。梯度大小设置为路径遍历的难度。用于imgradient3计算梯度幅度。到中心距离最小的曲面点之间的路径path1，被认为是函数增长最慢的方向F。但是path1是路径的一半。另一半 , path2, 是通过消除包含第一个表面点的子树并使用结果图再次计算最短路径树来计算的。功能im2col3定义为将 3D 数组转换为图形数据结构。
这是未经测试的代码：

% im2col for 3D image
function col = im2col3 (im)
  col = zeros(prod(size(im)-2), 26);
  n = 1;
  for kk = 1:3
    for jj = 1:3
      for ii = 1:3
        if ~(ii == 2 && jj == 2 && kk == 2)
          col(:, n) = reshape(im(ii:end-(3-ii), jj:end-(3-jj), kk:end-(3-kk)),[],1);
          n = n + 1;
        end
      end
    end
  end
end

% gradient magnitude
[Gmag,~,~] = imgradient3(F);

% convert array to graph
idx = reshape(1:numel(F), size(F));
t = im2col3(idx);
w = Gmag(t);
s = repmat(reshape(idx(2:end-1,2:end-1,2:end-1),[],1), size(t, 2));
G = graph(s, t, w);

% form the sphere
[~, sphere_center] = min(F(:));
[r, c, z] = ind2sub(size(F), sphere_center);
sphere_radius_2 = min([r c z])  ^ 2;
sphere_logical = ((1:size(F, 1)).'- r) .^ 2 ...
               + ((1:size(F, 2))- c) .^ 2 ...
               + (reshape(1:size(F, 3), 1, 1, [])- z) .^ 2 ...
               < sphere_radius_2;

se = strel('cube',3);
sphere_surface = xor(imerode(sphere_logical, se), sphere_logical);
sphere_nodes = idx(sphere_surface);

% compute the half of the path
[T, D] = shortestpathtree(G, sphere_center, sphere_nodes,'OutputForm','cell');
[mn1, im1] = min(D);

path1 = T{im1};

% eliminate the sub-tree
subtree_root = path1(2);
subtree_nodes = bfsearch(T, subtree_root);
G1 = rmnodes (G, subtree_nodes);
sphere_nodes = setdiff (sphere_nodes, subtree_nodes);

% computing other half
[T1, D1] = shortestpathtree(G1, sphere_center, sphere_nodes,'OutputForm','cell');
[mn2, im2] = min(D1);

path2 = T1{im2};