如何检测三维图中"脊"的出现位置？

Question

请问如何检测得到下图中这些“脊（隆起）”在图中的位置。

产生该图像的代码是：

c1 = 3/5;
c2 = -3/10;
ck = 1/4;
k1 = 3;
k2 = 3/10;
H1[s_] := 1/(s^2 + c1*s + k1)
H2[s1_, s2_] := (2*k2 + ck (s1 + s2) + 2 c2*s1*s2)*H1[s1]*H1[s2]*
  H1[s1 + s2]/2
h2[t1_, t2_] = 
  InverseLaplaceTransform[H2[s1, s2], {s1, s2}, {t1, t2}];

Plot3D[h2[t1, t2], {t1, 0, 16}, {t2, 0, 16}, PlotRange -> All, 
 Exclusions -> None, PlotPoints -> {80, 80}, 
 ColorFunction -> "Rainbow"]

---

根据@苹果 的方法，得到的如下图：

我实际的离散数据如下图（数据如何上传？）：

我实际想得到的是像这样：

即想得到这些“山脊”的位置。苹果的方法得到的是孤立的点，实际上“山脊”上的点并不同时满足在两个方向的偏导为零。@Lozmlve*永

Answer 1

V0

目测是没有解析解的，那就数值求解吧。

solve[{x0_, y0_}] /; (y0 != x0) := 
  Quiet@Check[
    FindRoot[{D[h2[t1, t2], t1] == 0, 
      D[h2[t1, t2], t2] == 0}, {{t1, x0}, {t2, y0}}], {}];
solve[{x0_, y0_}] /; (y0 == x0) := {};
sol = Partition[
   Flatten@Table[solve[{x0, y0}], {x0, 0, 16, 0.2}, {y0, 0, 16, 0.2}],
    2];

p = Select[{t1, t2} /. DeleteDuplicates[sol], 
   0.01 <= #1 <= 16 && 0.01 <= #2 <= 16 & @@ # &];

这是找到了所有的导数为0的地方，还需要进行一次筛选，才能够找到极大值。

V1

对“脊”的定义不明确。所以随手写了一个。思路是t2固定时，寻找极大值位置的t1。

c1 = 3/5;
c2 = -3/10;
ck = 1/4;
k1 = 3;
k2 = 3/10;
H1[s_] := 1/(s^2 + c1*s + k1)
H2[s1_, s2_] := (2*k2 + ck (s1 + s2) + 2 c2*s1*s2)*H1[s1]*H1[s2]*H1[s1 + s2]/2
h2[t1_, t2_] = InverseLaplaceTransform[H2[s1, s2], {s1, s2}, {t1, t2}];

p1 = Plot3D[h2[t1, t2], {t1, 0, 16}, {t2, 0, 16}, PlotRange -> All, 
   Exclusions -> None, PlotPoints -> {80, 80}, 
   ColorFunction -> "Rainbow", ImageSize -> 100];
  
expr = h2[t1, t2] /. HeavisideTheta[_, _] -> 0;(*只考虑t2<=t1的区间*)

cf = With[{expr = expr}, Compile[{{t1, _Real}, {t2, _Real}}, expr]];
data = Table[cf[t1, t2], {t2, 0, 16, 0.01}, {t1, 0, 16, 0.01}]; 

pos = Monitor[
   Join @@ Table[{0.01 (First@# - 1), 0.01 (i - 1), Last@#} & /@ 
      FindPeaks[data[[i]]], {i, Length@data}], i];
pos2 = Select[pos, #1 >= #2 & @@ # &];

Show[p1, Graphics3D[{Red, PointSize[0.01], Point[pos2]}], 
 Graphics3D[{Blue, Point[{#, 0.5, cf[#, 0.5]} & /@ Range[0, 16, 0.1]]}], 
 ImageSize -> 300]

Comments

谢谢。我给出这个具有解析式的例子是为了方便说明，实际上我只有离散数据，当然也可以用差分替代你这里的解析求导。不知道从图像特征提取角度有没有更好的方法，我试了下RidgeFilter（脊滤波），不太成功。还是很感谢。

---

已更新V1。。。。。。。

Answer 2

产生你说的离散点

p = Table[{t1, t2, h2[t1, t2]}, {t1, 0, 16, .1}, {t2, 0, 16, .1}];

图像处理

MaxDetect[p // Image]

注意轴的方向，画图坐标和Graphics之间的区别

---

不是很清楚你具体要什么

p = Table[h2[t1, t2], {t1, 0, 16, .1}, {t2, 0, 16, .1}];
pic = Colorize[p // Image, ColorFunction -> GrayLevel];

这种？

r1 = MorphologicalBinarize[pic, {0.08, 0.15}]

还是这种？

r2 = Thinning@r1

Comments

float/double精度怎么处理。

---

@苹果 哦，这个看他的那个离散的程序，他在你的评论里说他的是离散数据

---

。。。之前没看到他的回复。

---

谢谢，请看我的补充内容。

---

那怎么定义“山脊”。。。

---

这确实是个问题，人一眼就能看明白，但如何将其进行数学表达，我也还在思考。

---

我之前是想借鉴CrossingDetect帮助文档中的第二个例子，但却不知如何将三维图像转化成这种二维图像。RidgeFilter帮助文档的"范围"栏中的第二个算例和我这个有类似之处，但似乎对于参数的选择很敏感。