10.20日志—claude的讲解文件关于code文件夹

1、告诉AI

分析code目录下工程,各个文件的调用关系,各个函数的意义和调用关系,show me the code for that,可以输出md文档形式,也可以输出带有图画的htmll形式(便于我理解

2、讲解文件如下:

点击查看代码
<!DOCTYPE html>
<html lang="zh-CN">
<head>
<meta charset="UTF-8">
<meta name="viewport" content="width=device-width, initial-scale=1.0">
<title>code/ 工程分析 — 文件、函数与调用关系</title>
<style>
  :root{
    --bg:#f1f5f9; --card:#ffffff; --ink:#0f172a; --muted:#64748b;
    --line:#e2e8f0; --blue:#2563eb; --purple:#7c3aed; --green:#16a34a;
    --amber:#d97706; --red:#dc2626; --nav:#0f172a; --navink:#cbd5e1;
  }
  *{box-sizing:border-box}
  html{scroll-behavior:smooth}
  body{margin:0;background:var(--bg);color:var(--ink);font-family:"Segoe UI","Microsoft YaHei","PingFang SC",system-ui,-apple-system,sans-serif;line-height:1.72;font-size:15px}
  .layout{display:grid;grid-template-columns:276px minmax(0,1fr);max-width:1560px;margin:0 auto}
  aside{position:sticky;top:0;height:100vh;overflow-y:auto;background:var(--nav);color:var(--navink);padding:26px 16px 40px}
  aside .brand{font-size:18px;font-weight:800;color:#fff;line-height:1.35;margin:0 8px 6px}
  aside .sub{font-size:12px;color:#94a3b8;margin:0 8px 22px}
  aside a{display:block;color:var(--navink);text-decoration:none;padding:7px 10px;border-radius:8px;font-size:13.5px;margin:1px 0}
  aside a:hover{background:#1e293b;color:#fff}
  aside .group{font-size:11px;letter-spacing:.12em;color:#64748b;text-transform:uppercase;margin:18px 10px 6px}
  main{padding:38px 50px 90px;min-width:0}
  h1{font-size:31px;line-height:1.25;margin:0 0 8px;letter-spacing:-.01em}
  .subtitle{color:var(--muted);font-size:15px;margin:0 0 26px}
  h2{font-size:25px;margin:54px 0 6px;padding-bottom:10px;border-bottom:2px solid var(--line);letter-spacing:-.01em}
  h3{font-size:19px;margin:28px 0 8px}
  h4{font-size:15.5px;margin:20px 0 6px}
  p{margin:10px 0}
  .card{background:var(--card);border:1px solid var(--line);border-radius:16px;padding:22px 26px;margin:18px 0;box-shadow:0 1px 2px rgba(15,23,42,.04)}
  .card.tight{padding:16px 20px}
  .lead{font-size:16.5px;color:#334155}
  .muted{color:var(--muted)}
  .small{font-size:13px}
  .badge{display:inline-block;font-size:11.5px;font-weight:700;padding:2px 9px;border-radius:999px;vertical-align:middle;letter-spacing:.02em;white-space:nowrap}
  .b-script{background:#dbeafe;color:#1e40af;border:1px solid #bfdbfe}
  .b-mex{background:#ede9fe;color:#5b21b6;border:1px solid #ddd6fe}
  .b-helper{background:#dcfce7;color:#166534;border:1px solid #bbf7d0}
  .b-display{background:#fef3c7;color:#92400e;border:1px solid #fde68a}
  .b-legacy{background:#f1f5f9;color:#475569;border:1px solid #cbd5e1}
  .b-builtin{background:#fce7f3;color:#9d174d;border:1px solid #fbcfe8}
  .b-unused{background:#fee2e2;color:#991b1b;border:1px solid #fecaca}
  table{width:100%;border-collapse:collapse;margin:14px 0;font-size:13.8px}
  th,td{border-bottom:1px solid var(--line);padding:9px 11px;text-align:left;vertical-align:top}
  th{background:#f8fafc;font-weight:700;color:#334155;position:relative}
  tr:hover td{background:#fbfdff}
  td code, p code, li code, .kv code{background:#eef2ff;color:#3730a3;border-radius:5px;padding:1px 6px;font-family:Consolas,"Cascadia Code","Courier New",monospace;font-size:12.8px;white-space:nowrap}
  pre{margin:0;background:#0b1220;color:#e2e8f0;padding:16px 18px;border-radius:0 0 10px 10px;overflow:auto;font-family:Consolas,"Cascadia Code","Courier New",monospace;font-size:12.6px;line-height:1.55;tab-size:4}
  pre code{background:none;color:inherit;padding:0;white-space:pre}
  details{margin-top:12px;border:1px solid var(--line);border-radius:10px;background:#fff}
  summary{cursor:pointer;padding:10px 14px;background:#f8fafc;font-weight:650;border-radius:10px;font-size:13.5px;color:#334155;user-select:none}
  summary:hover{background:#f1f5f9}
  details[open] summary{border-bottom:1px solid var(--line);border-radius:10px 10px 0 0}
  .sig{font-family:Consolas,"Cascadia Code","Courier New",monospace;font-size:13px;background:#f8fafc;border-left:4px solid var(--blue);padding:10px 14px;border-radius:0 8px 8px 0;overflow-x:auto;white-space:nowrap}
  .kv{display:grid;grid-template-columns:112px 1fr;gap:7px 14px;font-size:13.8px;margin:14px 0}
  .kv > div:nth-child(odd){color:var(--muted);font-weight:600}
  .fn{background:#fff;border:1px solid var(--line);border-radius:16px;padding:20px 24px;margin:20px 0;box-shadow:0 1px 2px rgba(15,23,42,.04)}
  .fn h3{margin:0 0 10px;display:flex;align-items:center;gap:10px;flex-wrap:wrap}
  .fn h3 .path{font-size:12.5px;color:var(--muted);font-weight:500;font-family:Consolas,monospace}
  .diagram{overflow-x:auto;background:#fff;border:1px solid var(--line);border-radius:16px;padding:14px;margin:18px 0;box-shadow:0 1px 2px rgba(15,23,42,.04)}
  .diagram svg{display:block;min-width:900px;width:100%;height:auto}
  .diagram .caption{font-size:12.5px;color:var(--muted);text-align:center;margin-top:8px}
  .equation{text-align:center;font-size:19px;padding:20px 16px;background:#f8fafc;border:1px solid var(--line);border-radius:12px;font-family:"Cambria Math","Times New Roman",serif;overflow-x:auto}
  .callout{border-left:4px solid var(--amber);background:#fffbeb;padding:12px 16px;border-radius:8px;margin:14px 0;font-size:14px}
  .callout.info{border-color:var(--blue);background:#eff6ff}
  .callout.danger{border-color:var(--red);background:#fef2f2}
  .callout.ok{border-color:var(--green);background:#f0fdf4}
  .callout b{display:block;margin-bottom:3px}
  .pipeline{display:grid;grid-template-columns:repeat(3,minmax(0,1fr));gap:12px;margin:16px 0;counter-reset:step}
  .step{position:relative;min-width:0;background:#fff;border:1px solid #cbd5e1;border-radius:12px;padding:11px 12px 11px 42px;font-size:12.8px}
  .step::before{counter-increment:step;content:counter(step);position:absolute;left:11px;top:11px;width:22px;height:22px;border-radius:50%;background:var(--blue);color:#fff;font-size:11.5px;font-weight:700;display:flex;align-items:center;justify-content:center}
  .step .fnname{font-family:Consolas,monospace;font-weight:700;color:#1d4ed8;font-size:12.8px;word-break:break-word}
  .step .shape{display:inline-block;background:#e0e7ff;color:#3730a3;border-radius:999px;padding:1px 8px;font-size:11.5px;margin:4px 0}
  .step .note{color:var(--muted);font-size:12px;line-height:1.45}
  .arrow{display:none}
  .grid2{display:grid;grid-template-columns:repeat(2,minmax(0,1fr));gap:16px}
  .grid3{display:grid;grid-template-columns:repeat(3,minmax(0,1fr));gap:14px}
  .mini{background:#f8fafc;border:1px solid var(--line);border-radius:12px;padding:12px 14px;font-size:13px}
  .mini b{color:#334155}
  .toc-note{font-size:12.5px;color:var(--muted);margin:8px 8px 0}
  .top{position:fixed;right:22px;bottom:22px;background:#0f172a;color:#fff;border-radius:999px;padding:9px 15px;text-decoration:none;font-size:12.5px;box-shadow:0 6px 18px rgba(15,23,42,.25);opacity:.85}
  .top:hover{opacity:1}
  .filetree{background:#0b1220;color:#cbd5e1;border-radius:12px;padding:16px 18px;font-family:Consolas,"Cascadia Code",monospace;font-size:12.8px;line-height:1.65;overflow-x:auto;white-space:pre}
  .filetree .c{color:#64748b}
  .filetree .k{color:#93c5fd}
  .kbd{font-family:Consolas,monospace;background:#f1f5f9;border:1px solid #cbd5e1;border-bottom-width:2px;border-radius:5px;padding:0 6px;font-size:12px}
  @media (max-width:1050px){
    .layout{grid-template-columns:1fr}
    aside{position:relative;height:auto;display:flex;flex-wrap:wrap;gap:2px;align-items:center;padding:14px}
    aside .brand,aside .sub{width:100%}
    aside .group{display:none}
    aside a{display:inline-block}
    main{padding:26px}
    .grid2,.grid3,.pipeline{grid-template-columns:1fr}
  }
  @media print{
    aside,.top{display:none}
    .layout{display:block}
    main{padding:0}
    .card,.fn,.diagram{break-inside:avoid;box-shadow:none}
    details{display:block}
    details > summary{display:none}
    details > pre{display:block}
  }
</style>
</head>
<body>
<div class="layout">
<aside>
  <div class="brand">code/ 工程分析</div>
  <div class="sub">文件 · 函数 · 调用关系 · 源码</div>
  <div class="group">总览</div>
  <a href="#overview">一、工程概览</a>
  <a href="#files">二、文件清单与职责</a>
  <a href="#callgraph">三、调用关系图</a>
  <a href="#dataflow">四、数据流水线</a>
  <a href="#loop">五、水平集演化循环</a>
  <div class="group">函数详解</div>
  <a href="#functions">六、函数逐个详解</a>
  <a href="#mex">七、MEX / C++ 内部结构</a>
  <div class="group">脚本与注意事项</div>
  <a href="#variants">八、demo 脚本对比</a>
  <a href="#legacy">九、遗留与未使用代码</a>
  <a href="#gotchas">十、坑点与注意事项</a>
  <a href="#run">十一、运行速查</a>
</aside>

<main>
<h1>code/ 工程分析</h1>
<p class="subtitle">基于字典的纹理图像分割 · 文件职责、函数签名、调用关系与源码对照<br>
生成日期:2026-10-19 · 分析对象:<code>E:\xbb\Workspace\MATLAB\ZiDianXueXi\code</code></p>

<section id="overview">
<h2>一、工程概览</h2>
<div class="card">
  <p class="lead"><b>一句话:</b>先用图像自身训练一棵 k-means 纹理字典树,把每个像素的 <i>M×M</i> 邻域翻译成字典词;再用双邻接矩阵把“当前区域标签”在字典空间与像素空间之间来回扩散,得到纹理概率图;最后用两个高斯模型 + 曲率项驱动水平集 φ 演化,φ=0 的等值线就是分割边界。</p>
  <p>这是 Dahl &amp; Dahl, <i>Dictionary Based Image Segmentation</i>(SCIA 2015)思路的 MATLAB 实现,不是深度学习工程。整个工程没有训练集、没有模型文件,算法所需的“字典”完全来自当前输入图像。</p>
  <div class="grid3">
    <div class="mini"><b>① 纹理表示</b><br><code>build_km_tree</code> + <code>search_km_tree</code> 把图像变成“字典树 + 每像素词号 A”。</div>
    <div class="mini"><b>② 概率扩散</b><br><code>biadjacency_matrix</code> + <code>transition_matrix</code> 给出 <code>T1</code>、<code>T2</code>,用 <code>T2*(T1*mask)</code> 把区域标签按纹理相似度扩散。</div>
    <div class="mini"><b>③ 曲线演化</b><br>加权拟合前景/背景两个高斯,计算对数似然比;叠加曲率(长度项)与面积项,迭代更新 φ。</div>
  </div>
  <div class="callout info"><b>依赖环境</b>MATLAB R2025a(本机路径 <code>E:\xbb\ProgramFiles\MATLAB\R2025a\bin\matlab.exe</code>)+ 3 个自编译 MEX(<code>build_km_tree</code>、<code>search_km_tree</code>、<code>biadjacency_matrix</code>)。旧版 <code>.mexw64</code> 已在本机 R2025a 验证可加载运行;换 MATLAB 版本后若报 “Invalid MEX-file”,在 <code>texture_functions/</code> 下运行 <code>compile_mex_functions</code> 重编。</div>
</div>

<h3>目录树(带职责标注)</h3>
<div class="filetree"><span class="k">code/</span>
├─ <span class="k">demo.m</span>                  <span class="c"># 主脚本:手动 roipoly 初始化,100 步</span>
├─ <span class="k">demo_v2.m</span> / <span class="k">demo1.m</span>       <span class="c"># 主脚本:圆心/半径自动初始化,100 步(两者仅格式差异)</span>
├─ <span class="k">demo3.m</span>                  <span class="c"># 变体:test_case=13,21 步,额外显示概率图 P</span>
├─ demo(SSCE-…).m            <span class="c"># legacy:initial_mask + mask2sdf 初始化,10 步</span>
├─ demo(ZGC-…).m             <span class="c"># legacy:同上,20 步</span>
├─ <span class="k">load_test_case.m</span>         <span class="c"># 测试用例参数表(0–13)</span>
├─ load_test_case(ZGC-…).m   <span class="c"># 旧参数表(支持 multiphase 截断)</span>
├─ load_test_case.asv        <span class="c"># 编辑器自动备份,非运行代码</span>
├─ <span class="k">deterministicKmeans.m</span>    <span class="c"># 确定性 k-means + 高斯拟合(initial≠0 时可选的初始化)</span>
├─ <span class="k">convertJxN.m</span>             <span class="c"># M×N×J → J×MN(根目录与 auxiliary_functions 各一份)</span>
├─ <span class="k">texture_functions/</span>
│  ├─ build_km_tree.cpp / .mexw64    <span class="c"># 建字典树</span>
│  ├─ search_km_tree.cpp / .mexw64   <span class="c"># 像素 → 字典词</span>
│  ├─ biadjacency_matrix.cpp / .mexw64 <span class="c"># 像素-字典双邻接矩阵 B</span>
│  ├─ transition_matrix.m            <span class="c"># B → T1, T2</span>
│  ├─ compile_mex_functions.m        <span class="c"># 编译上面 3 个 MEX</span>
│  └─ .idea/                         <span class="c"># IDE 工程文件,可忽略</span>
├─ <span class="k">level_set_functions/</span>
│  ├─ curvature.m            <span class="c"># 曲率 κ(φ)</span>
│  ├─ gaussianValue.m        <span class="c"># 多维高斯概率密度</span>
│  ├─ initial_mask.m         <span class="c"># legacy:圆盘并集掩膜</span>
│  ├─ mask2sdf.m             <span class="c"># legacy:掩膜 → 符号距离场</span>
│  ├─ relative_probability.m <span class="c"># 多相相对概率(当前未调用)</span>
│  ├─ show_contour.m         <span class="c"># 初始化图显示(带 cla)</span>
│  └─ show_phi.m             <span class="c"># legacy:φ 的 2D 热图</span>
└─ <span class="k">auxiliary_functions/</span>
   ├─ showContour.m          <span class="c"># 每轮迭代显示(无 cla)</span>
   ├─ showLevleSetFunction.m <span class="c"># 最终 φ 的 3D 曲面</span>
   ├─ blue_white_red.m       <span class="c"># 蓝白红色图(被 show_phi 使用)</span>
   └─ convertJxN.m           <span class="c"># 与根目录重复</span>
</div>
</section>

<section id="files">
<h2>二、文件清单与职责</h2>
<div class="card tight">
<table>
<thead><tr><th style="width:27%">文件</th><th style="width:11%">类型</th><th>作用</th><th style="width:20%">被谁调用 / 使用</th></tr></thead>
<tbody>
<tr><td><code>demo.m</code></td><td><span class="badge b-script">脚本</span></td><td>完整主流程;<code>roipoly</code> 手动圈选初始区域;100 步演化</td><td>用户手动运行</td></tr>
<tr><td><code>demo_v2.m</code> / <code>demo1.m</code></td><td><span class="badge b-script">脚本</span></td><td>完整主流程;按 <code>center</code>/<code>radius</code> 自动圆初始化;100 步(两者仅换行差异)</td><td>用户手动运行</td></tr>
<tr><td><code>demo3.m</code></td><td><span class="badge b-script">脚本</span></td><td><code>test_case=13</code>;循环硬编码 21 步;每轮额外 <code>imshow(P)</code>;结尾画 <code>log(P/(1-P))</code></td><td>用户手动运行</td></tr>
<tr><td><code>demo(SSCE-…).m</code></td><td><span class="badge b-legacy">脚本 legacy</span></td><td><code>test_case=3</code>;<code>initial_mask</code> + <code>mask2sdf</code> 初始化;10 步;启用面积归一化 <code>alpha</code>;结尾 <code>show_phi</code></td><td>旧实验脚本</td></tr>
<tr><td><code>demo(ZGC-…).m</code></td><td><span class="badge b-legacy">脚本 legacy</span></td><td><code>test_case=5</code>;同上;20 步;有“空区域”保护</td><td>旧实验脚本</td></tr>
<tr><td><code>load_test_case.m</code></td><td><span class="badge b-helper">函数</span></td><td>按 <code>test_case</code> 返回图像路径、圆心、半径、<code>patch_size</code>、<code>epsilon</code>、<code>gamma</code>、<code>beta</code>、<code>nr_step</code></td><td>所有 demo 脚本</td></tr>
<tr><td><code>load_test_case(ZGC-…).m</code></td><td><span class="badge b-legacy">函数 legacy</span></td><td>旧参数表;默认 <code>epsilon=0.01</code>、<code>beta=1</code>;末尾把多相 center/radius 截断为第一个圆</td><td>两个 legacy demo</td></tr>
<tr><td><code>load_test_case.asv</code></td><td><span class="badge b-legacy">备份</span></td><td>MATLAB 编辑器自动保存的旧版本,不参与运行</td><td>—</td></tr>
<tr><td><code>deterministicKmeans.m</code></td><td><span class="badge b-helper">函数</span></td><td>用确定性方式初始化 k-means,返回质心、权重、归属、协方差;可用于自动初始化</td><td>demo 的 <code>initial≠0</code> 分支(默认不执行)</td></tr>
<tr><td><code>convertJxN.m</code>(两处)</td><td><span class="badge b-helper">函数</span></td><td>把 <code>M×N×J</code> 图像转成 <code>J×MN</code>,每列一个像素</td><td>所有 demo 的彩色图分支</td></tr>
<tr><td><code>build_km_tree.cpp/.mexw64</code></td><td><span class="badge b-mex">MEX</span></td><td>随机采样 <code>M×M</code> 块,逐层 k-means,构建分叉字典树</td><td>所有 demo</td></tr>
<tr><td><code>search_km_tree.cpp/.mexw64</code></td><td><span class="badge b-mex">MEX</span></td><td>对每个像素的邻域块,从树根贪心下降到叶子,输出词号图 A</td><td>所有 demo</td></tr>
<tr><td><code>biadjacency_matrix.cpp/.mexw64</code></td><td><span class="badge b-mex">MEX</span></td><td>把“像素 ↔ (词, 块内位置)”关系写成稀疏双邻接矩阵 B</td><td>所有 demo</td></tr>
<tr><td><code>transition_matrix.m</code></td><td><span class="badge b-helper">函数</span></td><td>由 B 计算标签→字典的 <code>T1</code> 和字典→像素概率的 <code>T2</code></td><td>所有 demo</td></tr>
<tr><td><code>compile_mex_functions.m</code></td><td><span class="badge b-script">脚本</span></td><td>依次 <code>mex</code> 三个 C++ 文件</td><td>手动运行(换 MATLAB 版本时)</td></tr>
<tr><td><code>curvature.m</code></td><td><span class="badge b-helper">函数</span></td><td>用归一化梯度的散度近似水平集曲率 <code>κ(φ)</code></td><td>演化循环(每轮)</td></tr>
<tr><td><code>gaussianValue.m</code></td><td><span class="badge b-helper">函数</span></td><td>计算每个像素在给定高斯 <code>N(m,C)</code> 下的概率密度</td><td>演化循环(每轮两次)</td></tr>
<tr><td><code>initial_mask.m</code></td><td><span class="badge b-legacy">函数 legacy</span></td><td>生成一个或多个圆盘的并集逻辑掩膜</td><td>两个 legacy demo</td></tr>
<tr><td><code>mask2sdf.m</code></td><td><span class="badge b-legacy">函数 legacy</span></td><td>掩膜 → 符号距离场(注意符号约定,见坑点)</td><td>两个 legacy demo</td></tr>
<tr><td><code>relative_probability.m</code></td><td><span class="badge b-unused">未使用</span></td><td>多相水平集的相对概率 <code>P_i / (P_i + max_{j≠i} P_j)</code></td><td>当前没有任何脚本调用</td></tr>
<tr><td><code>show_contour.m</code></td><td><span class="badge b-display">显示</span></td><td>初始化时显示底图 + φ=0 轮廓;会先 <code>cla</code></td><td>所有 demo(初始化一次)</td></tr>
<tr><td><code>showContour.m</code></td><td><span class="badge b-display">显示</span></td><td>每轮显示底图 + 当前 φ=0 轮廓;<b>不</b> <code>cla</code></td><td>所有 demo(每轮)</td></tr>
<tr><td><code>showLevleSetFunction.m</code></td><td><span class="badge b-display">显示</span></td><td>把 φ 画成 3D 曲面,并叠加红色零等值线</td><td>所有 demo(结尾)</td></tr>
<tr><td><code>show_phi.m</code></td><td><span class="badge b-legacy">显示 legacy</span></td><td>φ 的 2D 热图 + 绿色零等值线</td><td>两个 legacy demo</td></tr>
<tr><td><code>blue_white_red.m</code></td><td><span class="badge b-display">显示</span></td><td>蓝-白-红色图生成</td><td>仅 <code>show_phi.m</code></td></tr>
</tbody>
</table>
</div>
</section>

<section id="callgraph">
<h2>三、调用关系图</h2>
<p>所有 <code>demo*.m</code> 都是<b>脚本</b>(不是函数),它们直接按顺序调用下面这些函数。下图按“入口脚本 → 三组函数”组织;虚线表示仅 legacy 脚本或可选分支使用。</p>
<div class="diagram">
<svg viewBox="0 0 1120 690" role="img" aria-label="code 目录调用关系图">
  <defs>
    <marker id="arrow" viewBox="0 0 10 10" refX="9" refY="5" markerWidth="7" markerHeight="7" orient="auto-start-reverse">
      <path d="M0,0 L10,5 L0,10 z" fill="#64748b"/>
    </marker>
    <marker id="arrowLegacy" viewBox="0 0 10 10" refX="9" refY="5" markerWidth="7" markerHeight="7" orient="auto-start-reverse">
      <path d="M0,0 L10,5 L0,10 z" fill="#94a3b8"/>
    </marker>
  </defs>
  <style>
    .node{fill:#fff;stroke:#94a3b8;stroke-width:1.4;rx:10}
    .node.script{fill:#dbeafe;stroke:#3b82f6}
    .node.mex{fill:#ede9fe;stroke:#8b5cf6}
    .node.helper{fill:#dcfce7;stroke:#22c55e}
    .node.display{fill:#fef3c7;stroke:#f59e0b}
    .node.legacy{fill:#f8fafc;stroke:#94a3b8;stroke-dasharray:5 4}
    .node .t{font-size:13.5px;font-weight:700;fill:#0f172a}
    .node .s{font-size:11px;fill:#64748b}
    .panel{fill:#f8fafc;stroke:#e2e8f0;stroke-width:1.5;rx:14}
    .panelTitle{font-size:14px;font-weight:700;fill:#334155}
    .edge{stroke:#64748b;stroke-width:1.5;fill:none}
    .edge.dash{stroke:#94a3b8;stroke-dasharray:6 5}
    .lbl{font-size:11px;fill:#64748b}
    .bus{stroke:#94a3b8;stroke-width:1.2;fill:none}
  </style>

  <!-- panels -->
  <rect class="panel" x="20" y="160" width="530" height="250"/>
  <rect class="panel" x="570" y="160" width="530" height="250"/>
  <rect class="panel" x="20" y="445" width="1080" height="175"/>

  <!-- entry -->
  <rect class="node script" x="380" y="20" width="360" height="72"/>
  <text class="t" x="560" y="46" text-anchor="middle">入口脚本(Scripts)</text>
  <text class="s" x="560" y="65" text-anchor="middle">demo.m · demo_v2.m · demo1.m · demo3.m</text>
  <text class="s" x="560" y="82" text-anchor="middle">demo(SSCE…).m · demo(ZGC…).m(legacy)</text>

  <!-- buses -->
  <path class="bus" d="M560 92 V132"/>
  <path class="bus" d="M30 132 H1090"/>
  <path class="bus" d="M30 132 V270 H550"/>

  <!-- row 1 drops -->
  <path class="edge" d="M125 132 V198" marker-end="url(#arrow)"/>
  <path class="edge" d="M300 132 V198" marker-end="url(#arrow)"/>
  <path class="edge" d="M465 132 V198" marker-end="url(#arrow)"/>
  <!-- panel 2 drops -->
  <path class="edge" d="M680 132 V198" marker-end="url(#arrow)"/>
  <path class="edge" d="M865 132 V198" marker-end="url(#arrow)"/>
  <path class="edge" d="M1015 132 V198" marker-end="url(#arrow)"/>
  <!-- row 2 drops (second bus) -->
  <path class="edge" d="M120 270 V288" marker-end="url(#arrow)"/>
  <path class="edge" d="M310 270 V288" marker-end="url(#arrow)"/>
  <path class="edge" d="M510 270 V288" marker-end="url(#arrow)"/>
  <path class="edge" d="M405 270 V358" marker-end="url(#arrow)"/>
  <!-- panel 3 feed -->
  <path class="edge" d="M560 132 V443" marker-end="url(#arrow)"/>

  <!-- panel 1 boxes -->
  <rect class="node helper" x="35" y="200" width="180" height="46"/>
  <text class="t" x="125" y="219" text-anchor="middle">load_test_case</text>
  <text class="s" x="125" y="235" text-anchor="middle">用例参数表(路径/半径)</text>

  <rect class="node helper" x="230" y="200" width="140" height="46"/>
  <text class="t" x="300" y="219" text-anchor="middle">convertJxN</text>
  <text class="s" x="300" y="235" text-anchor="middle">彩色图 → X</text>

  <rect class="node helper" x="385" y="200" width="160" height="46"/>
  <text class="t" x="465" y="219" text-anchor="middle">transition_matrix</text>
  <text class="s" x="465" y="235" text-anchor="middle">B → T1, T2</text>

  <rect class="node mex" x="40" y="290" width="160" height="46"/>
  <text class="t" x="120" y="309" text-anchor="middle">build_km_tree</text>
  <text class="s" x="120" y="325" text-anchor="middle">MEX · 建字典树</text>

  <rect class="node mex" x="230" y="290" width="160" height="46"/>
  <text class="t" x="310" y="309" text-anchor="middle">search_km_tree</text>
  <text class="s" x="310" y="325" text-anchor="middle">MEX · 像素 → 词</text>

  <rect class="node mex" x="420" y="290" width="180" height="46"/>
  <text class="t" x="510" y="309" text-anchor="middle">biadjacency_matrix</text>
  <text class="s" x="510" y="325" text-anchor="middle">MEX · 稀疏矩阵 B</text>

  <rect class="node legacy" x="280" y="360" width="260" height="46"/>
  <text class="t" x="410" y="379" text-anchor="middle">deterministicKmeans</text>
  <text class="s" x="410" y="395" text-anchor="middle">可选:initial ≠ 0 时初始化</text>

  <!-- panel 2 boxes -->
  <rect class="node helper" x="600" y="200" width="160" height="46"/>
  <text class="t" x="680" y="219" text-anchor="middle">gaussianValue</text>
  <text class="s" x="680" y="235" text-anchor="middle">N(x; m, C)</text>

  <rect class="node helper" x="790" y="200" width="150" height="46"/>
  <text class="t" x="865" y="219" text-anchor="middle">curvature</text>
  <text class="s" x="865" y="235" text-anchor="middle">κ(φ)</text>

  <rect class="node display" x="960" y="200" width="110" height="46"/>
  <text class="t" x="1015" y="219" text-anchor="middle">showContour</text>
  <text class="s" x="1015" y="235" text-anchor="middle">每轮轮廓</text>

  <rect x="600" y="280" width="470" height="92" rx="10" fill="#fff" stroke="#e2e8f0"/>
  <text class="s" x="616" y="301" fill="#334155">循环体(demo.m 第 84–123 行):</text>
  <text class="s" x="616" y="320">in = phi &gt; 0 → DictProb = T1*in(:) → p = T2*DictProb</text>
  <text class="s" x="616" y="339">pr = [p; 1−p] → 加权拟合两个高斯 → log_tr</text>
  <text class="s" x="616" y="358">δ = (1/ε)·p·(1−p) → φ += Δt·δ·(log_tr + β·κ − γ)</text>

  <!-- panel 3 boxes -->
  <rect class="node display" x="40" y="480" width="140" height="44"/>
  <text class="t" x="110" y="498" text-anchor="middle">roipoly</text>
  <text class="s" x="110" y="514" text-anchor="middle">手动初始化(demo.m)</text>

  <rect class="node legacy" x="200" y="480" width="150" height="44"/>
  <text class="t" x="275" y="498" text-anchor="middle">initial_mask</text>
  <text class="s" x="275" y="514" text-anchor="middle">legacy 圆盘掩膜</text>

  <rect class="node legacy" x="370" y="480" width="140" height="44"/>
  <text class="t" x="440" y="498" text-anchor="middle">mask2sdf</text>
  <text class="s" x="440" y="514" text-anchor="middle">legacy 掩膜 → SDF</text>

  <rect class="node display" x="530" y="480" width="150" height="44"/>
  <text class="t" x="605" y="498" text-anchor="middle">show_contour</text>
  <text class="s" x="605" y="514" text-anchor="middle">初始化显示(cla)</text>

  <rect class="node display" x="700" y="480" width="210" height="44"/>
  <text class="t" x="805" y="498" text-anchor="middle">showLevleSetFunction</text>
  <text class="s" x="805" y="514" text-anchor="middle">最终 φ 的 3D 曲面</text>

  <rect class="node legacy" x="930" y="480" width="130" height="44"/>
  <text class="t" x="995" y="498" text-anchor="middle">show_phi</text>
  <text class="s" x="995" y="514" text-anchor="middle">legacy 2D 热图</text>

  <rect class="node legacy" x="930" y="550" width="130" height="44"/>
  <text class="t" x="995" y="568" text-anchor="middle">blue_white_red</text>
  <text class="s" x="995" y="584" text-anchor="middle">蓝白红色图</text>

  <!-- legacy internal arrows -->
  <path class="edge dash" d="M350 502 H370" marker-end="url(#arrowLegacy)"/>
  <path class="edge dash" d="M995 524 V550" marker-end="url(#arrowLegacy)"/>

  <!-- panel titles -->
  <text class="panelTitle" x="30" y="432">① 准备与纹理建模</text>
  <text class="panelTitle" x="580" y="432">② 水平集演化(每轮迭代)</text>
  <text class="panelTitle" x="30" y="468">③ 初始化与可视化</text>

  <!-- legend -->
  <path class="edge" d="M40 650 H100" marker-end="url(#arrow)"/>
  <text class="lbl" x="110" y="654">调用 / 数据依赖</text>
  <path class="edge dash" d="M250 650 H310" marker-end="url(#arrowLegacy)"/>
  <text class="lbl" x="320" y="654">仅 legacy 脚本或可选分支</text>
</svg>
<div class="caption">图 1 函数调用关系。脚本直接调用函数;同一层中的函数之间没有调用关系(数据依赖见下一节)。</div>
</div>

<div class="callout"><b>读图要点</b>
<ul style="margin:6px 0 0 18px;padding:0">
  <li>入口是 5 个脚本;其中 <code>demo.m</code>、<code>demo_v2.m</code>、<code>demo1.m</code>、<code>demo3.m</code> 是当前版本,两个带括号文件名的是 legacy 实验脚本。</li>
  <li><code>load_test_case</code> 只提供参数;真正的图像读取是脚本里的 <code>imread</code>。</li>
  <li>三个 MEX 是“纹理表示”的核心;<code>transition_matrix</code> 把 MEX 的输出接成可重复使用的扩散算子。</li>
  <li><code>gaussianValue</code>、<code>curvature</code> 只在迭代循环里调用;<code>showContour</code> 也是每轮调用。</li>
  <li><code>initial_mask</code>/<code>mask2sdf</code>/<code>show_phi</code>/<code>blue_white_red</code> 只被 legacy 脚本使用;<code>relative_probability</code> 完全没有被调用。</li>
</ul></div>
</section>
<section id="dataflow">
<h2>四、数据流水线</h2>
<p>主流程是严格的单向流水线;唯一的反馈是水平集演化循环(下一节)。下图按编号顺序把每一步的函数、输出变量和形状放在一起。</p>

<div class="pipeline">
  <div class="step"><div class="fnname">imread</div><div class="shape">im: r×c×l</div><div class="note">读取 <code>test_images/…</code>,uint8</div></div>
  <div class="arrow">→</div>
  <div class="step"><div class="fnname">double / convertJxN</div><div class="shape">im_double / X</div><div class="note">X 为 d×N,每列一个像素;灰度图 X 为 1×N</div></div>
  <div class="arrow">→</div>
  <div class="step"><div class="fnname">build_km_tree</div><div class="shape">tree: (d·M²)×K</div><div class="note">默认 27×780;K = Σ b<sup>i</sup></div></div>
  <div class="arrow">→</div>
  <div class="step"><div class="fnname">search_km_tree</div><div class="shape">A: r×c</div><div class="note">词号 1…K;边界 (M−1)/2 为 0</div></div>
  <div class="arrow">→</div>
  <div class="step"><div class="fnname">biadjacency_matrix</div><div class="shape">B: N×(M²K)</div><div class="note">稀疏 logical;N = r·c</div></div>
  <div class="arrow">→</div>
  <div class="step"><div class="fnname">transition_matrix</div><div class="shape">T1, T2</div><div class="note">T1: (M²K)×N;T2: N×(M²K)</div></div>
  <div class="arrow">→</div>
  <div class="step"><div class="fnname">初始化 φ</div><div class="shape">φ: r×c</div><div class="note">手动 roipoly 或圆心/半径 → ±2</div></div>
  <div class="arrow">→</div>
  <div class="step"><div class="fnname">演化循环</div><div class="shape">φ: r×c</div><div class="note">gaussianValue + curvature,nr_step 次</div></div>
  <div class="arrow">→</div>
  <div class="step"><div class="fnname">showLevleSetFunction</div><div class="shape">φ = 0</div><div class="note">最终分割边界(零等值线)</div></div>
</div>

<h3>关键变量的形状与含义</h3>
<div class="card tight">
<table>
<thead><tr><th style="width:16%">变量</th><th style="width:22%">形状 / 类型</th><th>含义</th></tr></thead>
<tbody>
<tr><td><code>im</code></td><td><code>r×c×l</code>,uint8</td><td>原图;<code>l=1</code> 灰度,<code>l=3</code> 彩色</td></tr>
<tr><td><code>im_double</code></td><td><code>r×c×l</code>,double</td><td>MEX 输入要求 double(0–255 量级)</td></tr>
<tr><td><code>X</code></td><td><code>d×N</code>,double</td><td>每列一个像素的通道值;<code>d=l</code>,<code>N=r·c</code></td></tr>
<tr><td><code>tree</code></td><td><code>(d·M²)×K</code>,double</td><td>字典原子;每一列是一个节点中心,默认 <code>27×780</code></td></tr>
<tr><td><code>A</code></td><td><code>r×c</code>,double</td><td>每个像素的字典词编号 <code>1…K</code>;图像边界 <code>(M−1)/2</code> 像素为 0</td></tr>
<tr><td><code>B</code></td><td><code>N×(M²K)</code>,sparse logical</td><td>“像素 ↔ (词, 块内位置)”双邻接矩阵</td></tr>
<tr><td><code>T1</code></td><td><code>(M²K)×N</code>,sparse double</td><td>标签图 → 字典空间;<code>DictProb = T1*in(:)</code></td></tr>
<tr><td><code>T2</code></td><td><code>N×(M²K)</code>,sparse double</td><td>字典空间 → 像素概率;<code>p = T2*DictProb</code></td></tr>
<tr><td><code>phi</code></td><td><code>r×c</code>,double</td><td>水平集:<code>φ&gt;0</code> 内部,<code>φ&lt;0</code> 外部,<code>φ=0</code> 是边界</td></tr>
<tr><td><code>pr</code></td><td><code>2×N</code>,double</td><td>第 1 行前景软权重 <code>p</code>,第 2 行背景 <code>1−p</code></td></tr>
<tr><td><code>m</code></td><td><code>d×2</code>,double</td><td>前景/背景的加权平均颜色</td></tr>
<tr><td><code>covar</code></td><td><code>d×d×2</code>,double</td><td>前景/背景的加权协方差矩阵</td></tr>
<tr><td><code>log_tr</code></td><td><code>r×c</code>,double</td><td>每个像素的对数似然比:更像前景还是背景</td></tr>
<tr><td><code>delta</code></td><td><code>r×c</code>,double</td><td>边界带权重(Dirac δ 的近似),只在 <code>p≈0.5</code> 附近较大</td></tr>
<tr><td><code>P</code></td><td><code>r×c×nr_step</code></td><td>每轮保存的概率图;当前只有 <code>demo3.m</code> 用来显示</td></tr>
</tbody>
</table>
</div>

<h3>为什么 <code>T2*(T1*mask)</code> 能扩散区域标签</h3>
<div class="card">
<p><code>B</code> 的构造方式(C++ 中 0-based 列号)是:对每个内部像素 <code>ic</code>,把它 <code>M×M</code> 邻域中的每个像素 <code>i</code>,连到“<code>ic</code> 的词号 + 该像素在块内的位置”这一列:</p>
<div class="sig">j = (c+dx) + (c+dy)*M + (A[ic]−1)*M*M&nbsp;&nbsp;&nbsp;&nbsp;c = (M−1)/2</div>
<p>于是每一列代表一个“字典词 + 块内相对位置”。<code>transition_matrix</code> 做的两步归一化是:</p>
<div class="equation">T1 = D<sub>y</sub><sup>−1</sup> B<sup>T</sup>,&nbsp;&nbsp; T2 = D<sub>x</sub><sup>−1</sup> B</div>
<p>其中 <code>Dxinv</code> 是 <code>1/(B 的行和+eps)</code>,<code>Dyinv</code> 是 <code>1/(B 的列和+eps)</code>。因此:</p>
<ul>
  <li><code>DictProb = T1*in(:)</code>:对每个字典条目,统计与它相连的前景像素比例(加权平均)。</li>
  <li><code>p = T2*DictProb</code>:对每个像素,把它邻域涉及的字典条目的概率再做一次加权平均。</li>
</ul>
<p>结果就是:<b>纹理相同</b>的像素共享字典词,标签会互相“渗透”;<b>纹理不同</b>的区域字典词不同,概率保持在 0/1 附近。这就是不依赖灰度、而依赖纹理来分割的关键。</p>
<div class="callout info"><b>一句话记忆</b>先把掩膜投影到字典空间(T1),再投影回像素空间(T2);两次投影都是按连接数归一化的平均,所以叫“扩散”。</div>
</div>
</section>

<section id="loop">
<h2>五、水平集演化循环</h2>
<div class="equation">φ ← φ + Δt · δ(φ) · [ log<sub>10</sub>( N(x; m<sub>in</sub>, C<sub>in</sub>) / N(x; m<sub>out</sub>, C<sub>out</sub>) ) + β·κ(φ) − γ ]</div>
<p style="text-align:center;color:#64748b;font-size:13.5px;margin-top:-4px">代码中的 <code>delta</code> 取 <code>(1/ε)·pr<sub>1</sub>·pr<sub>2</sub></code>;<code>log_tr</code> 是区域项,<code>curvature(phi)</code> 是长度项,<code>−gamma</code> 是面积项。</p>

<div class="diagram">
<svg viewBox="0 0 1120 355" role="img" aria-label="水平集迭代循环图">
  <defs>
    <marker id="arrow2" viewBox="0 0 10 10" refX="9" refY="5" markerWidth="7" markerHeight="7" orient="auto-start-reverse">
      <path d="M0,0 L10,5 L0,10 z" fill="#64748b"/>
    </marker>
  </defs>
  <style>
    .lnode{fill:#fff;stroke:#94a3b8;stroke-width:1.5;rx:12}
    .lnode.a{fill:#dbeafe;stroke:#3b82f6}
    .lnode.b{fill:#ede9fe;stroke:#8b5cf6}
    .lnode.c{fill:#dcfce7;stroke:#22c55e}
    .lnode.d{fill:#fef3c7;stroke:#f59e0b}
    .lt{font-size:13.5px;font-weight:700;fill:#0f172a}
    .ls{font-size:11.5px;fill:#64748b}
    .le{stroke:#64748b;stroke-width:1.6;fill:none}
    .ll{font-size:11px;fill:#64748b}
  </style>

  <rect class="lnode a" x="40" y="40" width="220" height="74"/>
  <text class="lt" x="150" y="65" text-anchor="middle">① in = phi &gt; 0</text>
  <text class="ls" x="150" y="85" text-anchor="middle">当前分割掩膜 r×c</text>
  <text class="ls" x="150" y="101" text-anchor="middle">逻辑值,内部为 true</text>

  <rect class="lnode b" x="390" y="40" width="220" height="74"/>
  <text class="lt" x="500" y="65" text-anchor="middle">② DictProb = T1*in(:)</text>
  <text class="ls" x="500" y="85" text-anchor="middle">掩膜 → 字典空间</text>
  <text class="ls" x="500" y="101" text-anchor="middle">(M²K)×1</text>

  <rect class="lnode b" x="740" y="40" width="220" height="74"/>
  <text class="lt" x="850" y="65" text-anchor="middle">③ p = T2*DictProb</text>
  <text class="ls" x="850" y="85" text-anchor="middle">字典空间 → 像素概率</text>
  <text class="ls" x="850" y="101" text-anchor="middle">N×1,值域约 [0,1]</text>

  <rect class="lnode c" x="740" y="200" width="220" height="74"/>
  <text class="lt" x="850" y="225" text-anchor="middle">④ 加权高斯拟合</text>
  <text class="ls" x="850" y="245" text-anchor="middle">pr = [p; 1−p]</text>
  <text class="ls" x="850" y="261" text-anchor="middle">m (d×2), covar (d×d×2)</text>

  <rect class="lnode c" x="390" y="200" width="220" height="74"/>
  <text class="lt" x="500" y="225" text-anchor="middle">⑤ log_tr = log10(N₁/N₂)</text>
  <text class="ls" x="500" y="245" text-anchor="middle">逐像素对数似然比</text>
  <text class="ls" x="500" y="261" text-anchor="middle">r×c</text>

  <rect class="lnode d" x="40" y="200" width="220" height="74"/>
  <text class="lt" x="150" y="225" text-anchor="middle">⑥ φ 更新</text>
  <text class="ls" x="150" y="245" text-anchor="middle">δ = (1/ε)·pr₁·pr₂</text>
  <text class="ls" x="150" y="261" text-anchor="middle">φ += Δt·δ·(log_tr+βκ−γ)</text>

  <!-- edges -->
  <path class="le" d="M260 77 H390" marker-end="url(#arrow2)"/>
  <path class="le" d="M610 77 H740" marker-end="url(#arrow2)"/>
  <path class="le" d="M850 114 V200" marker-end="url(#arrow2)"/>
  <path class="le" d="M740 237 H610" marker-end="url(#arrow2)"/>
  <path class="le" d="M390 237 H260" marker-end="url(#arrow2)"/>
  <path class="le" d="M150 200 V114" marker-end="url(#arrow2)"/>
  <text class="ll" x="160" y="160">下一轮</text>

  <!-- phase labels -->
  <text class="ll" x="300" y="28">纹理扩散(标签 → 概率)</text>
  <text class="ll" x="760" y="188">E 步:给定分割,估计概率/高斯</text>
  <text class="ll" x="300" y="188">M 步:给定概率,更新 φ</text>
  <text class="ls" x="560" y="330" text-anchor="middle">δ 在 p≈0.5 的边界带最大;κ = curvature(φ) 提供平滑;γ 控制收缩/扩张。</text>
</svg>
<div class="caption">图 2 每次迭代的计算链。代码顺序是 ①→②→③→④→⑤→⑥,然后回到 ①。</div>
</div>

<h3>每一项的代码与含义</h3>
<div class="card tight">
<table>
<thead><tr><th style="width:24%">代码</th><th style="width:15%">对应项</th><th>作用</th></tr></thead>
<tbody>
<tr><td><code>in = phi&gt;0</code></td><td>当前分割</td><td>φ 的正区域视为目标内部;得到逻辑掩膜</td></tr>
<tr><td><code>DictProb = T1*in(:)</code></td><td>纹理投影 1</td><td>把像素标签投影到字典空间(按连接数平均)</td></tr>
<tr><td><code>p = T2*DictProb</code></td><td>纹理投影 2</td><td>把字典概率投回每个像素,得到“属于内部”的概率</td></tr>
<tr><td><code>pr(1,:)=p; pr(2,:)=1-p</code></td><td>软权重</td><td>作为两个高斯拟合的样本权重</td></tr>
<tr><td><code>m(:,j) = X*pr(j,:)'/sum(pr(j,:))</code></td><td>加权均值</td><td>前景/背景的平均颜色</td></tr>
<tr><td><code>newCovar = …/sum(pr(j,:))</code></td><td>加权协方差</td><td>颜色的散布程度;加 <code>eps*eye(d)</code> 防止奇异</td></tr>
<tr><td><code>rcond(newCovar) &gt; 1e3*eps</code></td><td>数值保护</td><td>太病态就沿用上一轮的协方差</td></tr>
<tr><td><code>log_t = log10(gaussianValue(...)./gaussianValue(...))</code></td><td>区域项</td><td>每个像素更像前景还是背景;&gt;0 倾向内部</td></tr>
<tr><td><code>delta = (1/epsilon)*pr(1,:).*pr(2,:)</code></td><td>Dirac δ 近似</td><td>只在边界(p≈0.5)附近取值大,内部/外部深处接近 0</td></tr>
<tr><td><code>beta*curvature(phi)</code></td><td>长度/平滑项</td><td>抑制尖角与毛刺,保持曲线平滑</td></tr>
<tr><td><code>-gamma</code></td><td>面积/气球项</td><td><code>gamma=0</code> 时不起作用;正负控制收缩或扩张</td></tr>
<tr><td><code>phi = phi + timestep*delta.*(...)</code></td><td>演化</td><td>更新水平集;<code>timestep</code> 控制步长</td></tr>
</tbody>
</table>
</div>

<div class="callout"><b>“EM 式”理解</b>可以把每轮看成:E 步用当前 φ 估计每个像素的前景/背景软权重,并用这些权重拟合两个高斯;M 步用对数似然比和曲率更新 φ。注意代码里并没有严格的概率归一化,<code>T1/T2</code> 是加权平均算子,所以 <code>p</code> 只是“扩散后的软标签”,不是严格的概率。</div>
</section>

<section id="functions">
<h2>六、函数逐个详解</h2>
<p>每个条目包含:签名、输入/输出、调用关系、作用说明,以及可展开的源码。源码来自工程文件本身;注释保留原样(部分文件是 GBK 编码,已在生成时正确解码)。</p>

<article class="fn">
<h3><span class="badge b-script">脚本</span> demo.m <span class="path">code/demo.m</span></h3>
<div class="sig">脚本,无函数签名;在 code/ 下直接运行 demo</div>
<div class="kv">
  <div>作用</div><div>完整演示“纹理字典 + 水平集”分割:读取测试图 → 建字典 → 建扩散算子 → 手动圈选初始区域 → 迭代 100 次 → 显示结果。</div>
  <div>初始化</div><div><code>binary_image = roipoly;</code> 需要交互式 MATLAB 窗口;<code>phi = 2*binary_image - 1</code>,再 <code>sign(phi)*2</code> 变成 ±2。</div>
  <div>调用</div><div><code>load_test_case</code>、<code>convertJxN</code>、<code>build_km_tree</code>、<code>search_km_tree</code>、<code>biadjacency_matrix</code>、<code>transition_matrix</code>、<code>gaussianValue</code>、<code>curvature</code>、<code>show_contour</code>、<code>showContour</code>、<code>showLevleSetFunction</code></div>
  <div>关键行</div><div>演化主循环在第 84–123 行;核心更新在第 114 行。</div>
</div>
<p>脚本里保留了一些未使用的变量和注释掉的自动初始化方案(第 53–68 行)。<code>demo_v2.m</code>/<code>demo1.m</code> 把这段自动初始化启用,因此可以非交互运行;<code>demo3.m</code> 是 <code>test_case=13</code> 的 21 步变体。</p>
<details><summary>查看 demo.m 完整源码</summary><pre>
clear, close all
addpath texture_functions level_set_functions auxiliary_functions

%% SETTINGS
% loading image and initialization settings
test_case = 1; % chose one of 12 test cases, see with show_test_cases
[input_image,center,radius,patch_size,epsilon,gamma,beta,nr_step] = load_test_case(test_case);
% M used here is 3, 5, 9 or 15: larger works better for textured images

% dictionary settings
nToClust = 5000; % nr random patches for building the dictionary tree
branching_factor = 5; % branching factor for dictionary tree
number_layers = 4; % nr. layers in dictionary tree
    % branching_factor and number_layers relate to dictionary size
normalize = true; % normalization makes texture invariant to global intensity changes

% curve evolution settings
sigma = 1; % sigma for regularization using gaussian smoothing
 % nr. evolution steps 
timestep=1;
initial=0;
%% PRE-PROCESSING
% texture
im = imread(input_image);
im_double = double(im);
[r,c,l] = size(im_double);
if l==1 % just for equal visualization of rgb and gray test cases
    X = double(im(:)'); 
    im = repmat(im,[1,1,3]);    
else
    X=convertJxN(im_double);
end
% building dictionary
tree = build_km_tree(im_double,patch_size,branching_factor,...
    nToClust,number_layers,normalize);
% dictionary assignment
A = search_km_tree(im_double,tree,branching_factor,normalize);
% texture representation
B=biadjacency_matrix(A,patch_size);
[T1,T2] = transition_matrix(B);

means = mean(X, 2)'; % Compute manually now.
vars = std(X'); % Compute manually now.
X = (X - means'*ones(1,size(X,2))) ./ (vars'*ones(1,size(X, 2)));
if max(abs(std(X') - ones(size(std(X'))))) &gt; 0.05
    disp('WARNING: Not all variates are normalized properly!!');
end
gaussian_filter = fspecial('gaussian', [6*round(sigma)+1,1], sigma);
    % gaussian filter for regularization of the level set 
    curvem = []; covar = [];
% level set initialization
% phi=ones(r,c);
% if initial==0    
%     [xx yy] = meshgrid(1:c,1:r);
%     phi = (sqrt(((xx -center(1)).^2 + (yy - center(2)).^2 )) -radius);
%     phi = sign(phi).*2;    
% else
%   z = zeros(2, r*c);
%   h = ones(7,7);
%   [m prior junk covar] = deterministicKmeans(X(:,1:r*c),2);
%   for i=1 : 2
%      z(i,:)=prior(i)*gaussianValue(X, m(:,i), covar(:,:,i));
%      z(i,:)=reshape(filter2(h, reshape(z(i,:),r,c)),1,r*c);
%   end    
%   phi=reshape(log10(z(1,:)./z(2,:)),r,c);
%   phi = filter2(gaussian_filter,filter2(gaussian_filter',phi));
% end
figure, imagesc(uint8(im_double), [0, 255]),  axis off, axis equal;
hold on;
binary_image = roipoly; % 手动初始化轮廓
phi = 2 * binary_image - ones(r,c);
contour(phi, [0 0], 'blue', 'linewidth', 2);
phi = sign(phi) .* 2;
win=5;
h = ones(win,win);
N=r*c;
d=l;
fig1 = figure; 
show_contour(im,phi), title('initialization')
%% EVOLUTION, VARIANT FROM SCIA 2015
fig2 = figure;
 P=ones(r,c,nr_step);
for i = 1:nr_step   
    in = phi&gt;0; % current segmentation
    DictProb = T1*in(:);
%     alpha = sum(in(:))/sum(~in(:)); % area(in)/area(out)
%     DictProb = DictProb./(DictProb+alpha*(1-DictProb)); % area normalization
    p=T2*DictProb;
%     pr(1,:)=p;
%     pr(2,:)=1-p;
%     p1=T2*DP;    
    pr(1,:)=p;
    pr(2,:)=1-p;
    for j = 1:2
        if(sum(pr(j, :)) &gt; 0) %To catch empty clusters
            m(:, j) = (X * pr(j, :)') / sum(pr(j, :));
        end
        newCovar = inv(sum(pr(j, :))) * ( (X - m(:, j)*ones(1,N))*sparse(1:N, 1:N, pr(j,:), N, N) * ...
            (X - m(:, j)*ones(1,N))' );
        newCovar = newCovar + eps*eye(d);
        if rcond(newCovar) &gt; 1e3 * eps
            covar(:, :, j) = newCovar;
        end
        clear newCovar;
    end
    P(:,:,i) = reshape(p,size(A)); % probabilities    
    log_t=log10(gaussianValue(X, m(:,1), covar(:,:,1))./...
             gaussianValue(X,m(:,2), covar(:,:,2))) ;
    log_tr= reshape(log_t',r,c);
    delta=reshape((1/epsilon)*pr(1,:).*pr(2,:),r,c);
%     delta=(1/(4*epsilon))*((sech(phi/(2*epsilon))).^2);
%     delta=(epsilon/pi)./(epsilon^2.+phi.^2);
    phi=phi+timestep*delta.*(log_tr+beta*curvature(phi)-gamma); 
%     phi = filter2(gaussian_filter,filter2(gaussian_filter',phi));
   figure(2);
   showContour(im,phi,i);
%    saveas(gcf,['F:\微云\博士科研\Project\TextureCurveEvolution--ELECTRONICS LETTERS\结果\170104\','phi',num2str(i),'.jpg']);
%    figure(3);
%    imshow(P(:,:,i),[]);
%    saveas(gcf,['F:\微云\博士科研\Project\TextureCurveEvolution--ELECTRONICS LETTERS\结果\170104\','PA',num2str(i),'.jpg']);
   
end

figure(4);
showLevleSetFunction(phi)</pre></details>
</article>

<article class="fn">
<h3><span class="badge b-helper">函数</span> load_test_case <span class="path">code/load_test_case.m</span></h3>
<div class="sig">[image, center, radius, patch_size, epsilon, gamma, beta, nr_step] = load_test_case(test_case, multiphase)</div>
<div class="kv">
  <div>输入</div><div><code>test_case</code>:0–13 的整数;<code>multiphase</code>:文档中说明“true 返回多条曲线”,但当前实现没有使用它。</div>
  <div>输出</div><div><code>image</code> 图片路径;<code>center</code> 圆心(<b>[x, y] = [列, 行]</b>);<code>radius</code> 半径;其余为算法参数。</div>
  <div>调用</div><div>无(只做 <code>switch</code>)。</div>
  <div>被调用</div><div>所有 demo 脚本。</div>
</div>
<p><b>坐标约定:</b>脚本里用 <code>[xx,yy] = meshgrid(1:c,1:r)</code>,所以 <code>xx</code> 是列坐标、<code>yy</code> 是行坐标;<code>center(1)</code> 是 x(列),<code>center(2)</code> 是 y(行)。例如 case 1 的 <code>center=[150;250]</code> 表示“第 150 列、第 250 行”,不是“第 150 行”。</p>
<div class="callout danger"><b>当前版本的限制</b>当前 <code>load_test_case.m</code> 直接返回多圆情形的 <code>center</code>(k×2)和 <code>radius</code>(k 元向量),但没有像旧版那样截断成第一个圆。于是 <code>demo_v2.m</code>/<code>demo1.m</code> 里的 <code>sqrt(...)-radius</code> 在多圆用例(case 2–12)会因维度不匹配报错;只有 <code>radius</code> 为标量的 case 1 和 case 13 能直接跑。旧版 <code>load_test_case(ZGC-…).m</code> 末尾有 <code>center=center(1,:); radius=radius(1);</code> 的截断逻辑。</p>
<details><summary>查看 load_test_case.m 完整源码</summary><pre>function [image,center,radius,patch_size,epsilon,gamma,beta,nr_step] = load_test_case(test_case,multiphase)
% Loads one of test cases for defomable texture segmentation.  可变形纹理分割
% Returns input image and settings for circular curve initialization. 圆曲线初始化
% If multiphase is true returns multiple curves, otherwise only first. 如果多相为真,则返回多条曲线,否则仅返回第一条
        epsilon=0.1;
        gamma=0;
        beta=1.5;
        nr_step = 100;  %演化迭代次数
switch test_case
    case 0       
        image = 'test_images/2.jpg';
        center = [150,200; 250,350]; %初始圆心,demo未使用(仅 demo_v2 等)
        radius = [50,50];   %半径,demo未使用(仅 demo_v2 等)
        patch_size = 3;   %纹理小块边长 M(奇数)
        epsilon=0.005;  %epsilon 在很多数学表示"误差"或"极小量",在水平集中表示正则化参数
        gamma= 0;
        beta=2;
    case 1        
        image = 'test_images/134052.jpg';
        center = [150; 250];
        radius = 50;
        patch_size = 3;
        epsilon=0.005;
        gamma= 0;
        beta=2;
    case 2
        image = 'test_images/124084.jpg';
        center = [150,200; 180,300; 285,185];
        radius = [50,20,20];
        patch_size = 3;
    case 3
        image = 'test_images/12003.jpg';
        center = [150,200; 175,350; 225,465];
        radius = [50,20,20];
        patch_size = 3;
    case 4
        image = 'test_images/108073.jpg';
        center = [150,200; 250,210];
        radius = [50,20];
        patch_size = 5;
    case 6
        image = 'test_images/164074.jpg';
        center = [150,200; 280,410; 250,25];
        radius = [50,20,20];
        patch_size = 5;
    case 5
        image = 'test_images/105053.jpg';
        center = [170,170; 280,250];
        radius = [20,30];
        patch_size = 5; 
    case 7
        image = 'test_images/test_A1.png';
        center = [250,250; 375,375; 120,120];
        radius = [50,50,100];
        patch_size = 11; 
        epsilon=0.001;
        gamma=0;
        beta=1;
    case 8
        image = 'test_images/test_A2.png';
        center = [250,250; 375,375; 120,120];
        radius = [50,50,100];
        patch_size = 9; 
    case 9
        image = 'test_images/test_B1.png';
        center = [250,250; 550,250];
        radius = [50,50];
        patch_size = 9; 
    case 10
        image = 'test_images/test_B2.png';
        center = [250,250; 550,250];
        radius = [50,50];
        patch_size = 9; 
    case 11
        image = 'test_images/test_C.png';
        center = [250,250; 550,250];
        radius = [50,50];
        patch_size = 9; 
    case 12
        image = 'test_images/randen15.png';
        center = [128,32; 128,128; 128,225; 32,128];
        radius = [30,30,30,30];
        patch_size = 15; 
    case 13
        image = 'test_images/0145.jpg';
        center = [239,150];
        radius = 50;
        patch_size = 5; 
        epsilon=0.2;
        gamma=-1;
        nr_step=18;
end

end</pre></details>
</article>

<article class="fn">
<h3><span class="badge b-helper">函数</span> convertJxN <span class="path">code/convertJxN.m(另有一份在 auxiliary_functions/)</span></h3>
<div class="sig">y = convertJxN(A)</div>
<div class="kv">
  <div>输入</div><div><code>A</code>:<code>M×N×J</code>(或更高维)的图像数组。</div>
  <div>输出</div><div><code>y</code>:<code>J×MN</code>,每列是一个像素的 J 个通道值。</div>
  <div>调用</div><div><code>shiftdim</code>(MATLAB 内置)。</div>
  <div>被调用</div><div>所有 demo 的彩色图分支。</div>
</div>
<p>核心是两次 <code>shiftdim</code>:先把通道维移到最前,再转置成 <code>J×MN</code>。它和根目录的 <code>code/convertJxN.m</code> 内容完全相同;由于运行时当前目录 <code>code/</code> 在搜索路径最前,根目录版本通常会先被找到,但两者行为一致。</p>
<details><summary>查看 convertJxN.m 完整源码</summary><pre>function y = convertJxN(A)
% Convert a MxNxJ matrix to JxMN. The conversion is done columnwise.
% Can handle matrices of higher dimensionality, eg 
% can convert MxNxPxJ to JxMNP.
% Example:
%       convertJxN(imread('lenna.jpg'))
%   will give a 3xMN matrix, w/ each column the rgb intensities of a
%   single pixel.
%
% G.Sfikas 29/5/2006
%
% Update 16 Mar 2007: Handles matrices of any dimensionality too, now.
% (Eg 4-D, 5-D etc).
% Fix 16 Mar 2007: y = shiftdim(t2(:,:)', +1) should be t2(:,:) (no
% transpose). _DOUBTFUL_
t2 = shiftdim(A, max(size(size(A))) - 1);
y = shiftdim(t2(:,:),+1)';
return;</pre></details>
</article>

<article class="fn">
<h3><span class="badge b-mex">MEX</span> build_km_tree <span class="path">code/texture_functions/build_km_tree.cpp</span></h3>
<div class="sig">T = build_km_tree(I, M, b, t, L, n)</div>
<div class="kv">
  <div>输入</div><div><code>I</code>:double 图像(2D 或 3D);<code>M</code>:奇数块边长;<code>b</code>:分叉数;<code>t</code>:随机采样块数;<code>L</code>:期望层数;<code>n</code>:是否对每个块做 L2 归一化(可选,默认 false)。</div>
  <div>输出</div><div><code>T</code>:<code>(d·M²)×K</code> 的字典矩阵。默认 <code>M=3, b=5, L=4, t=5000</code> 时 <code>K=5+25+125+625=780</code>,即 <code>27×780</code>。</div>
  <div>内部调用</div><div><code>randperm</code>、<code>sample_patch</code>、<code>k_means</code>(→ <code>set_values</code>/<code>add_values</code>/<code>get_dist</code>)、<code>get_to</code>、<code>sort</code>。</div>
  <div>被调用</div><div>所有 demo 脚本(每张图一次)。</div>
</div>
<p><b>构建过程:</b></p>
<ol>
  <li>从图像内部随机抽 <code>t</code> 个 <code>M×M</code> 块(<code>t</code> 不会超过内部像素数);<code>n=true</code> 时每个块除以自身的 L2 范数,使纹理对整体明暗不敏感。</li>
  <li>第 1 层把所有块聚成 <code>b</code> 类;第 2 层在每个父类的块内再聚成 <code>b</code> 类;依次往下。每个节点内部做 30 轮 k-means,初始中心取该组排序后的前 <code>b</code> 个块。</li>
  <li>字典大小 <code>K = Σ<sub>i=1..L</sub> b<sup>i</sup></code>;如果这个和超过训练块数 <code>t</code>,代码会减少实际层数。运行时打印 <code>Number of nodes in resulting tree: K in L layers.</code></li>
  <li>叶节点之后的位置在矩阵里保持 <code>-1</code>;<code>search_km_tree</code> 遇到 <code>-1</code> 就知道已经到叶子。</li>
</ol>
<div class="callout info"><b>为什么字典树有效</b>每个节点代表一类“局部纹理模式”。逐层细化相当于对纹理空间做层次量化:上层是粗纹理类别,下层是更细的变体。搜索时只需每层比较 <code>b</code> 个中心,复杂度是 <code>O(L·b·d·M²)</code>,比在 780 个原子中暴力搜索快得多。</div>
<details><summary>查看 C++ 关键源码(gateway + k-means)</summary><pre>// The gateway routine
void mexFunction( int nlhs, mxArray *plhs[],
        int nrhs, const mxArray *prhs[])
{
    // input image (I), patch size (M*M), number of nodes in the tree (n), branching
    // factor (b), and number of training patches (n_train). Outputs the km-tree (tree)
    double *I, *tree; // pointers to image and tree
    int b, M, L, n, ndim, n_train; // variables
    const int *dim; // image dimensinos
    int dtree[2]; // tree dimensions
    bool normalize = false; // decide if vectors of image patches should be normalized to unit length
    
    /*  Check for proper number of arguments. */
    /* NOTE: You do not need an else statement when using
     mexErrMsgTxt within an if statement. It will never
     get to the else statement if mexErrMsgTxt is executed.
     (mexErrMsgTxt breaks you out of the MEX-file.)
     */
    if(nrhs &lt; 5 || nrhs &gt; 6)
        mexErrMsgTxt("Five or six inputs required.");
    if(nlhs != 1)
        mexErrMsgTxt("One output required.");
    
    if ( !mxIsDouble(prhs[0]) )
        mexErrMsgTxt("Image should be double!");
    
    // Create a pointer to the input matrix.
    I = mxGetPr(prhs[0]);
    
    // input passing
    double *Md, *bd, *Ld, *n_train_d;
    Md = mxGetPr(prhs[1]);
    M = (int)Md[0];
    
    bd = mxGetPr(prhs[2]);
    b = (int)bd[0];
    
    // check if number of clusters is smaller than branching factor
    if ( n &lt; b ){
        n = b;
    }
    
    n_train_d = mxGetPr(prhs[3]);
    n_train = (int)n_train_d[0];
    
    // determine number of tree nodes
    Ld = mxGetPr(prhs[4]);
    L = (int)Ld[0]; // layers in tree
    n = 0;
    int n_tmp = 0;
    int max_n = (double)n_train/1.0;
//     int max_n = (double)n_train/5.0;
    for ( int i = 0; i &lt; L; i++ ){
        n_tmp += pow((double)b,(i+1));
        if ( n_tmp &gt; max_n ){
            L = i+1;
            break;
        }
        n = n_tmp;
    }
    printf("Number of nodes in resulting tree: %d in %d layers.\n", n, L);
    
    if ( nrhs == 6 ){
        bool *normalize_d;
        normalize_d = (bool *)mxGetData(prhs[5]);
        normalize = normalize_d[0];
    }
    
    // check input properties
    if ( 1 - (M % 2)  || M &lt; 1)
        mexErrMsgTxt("M must be odd and positive.");
    
    if ( n &lt; 1 )
        mexErrMsgTxt("n must be positive.");
    
    if ( b &lt; 1 )
        mexErrMsgTxt("b must be positive.");
    
    // Get the dimensions of the matrix input.
    ndim = mxGetNumberOfDimensions(prhs[0]);
    if (ndim != 2 &amp;&amp; ndim != 3)
        mexErrMsgTxt("search_km_tree only works for 2-dimensional or 3-dimensional images.");
    
    // image dimensions
    dim = mxGetDimensions(prhs[0]);
    
    
    // image struct
    im_st Im;
    Im.im_data = I;
    Im.rows = dim[0];
    Im.cols = dim[1];
    if ( ndim == 3 )
    {
        Im.layers = dim[2];
    }
    else
    {
        Im.layers = 1;
    }
    Im.n_pix = Im.rows*Im.cols;
    
    dtree[0] = Im.layers*M*M;
    dtree[1] = n;
    
    // Set the output pointer to the output matrix. Array initialized to zero.
    plhs[0] = mxCreateNumericArray(2, dtree, mxDOUBLE_CLASS, mxREAL);
    
    // Create a C pointer to a copy of the output matrix.
    tree = mxGetPr(plhs[0]);
    for (int i = 0; i &lt; dtree[0]*dtree[1]; i++ )
        *(tree + i) = -1;
    
    // tree struct
    tree_st Tree;
    Tree.tree_data = tree;
    Tree.n_dim = dtree[0];
    Tree.n_nodes = dtree[1];
    Tree.branch_fac = b;
    Tree.M = M;
    Tree.Mh = (int)(0.5*((double)M-1.0));
    
    // build the km-tree using C++ subroutine
    build_km_tree ( Im, Tree, n_train, normalize );
    </pre><pre style="margin-top:10px">void k_means( vector&lt;im_patch&gt;&amp; patches, tree_st&amp; tree, int f, int t, int node )
{
    // vectors holding the sum of values in the cluster and a vector containing the change
    vector&lt;double&gt; centSum(tree.branch_fac*tree.n_dim), centChange(tree.branch_fac);
    vector&lt;int&gt; centCount(tree.branch_fac); // vector for counting the number of points in a cluster
    double min_d, d, tmp;//, diff; // variables for clustering
    // variables for cluster id and index of previous cluseter, which will be overwritten by new cluster id
    int id = 0, id_in = patches[f].idx;
    
    if ( t-f &gt; tree.branch_fac ){ // check if there are enough point to carry out the clustering
        // initialize the center positions
        for (  int i = 0; i &lt; tree.branch_fac; i++ ){
            set_values(tree, patches[f+i], node+i);
        }
        
        // run clutering for 20 iterations - only little change happens after 10 iterations
        for ( int n_run = 0; n_run &lt; 30; n_run++){
            
            for ( int i = f; i &lt; t; i++ ){ // go through the patches from f to t
                min_d = get_dist(tree, patches[i], node); // initially set min distance and id to the first
                id = 0;
                for ( int j = 1; j &lt; tree.branch_fac; j++ ){ // run throgh the other points
                    d = get_dist(tree, patches[i], node + j); // get the distance
                    if ( d &lt; min_d ){ // if the this cluster is closer set this as min distance
                        min_d = d;
                        id = j;
                    }
                }
                add_values(centSum, patches[i], id, tree.n_dim); // add the patch to the closest cluster
                centCount[id]++; // count the extra patch
                // update the idx to the child idx - note that all layers start with idx = 0
                patches[i].idx = (id + id_in*tree.branch_fac);
            }
            
            // update the clusters in the tree and calculate the center change (not used for anything)
            id = node*tree.n_dim;
            int id_c = 0;
            
            for ( int i = 0; i &lt; tree.branch_fac; i++ ){ // go through all new clusters
                if ( centCount[i] == 0 ){
                    centCount[i] = 1;
                }
                for ( int j = 0; j &lt; tree.n_dim; j++ ){ // go through cluster pixels
                    tmp = centSum[id_c]/centCount[i];
                    //diff = (tmp - *(tree.tree_data + id)); // for calculating center change
                    //centChange[i] += diff*diff;
                    *(tree.tree_data + id) = tmp;
                    id_c++;
                    id++;
                }
            }
            
            // set counter and sum to zero
            fill(centSum.begin(), centSum.end(), 0);
            fill(centCount.begin(), centCount.end(),0);
            fill(centChange.begin(), centChange.end(), 0);
        }
    }
}
</pre></details>
</article>

<article class="fn">
<h3><span class="badge b-mex">MEX</span> search_km_tree <span class="path">code/texture_functions/search_km_tree.cpp</span></h3>
<div class="sig">A = search_km_tree(I, T, b, n)</div>
<div class="kv">
  <div>输入</div><div><code>I</code>:double 图像;<code>T</code>:<code>build_km_tree</code> 的输出;<code>b</code>:分叉数(必须与建树时一致);<code>n</code>:归一化开关(也必须一致)。</div>
  <div>输出</div><div><code>A</code>:<code>r×c</code> 词号图,值 <code>1…K</code>;边界 <code>(M−1)/2</code> 像素为 0。</div>
  <div>内部调用</div><div><code>sample_patch</code>、<code>search_tree</code>、<code>get_dist</code>。</div>
  <div>被调用</div><div>所有 demo 脚本。</div>
</div>
<p><b>搜索过程:</b>对每个内部像素,取它的 <code>M×M</code> 块(同样可选归一化),从树的第 1 层开始,每层在 <code>b</code> 个候选中心中选距离最小的一个,然后下降到它的子节点;遇到 <code>-1</code> 叶标记就停止,返回最后匹配到的节点编号(+1 变成 MATLAB 的 1-based 词号)。</p>
<details><summary>查看 C++ 关键源码(search_tree / search_image)</summary><pre>int search_tree(im_st&amp; im, tree_st&amp; tree, int&amp; r, int&amp; c, bool&amp; normalize)
{
    int node = 0, node_min = -1, node_min_level, next_node; // variables for searching the tree
    double d_min = 10e100, d, d_min_level; 
    
    vector&lt;double&gt; patch = sample_patch(im, tree.M, c, r, normalize); // get the pixel values in a patch
    while ( node &lt; tree.n_nodes ){ // go through the tree
        if ( *(tree.tree_data + node*tree.n_dim) == -1 ){ // check if node is a leaf-node
            return node_min;
        }
        
        d_min_level = 10e100; // set minimum distance to high value
        for ( int i = 0; i &lt; tree.branch_fac; i++ ){ // go through nodes at level 
            next_node = node + i;
            d = get_dist(patch, tree, next_node);
            
            if ( d &lt; d_min_level ){ // set current node to the minimum distance
                d_min_level = d;
                node_min_level = next_node;
            }
        }
        if ( d_min_level &lt; d_min ){ // set overall minimum distance and minimum node
            d_min = d_min_level;
            node_min = node_min_level;
        }
        node = (node_min_level+1)*tree.branch_fac; // go to the child node
    }
    return node_min;
}

// The tree search function applied to the entire image - border is zero and interior is in 1,...,n
void search_image(im_st&amp; im, tree_st&amp; tree, double *A, bool&amp; normalize)
{
    int idx = tree.Mh*im.rows; // increase with empty rows at border
    for ( int i = tree.Mh; i &lt; im.cols-tree.Mh; i++ ){
        idx += tree.Mh; // first Mh pixels are border
        for ( int j = tree.Mh; j &lt; im.rows-tree.Mh; j++ ){           
            *(A + idx) = search_tree(im, tree, i, j, normalize) + 1; // find assignment
            idx++;
        }
        idx += tree.Mh; // last Mh pixels are border
    }
}</pre></details>
</article>

<article class="fn">
<h3><span class="badge b-mex">MEX</span> biadjacency_matrix <span class="path">code/texture_functions/biadjacency_matrix.cpp</span></h3>
<div class="sig">B = biadjacency_matrix(A, M) (文档还写了 B = biadjacency_matrix(A, M, K),但当前参数检查只接受 2 个输入)</div>
<div class="kv">
  <div>输入</div><div><code>A</code>:<code>search_km_tree</code> 输出的词号图;<code>M</code>:块边长。</div>
  <div>输出</div><div><code>B</code>:<code>N×(M²K)</code> 的稀疏 logical 矩阵;<code>N=r·c</code>。</div>
  <div>内部调用</div><div><code>std::sort</code>、<code>mxCreateSparseLogicalMatrix</code>。</div>
  <div>被调用</div><div>所有 demo 脚本。</div>
</div>
<p>对每个内部中心像素 <code>ic</code>,把它的 <code>M×M</code> 邻域中的每个像素 <code>i</code> 与列 <code>j=(c+dx)+(c+dy)*M+(A[ic]−1)*M²</code> 连一条边。所有值都是 <code>true</code>,所以 B 只记录“谁和谁相连”,权重由后面的归一化决定。</p>
<div class="callout danger"><b>参数检查的小 bug</b>注释和文档写了 3 参数版本 <code>(A,M,K)</code>,但 gateway 里是 <code>if ((nlhs != 1) || (nrhs != 2))</code>,第三个参数 <code>K</code> 永远不会被读取;实际调用请只用两个参数,<code>K</code> 会自动取 <code>max(A(:))</code>。</div>
<details><summary>查看 C++ 关键源码(B 的构造与稀疏填充)</summary><pre>void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
{
    /* Check for proper number of input and output arguments */    
    if ((nlhs != 1) || (nrhs != 2))
        mexErrMsgTxt("Usage: B = BIADJANCENCY_MATRIX(A,M) OR "
                "B = BIADJANCENCY_MATRIX(A,M,K).\n");
        
    /* Read input */
    double *A = mxGetPr(prhs[0]); // assignment image
    int X = mxGetM(prhs[0]); // image size X
	int Y = mxGetN(prhs[0]); // image size Y
    int M = (int)mxGetPr(prhs[1])[0]; // patch size
    int K; // number dict patches
    if (nrhs==3)
        K = (int)mxGetPr(prhs[2])[0]; 
    else{ // assumes number of dict patches is max(A)
        K = 0;
        for (int a=0; a&lt;X*Y; a++)
            if (A[a]&gt;K)
                K = A[a];
    }
    
    /* Compute some useful sizes */
    int c = (M-1)/2; // width of boundary having no assignment 
    int n = X*Y; // number of image pixels
    int m = M*M*K; // number of dict pixels
    int s = (X-M+1)*(Y-M+1)*M*M; // number image-dict links (elements in B)
    
    /* Finding elements of B as row-column indices */
    std::vector&lt;ij&gt; bij;
    bij.reserve(s); 
    int ic,i,j;   
    for (int y=0+c; y&lt;Y-c; y++){ // visiting patches centered around pixels
        for (int x=0+c; x&lt;X-c; x++){
            ic = x+y*X; // central patch pixel
            for (int dy=-c; dy&lt;=c; dy++){ // visiting pixels around central
                for (int dx=-c; dx&lt;=c; dx++){
                    i = (x+dx)+(y+dy)*X;
                    j = (c+dx)+(c+dy)*M+(A[ic]-1)*M*M;
                    bij.push_back(ij(i,j));
                }
            }
        }
    }
    
    /* Sorting elements in bij columnwise */
    std::sort (bij.begin(), bij.end());    
    
    /* Placeholder for output */
    plhs[0] = mxCreateSparseLogicalMatrix(n,m,s); // output mxArray, sparse logical matrix B
    if (plhs[0]==NULL)
        mexErrMsgTxt("Could not allocate enough memory!\n");
    
    /* Access fields of output mxArray via pointers  */
    mwIndex *ir = mxGetIr(plhs[0]); // row index (0 indexed) 
    mwIndex *jc = mxGetJc(plhs[0]); // cumulative number of elements per column 
    mxLogical *pr = mxGetLogicals(plhs[0]); // element values (will be all true)
        
    /* Converting row-column indices into row-cumulative column  */
    int k = 0; // for visiting elements of bij
    jc[0] = 0; // first element of cumulative sum is 0
    for (int bc=0; bc&lt;m; bc++){ // all columns of B        
        jc[bc+1] = jc[bc]; 
        while (k&lt;bij.size() &amp;&amp; bij[k].j==bc){
            jc[bc+1]++;
            ir[k] = bij[k].i;
            pr[k] = true;
            k++;
        }
    }
}</pre></details>
</article>

<article class="fn">
<h3><span class="badge b-helper">函数</span> transition_matrix <span class="path">code/texture_functions/transition_matrix.m</span></h3>
<div class="sig">[T1, T2] = transition_matrix(B) 或 T = transition_matrix(B)</div>
<div class="kv">
  <div>输入</div><div><code>B</code>:双邻接矩阵(稀疏或稠密)。</div>
  <div>输出</div><div><code>T1 = Dyinv*B'</code>(<code>(M²K)×N</code>);<code>T2 = Dxinv*B</code>(<code>N×(M²K)</code>)。一个输出时返回 <code>T = Dxinv*B*Dyinv*B' = T2*T1</code>。</div>
  <div>调用</div><div><code>sparse</code>、<code>sum</code>。</div>
  <div>被调用</div><div>所有 demo 脚本。</div>
</div>
<p><code>Dxinv</code> 用行和的倒数,<code>Dyinv</code> 用列和的倒数;加 <code>eps</code> 避免除零。两个矩阵都是按连接数归一化的平均算子,而不是互逆矩阵。</p>
<details><summary>查看 transition_matrix.m 完整源码</summary><pre>function [T1,T2] = transition_matrix(B)
%TRANSITION_MATRIX   Label to probability transition matrices.
% T = TRANSITION_MATRIX(B) builds transition matrix based on the bipartite
% matrix B, obtained from function B = bipartite_matrix(A,M). Transition
% matrix defines a linear mapping between images which have the same size
% as A, and may be used as a diffusion matrix or as a label image (L) to
% probability image (P) transition matrix, P = T*L.
% [T1,T2] = TRANSITION_MATRIX(B) returns intermediate transitions from a
% label image to a dictionary and from a dictionary to a probability image
% so that T = T2*T1. Depending on the size of the image and the dictionary,
% it might be more efficient to compute P = T2*(T1*L).
%
% Authors: vand@dtu.dk

[rc,nm] = size(B);

Dxinv = sparse(1:rc,1:rc,1./(sum(B,2)+eps),rc,rc);
Dyinv = sparse(1:nm,1:nm,1./(sum(B,1)+eps),nm,nm);

if nargout ==1
    T1 = Dxinv*B*Dyinv*B';
else
    T1 = Dyinv*B';
    T2 = Dxinv*B;
end

end</pre></details>
</article>

<article class="fn">
<h3><span class="badge b-helper">函数</span> curvature <span class="path">code/level_set_functions/curvature.m</span></h3>
<div class="sig">K = curvature(X)</div>
<div class="kv">
  <div>输入</div><div><code>X</code>:水平集函数(<code>r×c</code>)。</div>
  <div>输出</div><div><code>K</code>:逐像素曲率近似(<code>r×c</code>)。</div>
  <div>调用</div><div><code>gradient</code>(MATLAB 内置)。</div>
  <div>被调用</div><div>演化循环的每次迭代。</div>
</div>
<p>实现是“单位法向量的散度”:<code>n = ∇φ/|∇φ|</code>,再算 <code>∇·n = n<sub>xx</sub>+n<sub>yy</sub></code>。在零等值线附近,曲率大的地方(尖角、毛刺)会被平滑项拉平。</p>
<details><summary>查看 curvature.m 完整源码</summary><pre>function K = curvature(X)
[bdx,bdy]=gradient(X);
mag_bg=sqrt(bdx.^2+bdy.^2)+1e-10;
nx=bdx./mag_bg;
ny=bdy./mag_bg;
[nxx,nxy]=gradient(nx);
[nyx,nyy]=gradient(ny);
K=nxx+nyy;
return;</pre></details>
</article>

<article class="fn">
<h3><span class="badge b-helper">函数</span> gaussianValue <span class="path">code/level_set_functions/gaussianValue.m</span></h3>
<div class="sig">res = gaussianValue(X, m, covar)</div>
<div class="kv">
  <div>输入</div><div><code>X</code>:<code>d×N</code>,每列一个像素;<code>m</code>:<code>d×1</code> 均值;<code>covar</code>:<code>d×d</code> 协方差。</div>
  <div>输出</div><div><code>res</code>:<code>N×1</code> 的列向量,每个像素的概率密度。</div>
  <div>调用</div><div><code>inv</code>、<code>det</code>、<code>exp</code>。</div>
  <div>被调用</div><div>演化循环每轮两次(前景/背景),以及自动初始化分支。</div>
</div>
<p>公式是标准多维高斯:<code>(2π)<sup>−d/2</sup> |C|<sup>−1/2</sup> exp(−½ (x−m)<sup>T</sup>(C+εI)<sup>−1</sup>(x−m))</code>。代码对协方差额外加 <code>eps*eye(d)</code>,避免奇异。注意返回的是<b>列向量</b>,所以 demo 里先转置再 <code>reshape(log_t',r,c)</code>。</p>
<details><summary>查看 gaussianValue.m 完整源码</summary><pre>function res = gaussianValue(X, m, covar)
% Computes probability of dataset X (each _column_ is one datum),
% given a gaussian pdf, Normal(m, covar).
% Based on G.Sfikas 28 nov 2008 
N = size(X, 2);
d = size(X, 1);
bigM = m * ones(1, N);
mah = sum((X-bigM)' * (inv(covar)+eps*eye(size(X,1))) .* (X-bigM)', 2);
res = (2*pi)^(-0.5*d) * det(covar)^(-0.5) * exp(-0.5*mah);
return;</pre></details>
</article>
<article class="fn">
<h3><span class="badge b-display">显示</span> show_contour <span class="path">code/level_set_functions/show_contour.m</span></h3>
<div class="sig">show_contour(I, phi, lw)</div>
<div class="kv">
  <div>输入</div><div><code>I</code>:底图;<code>phi</code>:水平集(可以是 <code>r×c×k</code> 多相);<code>lw</code>:可选线宽。</div>
  <div>输出</div><div>无(画图)。</div>
  <div>调用</div><div><code>cla</code>、<code>imagesc</code>、<code>image</code>、<code>contour</code>、<code>colormap</code>。</div>
  <div>被调用</div><div>所有 demo 在初始化后调用一次。</div>
</div>
<p>支持灰度、RGB 和多层 φ;<code>cla</code> 会清空当前坐标轴,所以重复调用不会叠加轮廓。线宽默认是 <code>2-(size(phi,3)&gt;2)</code>(多相时细一点)。</p>
<details><summary>查看 show_contour.m 完整源码</summary><pre>function show_contour(I,phi,lw)
% Visualization of the zero level set
% Input: background image I, level set function phi, optional line width lw
% Author: vand@dtu.dk


if nargin &lt; 3
    lw = 2-(size(phi,3)&gt;2);
end

% should handle all our test cases
cla
if size(I,3)==1
    imagesc(I), colormap gray
elseif size(I,3)==3
    image(I)
else
    imagesc(mean(I,3)), colormap gray
end

axis image, hold on
col = 'gmcbrykw';
for i=1:size(phi,3)
    contour(phi(:,:,i),[0 0],col(mod(i-1,8)+1),'LineWidth',lw)
end

end
</pre></details>
</article>

<article class="fn">
<h3><span class="badge b-display">显示</span> showContour <span class="path">code/auxiliary_functions/showContour.m</span></h3>
<div class="sig">showContour(im, phi, i)</div>
<div class="kv">
  <div>输入</div><div><code>im</code>:底图;<code>phi</code>:当前水平集;<code>i</code>:迭代编号,用于标题。</div>
  <div>输出</div><div>无(画图)。</div>
  <div>调用</div><div><code>imagesc</code>、<code>axis</code>、<code>colormap</code>、<code>hold</code>、<code>contour</code>、<code>title</code>。</div>
  <div>被调用</div><div>所有 demo 的每一轮迭代。</div>
</div>
<p>与 <code>show_contour</code> 的区别:<b>不调用 <code>cla</code></b>,所以每轮的轮廓会层层叠加;而且使用 <code>imagesc(im,[0,1])</code>,传入的 <code>im</code> 是 0–255 的 uint8/double 时,按 CLim 语义大部分像素会饱和,底图可能偏白。只影响显示,不影响分割结果。</p>
<details><summary>查看 showContour.m 完整源码</summary><pre>function showContour(im,phi,i)
 imagesc(im,[0, 1]); 
    axis off; 
    axis equal; 
    colormap(gray);
    hold on;
    contour(phi,[0 0],'b','linewidth',2); 
%       contour(reshape(log10(z(1,:)./z(2,:)),imageSize),[0 0],'b','linewidth',2); 
    str=[ ' iterations',num2str(i)];
    title(str);
%     saveas(gcf,['F:\微云同步盘\博士科研\论文\fast-robust-bayesian-level-set\实验结果\0523\','86016128',num2str(iteration),'.jpg']);
    hold off</pre></details>
</article>

<article class="fn">
<h3><span class="badge b-display">显示</span> showLevleSetFunction <span class="path">code/auxiliary_functions/showLevleSetFunction.m</span></h3>
<div class="sig">showLevleSetFunction(phi)</div>
<div class="kv">
  <div>输入</div><div><code>phi</code>:水平集(<code>r×c</code>)。</div>
  <div>输出</div><div>无(3D 图)。</div>
  <div>调用</div><div><code>mesh</code>、<code>contour</code>、<code>view</code>、<code>set</code>。</div>
  <div>被调用</div><div>所有 demo 的结尾。</div>
</div>
<p>把 <code>phi</code> 旋转/翻转后画成 3D 曲面,红色线是 <code>phi=0</code> 的零等值线,也就是最终分割边界。函数名中的 <code>Levle</code> 是原作者的拼写(不是 <code>Level</code>),调用时必须保持一致。</p>
<details><summary>查看 showLevleSetFunction.m 完整源码</summary><pre>function showLevleSetFunction(phi)
mesh(rot90(flipud(phi)));
hold on;  contour(rot90(flipud(phi)), [0,0], 'r','LineWidth',2);
view([-80 35]);
axis on;
% [nrow, ncol]=imageSize;
% axis([1 ncol 1 nrow -5 5]);
set(gca,'ZTick',[-3:1:3]);
set(gca,'FontSize',14)</pre></details>
</article>

<article class="fn">
<h3><span class="badge b-legacy">函数 legacy</span> initial_mask <span class="path">code/level_set_functions/initial_mask.m</span></h3>
<div class="sig">phi = initial_mask(dim, radius, centers, distance)</div>
<div class="kv">
  <div>输入</div><div><code>dim=[r,c]</code>;<code>radius</code> 半径;<code>centers</code> 可为空、单个圆心的 <code>[x,y]</code> 或 k×2;<code>distance</code> 可选,给定时生成等间距圆心网格。</div>
  <div>输出</div><div><code>phi</code>:逻辑掩膜(<code>r×c</code>),多个圆盘取并集。</div>
  <div>调用</div><div><code>ndgrid</code>。</div>
  <div>被调用</div><div>两个 legacy demo。</div>
</div>
<p>默认圆心是图像中心;每个圆盘 <code>(X−cx)²+(Y−cy)² &lt; radius²</code> 取并集。</p>
<div class="callout"><b>坐标顺序提醒</b>函数内部用 <code>[X,Y]=ndgrid(1:dim(1),1:dim(2))</code>,<code>X</code> 沿第一维(行)、<code>Y</code> 沿第二维(列);而 <code>demo_v2</code> 用 <code>meshgrid(1:c,1:r)</code>,第一维是列。混用两套初始化时,<code>center</code> 的 <code>[x,y]</code> 顺序需要核对,否则圆可能被转置到错误位置。</div>
<details><summary>查看 initial_mask.m 完整源码</summary><pre>function phi = initial_mask(dim,radius,centers,distance)
% Initial circular mask for level set function
% Either centers or distance should be given (or neither)

if nargin&lt;3 || isempty(centers)
    centers = (dim+1)/2;
end

if nargin==4 % multiple equaly distanced circles
    cx = radius:distance:dim(1)-radius;
    cy = radius:distance:dim(2)-radius;    
    cx = cx + (1+dim(1)-cx(1)-cx(end))/2; % centering
    cy = cy + (1+dim(2)-cy(1)-cy(end))/2; 
    [CX,CY] = ndgrid(cx,cy);
    centers = [CX(:),CY(:)];
end

phi = false(dim);
[X,Y] = ndgrid(1:dim(1),1:dim(2));

for s = 1:size(centers,1)
    in = (X-centers(s,1)).^2+(Y-centers(s,2)).^2 &lt; radius^2;
    phi = phi | in;
end</pre></details>
</article>

<article class="fn">
<h3><span class="badge b-legacy">函数 legacy</span> mask2sdf <span class="path">code/level_set_functions/mask2sdf.m</span></h3>
<div class="sig">phi = mask2sdf(mask)</div>
<div class="kv">
  <div>输入</div><div><code>mask</code>:逻辑掩膜,可多层(<code>r×c×k</code>)。</div>
  <div>输出</div><div><code>phi</code>:符号距离场(近似)。</div>
  <div>调用</div><div><code>bwdist</code>。</div>
  <div>被调用</div><div>两个 legacy demo。</div>
</div>
<p>公式是 <code>bwdist(mask) − bwdist(~mask) + mask − 0.5</code>。实测在本机 R2025a 上,<b>掩膜内部为负、外部为正</b>;而当前 demo 的约定是 <code>phi&gt;0</code> 为内部。因此把 <code>initial_mask</code> + <code>mask2sdf</code> 接到当前演化循环时,<code>in=phi&gt;0</code> 会选中外部,符号需要反过来(例如 <code>phi=-mask2sdf(mask)</code>),或者把演化里的 <code>phi&gt;0</code> 改成 <code>phi&lt;0</code>。</p>
<div class="callout danger"><b>已实测</b>内部像素输出约 <code>−1.5</code>,外部约 <code>+1.74</code>。这是 legacy 脚本与当前主流程之间一个容易忽略的符号差异。</div>
<details><summary>查看 mask2sdf.m 完整源码</summary><pre>function phi=mask2sdf(mask)
% Mask to signed distance field

phi = zeros(size(mask));
for p=1:size(mask,3)
    phi(:,:,p) = bwdist(mask(:,:,p))-bwdist(~mask(:,:,p))+mask(:,:,p)-0.5;
end</pre></details>
</article>

<article class="fn">
<h3><span class="badge b-unused">未使用</span> relative_probability <span class="path">code/level_set_functions/relative_probability.m</span></h3>
<div class="sig">P_relative = relative_probability(P)</div>
<div class="kv">
  <div>输入</div><div><code>P</code>:<code>r×c×k</code> 的多相概率图。</div>
  <div>输出</div><div><code>P_relative</code>:同尺寸,<code>P_i/(P_i+max_{j≠i}P_j)</code>。</div>
  <div>调用</div><div><code>zeros</code>、<code>max</code>。</div>
  <div>被调用</div><div>当前没有脚本调用它;是为多相/多区域扩展预留的。</div>
</div>
<details><summary>查看 relative_probability.m 完整源码</summary><pre>function P_relative = relative_probability(P)
% Relative probabilities for multi-phase level sets.
% Argumentation in Dahl and Dahl, "Dictionary Based Image Segmentation", SCIA 2015

P_relative = zeros(size(P));
for i=1:size(P,3)
    P_in = P(:,:,i);
    P_out = max(P(:,:,[1:i-1,i+1:end]),[],3);
    P_relative(:,:,i) = P_in./(P_in+P_out);
end
</pre></details>
</article>

<article class="fn">
<h3><span class="badge b-legacy">显示 legacy</span> show_phi <span class="path">code/level_set_functions/show_phi.m</span></h3>
<div class="sig">show_phi(phi)</div>
<div class="kv">
  <div>输入</div><div><code>phi</code>:水平集。</div>
  <div>输出</div><div>无(2D 热图)。</div>
  <div>调用</div><div><code>imagesc</code>、<code>blue_white_red</code>、<code>colormap</code>、<code>contour</code>。</div>
  <div>被调用</div><div>两个 legacy demo 的结尾。</div>
</div>
<p>用蓝-白-红色图显示 φ 的正负,绿色线是零等值线。当前 4 个 demo 改用 <code>showLevleSetFunction</code> 的 3D 曲面,所以它只出现在旧脚本里。</p>
<details><summary>查看 show_phi.m 完整源码</summary><pre>function show_phi(phi)
% Visualization of the level set function phi
% Author: vand@dtu.dk


imagesc(phi,[min(phi(:)),max(phi(:))]), axis image, hold on
colormap(blue_white_red(0,256,[min(phi(:)),max(phi(:))]))
contour(phi,[0 0],'g','LineWidth',2)

end

</pre></details>
</article>

<article class="fn">
<h3><span class="badge b-display">显示</span> blue_white_red <span class="path">code/auxiliary_functions/blue_white_red.m</span></h3>
<div class="sig">cmap = blue_white_red(c, M, CLim)</div>
<div class="kv">
  <div>输入</div><div><code>c</code>:白色对应的值(默认 0);<code>M</code>:色图长度(默认当前色图长度);<code>CLim</code>:色限(默认当前坐标轴的 <code>CLim</code>)。</div>
  <div>输出</div><div><code>cmap</code>:<code>M×3</code> 色图。</div>
  <div>调用</div><div><code>colormap</code>、<code>get</code>、<code>linspace</code>、<code>interp1</code>。</div>
  <div>被调用</div><div>仅 <code>show_phi.m</code>。</div>
</div>
<details><summary>查看 blue_white_red.m 完整源码</summary><pre>function cmap = blue_white_red(c,M,CLim)
%BLUEWHITERED   Blue, white, and red color map.
%   BLUEWHITERED returns blue-white-red colormap, with colormap limits
%   given by the current colormap. White color corresponds  to the c value,
%   which by default is 0 for negative-to-positive colormaps, and mean of
%   colormap limits othervise. The length of the colormap is M, which by
%   default is the length of the current colormap. Optionaly limits of
%   colormap may also be given.
%   Author: vand@dtu.dk

% number of elements in colormap
if nargin &lt; 2
    M = size(colormap,1);
end

if nargin &lt; 1
    c = 0;
end

% finding limits of current colormap
if nargin &lt; 3
    CLim = get(gca, 'CLim');
end

% if current colormap does not include 0, we use only half of the
% red-white-blue colormap
low = CLim(1);
high = CLim(2);

if c&lt;=low || c&gt;= high
    c = (low + high)/2;
end

X = [low; (low+c)/2; c; (high+c)/2; high];
V = [0,0,0.5;
    0,0.5,1;
    1,1,1;
    1,0,0;
    0.5,0,0];
Xq = linspace(low,high,M);
cmap = interp1(X,V,Xq);

end
</pre></details>
</article>

<article class="fn">
<h3><span class="badge b-helper">函数</span> deterministicKmeans <span class="path">code/deterministicKmeans.m</span></h3>
<div class="sig">[centroids, weights, idx, covars] = deterministicKmeans(X, K)</div>
<div class="kv">
  <div>输入</div><div><code>X</code>:<code>d×N</code>,每列一个样本;<code>K</code>:聚类数。</div>
  <div>输出</div><div><code>centroids</code>(<code>d×K</code>)、<code>weights</code>(1×K,归一化比例)、<code>idx</code>(1×N 归属)、<code>covars</code>(<code>d×d×K</code>)。</div>
  <div>调用</div><div><code>kmeans</code>、<code>cov</code>、<code>rcond</code>、<code>eps</code>。</div>
  <div>被调用</div><div>demo 的 <code>initial≠0</code> 分支;默认 <code>initial=0</code> 不执行。<b>本机没有 Statistics and Machine Learning Toolbox,<code>kmeans</code> 不可用</b>,所以该分支目前跑不了。</div>
</div>
<p>它用确定性的方式初始化 <code>K</code> 个中心(沿每个维度取均值 ± 标准差),再调用 MATLAB 的 <code>kmeans</code>;最后用每簇样本计算协方差。注意第 39 行 <code>centroids(i,:) = ... + allDataMean(d)</code> 用的是 <code>d</code>(维度数)而不是 <code>i</code>,中心初始化可能因此偏移;这条自动初始化路径属于未充分验证的遗留代码。</p>
<details><summary>查看 deterministicKmeans.m 完整源码</summary><pre>function [centroids weights idx covars] = deterministicKmeans(X, K)
% deterministicKmeans(X, K)
% Computes K centroids for data in X (each column is one datum), using
% the k-means algorithm. Uses a deterministic way, depending only on X,
% to initialize the k-means algorithm.
%
% Important notes:
%       The kmeans function in MATLAB2006b assumes that input is
%       a matrix where each _row_ is one datum.
%
% Arguments:
% X             Data. Each data is one column. (dxN)
% K             Number of desired centroids. (scalar)
%
% centroids     Resulting centroids. Each centroid is one column. (dxK)
% weights       Percentages showing how many data 'belong' to each
%                   centroid.
% idx           Segmentation (assignment of a cluster to each datum)
% covars        Covariance matrices for each cluster.
% 
%
% See also:
%       kmeans, multistartKmeans
% G.Sfikas 24 Apr 2007, 31 Jan 2008 (segmentation)
%          28 Apr 2008 (check for condition of output covariances)
%
weights = inv(K) * ones(1, K);
[d N] = size(X);
covars = zeros(d, d, K);
allDataMean = mean(X, 2);
allDataCov = cov(X');
centroids = zeros(d, K);
deviations = sqrt(diag(allDataCov));
if K == 1
    centroids(:, 1) = allDataMean;
else
    for i = 1:d
        centroids(i, :) = -deviations(i):(deviations(i)*2)/(K-1):deviations(i);
        centroids(i, :) = centroids(i, :) + allDataMean(d);
    end
    [idx, centroids] = kmeans(X', K, 'start', centroids', ...
        'EmptyAction', 'singleton', 'Maxiter', 200);
    centroids = centroids';
end

for j = 1:K
    weights(j) = sum(idx == j);
    covars(:,:,j) = cov(X(:,idx == j)') + eps*eye(d);
    if(rcond(covars(:,:,j)) &lt;= 2*eps)
        covars(:,:,j) = eye(d);
    end
end
weights = weights / sum(weights);
return;</pre></details>
</article>

<article class="fn">
<h3><span class="badge b-script">脚本</span> compile_mex_functions <span class="path">code/texture_functions/compile_mex_functions.m</span></h3>
<div class="sig">脚本,无函数签名;在 texture_functions/ 下运行 compile_mex_functions</div>
<div class="kv">
  <div>作用</div><div>依次编译 <code>biadjacency_matrix.cpp</code>、<code>build_km_tree.cpp</code>、<code>search_km_tree.cpp</code>,输出同目录的 <code>.mexw64</code>。</div>
  <div>前提</div><div>已配置 C++ 编译器(<code>mex -setup C++</code>);当前目录必须是 <code>texture_functions/</code>,否则 <code>mex</code> 找不到同目录的 <code>.cpp</code>。</div>
</div>
<details><summary>查看 compile_mex_functions.m 完整源码</summary><pre>% compile mex files
mex -largeArrayDims biadjacency_matrix.cpp
mex build_km_tree.cpp % based on Euclidean distance
mex search_km_tree.cpp % based on Euclidean distance

/* mex -largeArrayDims biadjacency_matrix.cpp */
/* Compile the biadjacency_matrix.cpp file using the mex command with the -largeArrayDims option */
/* This command is used to compile the C++ code into a MATLAB executable function, allowing MATLAB to call the C++ code directly */
/* The -largeArrayDims option is used to support large array dimensions, overcoming the limitations of MATLAB's default array size */

/* mex build_km_tree.cpp % based on Euclidean distance */
/* Compile the build_km_tree.cpp file using the mex command */
/* The comment indicates that this code is based on the Euclidean distance calculation method, providing context for the compilation purpose */
/* The mex command compiles C/C++ code into a form that can be called directly in MATLAB, enhancing MATLAB's functionality with high-performance C/C++ code */

/* mex search_km_tree.cpp % based on Euclidean distance */
/* Compile the search_km_tree.cpp file using the mex command */
/* The comment also indicates that this code is based on the Euclidean distance calculation method, emphasizing the consistency of this code with the previous compilation task */
/* This command aims to make the C++ implementation of the k-means tree search algorithm accessible to MATLAB for efficient data search and processing */
</pre></details>
</article>
</section>

<section id="mex">
<h2>七、MEX / C++ 内部结构</h2>
<p>三个 MEX 都是“MATLAB 入口 + 一组静态 C++ 辅助函数”的结构。C++ 内部没有调用 MATLAB 工具箱,因此不依赖 Statistics Toolbox;k-means 是作者自己实现的。</p>

<div class="card">
<div class="filetree"><span class="k">build_km_tree</span>
├─ <span class="k">randperm</span>            <span class="c">// 部分 Fisher–Yates,随机选 t 个块索引</span>
├─ <span class="k">sample_patch</span>          <span class="c">// 取 M×M×layers 块;可选 L2 归一化</span>
├─ <span class="k">k_means</span>               <span class="c">// 每个节点内做 30 轮 b-means</span>
│  ├─ <span class="k">set_values</span>          <span class="c">// 把块复制到节点中心</span>
│  ├─ <span class="k">add_values</span>          <span class="c">// 累加簇内块,用于求均值</span>
│  └─ <span class="k">get_dist</span>            <span class="c">// 平方欧氏距离</span>
└─ <span class="k">get_to</span>                <span class="c">// 按父节点 idx 找到该组块的结束位置</span>

<span class="k">search_km_tree</span>
├─ <span class="k">sample_patch</span>          <span class="c">// 与建树时相同的归一化</span>
└─ <span class="k">search_tree</span>           <span class="c">// 每层选最近中心,贪心下降到叶</span>
   └─ <span class="k">get_dist</span>            <span class="c">// 平方欧氏距离</span>

<span class="k">biadjacency_matrix</span>
└─ <span class="k">std::sort</span> + 稀疏填充      <span class="c">// 收集 (像素, 列) 对,按列排序后写 CSR 结构</span></div>
</div>

<div class="card tight">
<table>
<thead><tr><th style="width:20%">C++ 函数</th><th style="width:16%">文件</th><th>作用</th><th style="width:18%">调用者</th></tr></thead>
<tbody>
<tr><td><code>mexFunction</code></td><td>build_km_tree.cpp</td><td>参数检查、算 K/L、分配 <code>tree</code>、调用 <code>build_km_tree</code></td><td>MATLAB</td></tr>
<tr><td><code>build_km_tree</code></td><td>build_km_tree.cpp</td><td>采样块、逐层聚类、填充树矩阵</td><td>mexFunction</td></tr>
<tr><td><code>sample_patch</code></td><td>两者都有</td><td>按中心坐标取 <code>M×M×layers</code> 块;<code>normalize</code> 时除以 L2 范数</td><td>build / search</td></tr>
<tr><td><code>randperm</code></td><td>build_km_tree.cpp</td><td>随机选择训练块索引(不重复)</td><td>build_km_tree</td></tr>
<tr><td><code>k_means</code></td><td>build_km_tree.cpp</td><td>对一个父节点下的块做 <code>b</code> 类聚类,30 轮</td><td>build_km_tree</td></tr>
<tr><td><code>set_values</code> / <code>add_values</code></td><td>build_km_tree.cpp</td><td>复制/累加块到中心</td><td>k_means</td></tr>
<tr><td><code>get_dist</code></td><td>两者都有</td><td>平方欧氏距离</td><td>k_means / search_tree</td></tr>
<tr><td><code>get_to</code></td><td>build_km_tree.cpp</td><td>找到某个父节点 idx 对应块组的末尾</td><td>build_km_tree</td></tr>
<tr><td><code>search_tree</code></td><td>search_km_tree.cpp</td><td>单像素贪心搜索:每层在 <code>b</code> 个节点中选最近,直到 <code>-1</code> 叶</td><td>search_image</td></tr>
<tr><td><code>search_image</code></td><td>search_km_tree.cpp</td><td>遍历内部像素,写词号;边界写 0</td><td>mexFunction</td></tr>
</tbody>
</table>
</div>

<div class="callout info"><b>两个必须保持一致的参数</b>建树和搜索时的 <code>b</code>(分叉数)与 <code>normalize</code> 必须一致;否则树的结构或距离尺度对不上,<code>A</code> 会退化成随机词号。 <code>M</code> 不需要显式传给 <code>search_km_tree</code>,它从 <code>tree</code> 的行数反推:<code>M=sqrt(size(T,1)/layers)</code>。</div>
</section>

<section id="variants">
<h2>八、demo 脚本对比</h2>
<div class="card tight">
<table>
<thead><tr><th style="width:19%">脚本</th><th style="width:8%">test_case</th><th style="width:17%">初始化</th><th style="width:10%">迭代</th><th style="width:10%">面积归一化</th><th>额外显示 / 备注</th></tr></thead>
<tbody>
<tr><td><code>demo.m</code></td><td>1</td><td><code>roipoly</code> 手动圈选</td><td>100(<code>nr_step</code>)</td><td>否</td><td>结尾 3D;需要交互式窗口;自动初始化整段被注释</td></tr>
<tr><td><code>demo_v2.m</code> / <code>demo1.m</code></td><td>1</td><td>圆心 + 半径 → <code>sign(...)*2</code></td><td>100</td><td>否</td><td>结尾 3D;两者仅换行/赋值格式不同;可 <code>-batch</code> 运行</td></tr>
<tr><td><code>demo3.m</code></td><td>13</td><td>圆心 + 半径</td><td><b>21(硬编码)</b></td><td>否</td><td>每轮 <code>imshow(P(:,:,i))</code>;结尾画 <code>log(P/(1-P))</code>;case 13 的 <code>nr_step=18</code> 被忽略</td></tr>
<tr><td><code>demo(SSCE-…).m</code></td><td>3</td><td><code>initial_mask</code> + <code>mask2sdf</code></td><td>10(硬编码)</td><td><b>是</b>(<code>alpha</code>)</td><td>结尾 <code>show_phi</code>;legacy;无空区域保护</td></tr>
<tr><td><code>demo(ZGC-…).m</code></td><td>5</td><td><code>initial_mask</code> + <code>mask2sdf</code></td><td>20(硬编码)</td><td><b>是</b>(<code>alpha</code>)</td><td>结尾 <code>show_phi</code>;legacy;有空区域保护</td></tr>
</tbody>
</table>
</div>

<h3>自动初始化的两种写法</h3>
<div class="grid2">
  <div class="card tight">
    <h4>当前版(demo_v2 / demo1 / demo3)</h4>
    <p class="small">用 <code>meshgrid(1:c,1:r)</code>,<code>center</code> 的第一个分量是列坐标 x。只适用于标量半径。</p>
    <pre>if initial==0
    
    [xx yy] = meshgrid(1:c,1:r);
    phi = (sqrt(((xx -center(1)).^2 + (yy - center(2)).^2 )) -radius);
    phi = sign(phi).*2;    </pre>
  </div>
  <div class="card tight">
    <h4>手动版(demo.m)</h4>
    <p class="small">用 <code>roipoly</code> 交互圈选;返回值是逻辑掩膜,再用 <code>2*mask-1</code> 与 <code>sign</code> 变成 ±2。</p>
    <pre>figure, imagesc(uint8(im_double), [0, 255]),  axis off, axis equal;
hold on;
binary_image = roipoly; % 手动初始化轮廓
phi = 2 * binary_image - ones(r,c);
contour(phi, [0 0], 'blue', 'linewidth', 2);
phi = sign(phi) .* 2;</pre>
  </div>
</div>

<h3>演化主循环源码(demo.m 第 84–123 行)</h3>
<details><summary>展开循环源码</summary><pre>for i = 1:nr_step   
    in = phi&gt;0; % current segmentation
    DictProb = T1*in(:);
%     alpha = sum(in(:))/sum(~in(:)); % area(in)/area(out)
%     DictProb = DictProb./(DictProb+alpha*(1-DictProb)); % area normalization
    p=T2*DictProb;
%     pr(1,:)=p;
%     pr(2,:)=1-p;
%     p1=T2*DP;    
    pr(1,:)=p;
    pr(2,:)=1-p;
    for j = 1:2
        if(sum(pr(j, :)) &gt; 0) %To catch empty clusters
            m(:, j) = (X * pr(j, :)') / sum(pr(j, :));
        end
        newCovar = inv(sum(pr(j, :))) * ( (X - m(:, j)*ones(1,N))*sparse(1:N, 1:N, pr(j,:), N, N) * ...
            (X - m(:, j)*ones(1,N))' );
        newCovar = newCovar + eps*eye(d);
        if rcond(newCovar) &gt; 1e3 * eps
            covar(:, :, j) = newCovar;
        end
        clear newCovar;
    end
    P(:,:,i) = reshape(p,size(A)); % probabilities    
    log_t=log10(gaussianValue(X, m(:,1), covar(:,:,1))./...
             gaussianValue(X,m(:,2), covar(:,:,2))) ;
    log_tr= reshape(log_t',r,c);
    delta=reshape((1/epsilon)*pr(1,:).*pr(2,:),r,c);
%     delta=(1/(4*epsilon))*((sech(phi/(2*epsilon))).^2);
%     delta=(epsilon/pi)./(epsilon^2.+phi.^2);
    phi=phi+timestep*delta.*(log_tr+beta*curvature(phi)-gamma); 
%     phi = filter2(gaussian_filter,filter2(gaussian_filter',phi));
   figure(2);
   showContour(im,phi,i);
%    saveas(gcf,['F:\微云\博士科研\Project\TextureCurveEvolution--ELECTRONICS LETTERS\结果\170104\','phi',num2str(i),'.jpg']);
%    figure(3);
%    imshow(P(:,:,i),[]);
%    saveas(gcf,['F:\微云\博士科研\Project\TextureCurveEvolution--ELECTRONICS LETTERS\结果\170104\','PA',num2str(i),'.jpg']);
   
end</pre></details>
</section>

<section id="legacy">
<h2>九、遗留与未使用代码</h2>
<div class="card">
<ul>
  <li><b>完全未调用:</b><code>relative_probability.m</code>(多相扩展预留),当前 5 个脚本都没有引用。</li>
  <li><b>仅 legacy 脚本使用:</b><code>initial_mask.m</code>、<code>mask2sdf.m</code>、<code>show_phi.m</code>、<code>blue_white_red.m</code>。</li>
  <li><b>可选但当前跑不了:</b><code>deterministicKmeans.m</code> 依赖 <code>kmeans</code>,本机未安装 Statistics and Machine Learning Toolbox(<code>exist('kmeans')==0</code>)。默认 <code>initial=0</code> 不会走到这里。</li>
  <li><b>demo.m 里的死代码:</b><code>gaussian_filter</code>(唯一使用处被注释)、<code>curvem</code>、<code>win</code>/<code>h</code>、<code>initial</code>(自动初始化整段被注释)、<code>P</code>(当前 demo.m 不显示)。</li>
  <li><b>重复文件:</b><code>convertJxN.m</code> 在根目录和 <code>auxiliary_functions/</code> 各一份(内容相同);<code>demo1.m</code> 与 <code>demo_v2.m</code> 基本相同;<code>load_test_case.asv</code> 是编辑器备份。</li>
  <li><b>旧自动保存:</b>两个带括号时间戳的 <code>demo</code>/<code>load_test_case</code> 文件是 MATLAB 自动保存的旧实验版本,不属于当前主流程。</li>
  <li><b>IDE 文件:</b><code>texture_functions/.idea/</code> 是 JetBrains 工程配置,可忽略。</li>
</ul>
</div>
</section>

<section id="gotchas">
<h2>十、坑点与注意事项</h2>
<div class="card tight">
<table>
<thead><tr><th style="width:7%">级别</th><th style="width:27%">问题</th><th>说明 / 建议</th></tr></thead>
<tbody>
<tr><td><span class="badge b-unused">高</span></td><td>工作目录必须是 <code>code/</code></td><td>所有图片路径都是 <code>test_images/…</code>。在别的目录运行 <code>demo</code> 会找不到图;DDLCN 也一样要求其根目录。</td></tr>
<tr><td><span class="badge b-unused">高</span></td><td><code>test_case=0</code> 的图片不存在</td><td><code>load_test_case.m</code> 指向 <code>test_images/2.jpg</code>,目录里没有这个文件。可运行的用例是 1–13。</td></tr>
<tr><td><span class="badge b-unused">高</span></td><td>多圆用例 + 自动初始化报错</td><td>已验证:<code>radius</code> 为向量时 <code>sqrt(...)-radius</code> 维度不匹配。当前自动初始化只适合 case 1 / case 13;其他用例用 <code>demo.m</code> 手动圈选,或先把 <code>radius</code> 截断成第一个圆。</td></tr>
<tr><td><span class="badge b-unused">高</span></td><td><code>mask2sdf</code> 符号相反</td><td>已验证:内部为负、外部为正,而主流程约定 <code>phi&gt;0</code> 为内部。接 legacy 初始化时要取反,否则前景/背景颠倒。</td></tr>
<tr><td><span class="badge b-unused">高</span></td><td>GBK 编码文件</td><td><code>demo.m</code>、<code>demo1.m</code>、<code>demo3.m</code>、<code>demo_v2.m</code>、<code>showContour.m</code> 是 GBK/legacy 编码,直接按 UTF-8 编辑会破坏中文注释。</td></tr>
<tr><td><span class="badge b-unused">高</span></td><td><code>kmeans</code> 不可用</td><td>本机没有 Statistics and Machine Learning Toolbox;<code>deterministicKmeans</code> 的自动初始化分支会报 “函数或变量 'kmeans' 无法识别”。主流程不受影响。</td></tr>
<tr><td><span class="badge b-legacy">中</span></td><td><code>showContour</code> 不清屏、色限固定 <code>[0,1]</code></td><td>每轮轮廓叠加;底图可能因 CLim 饱和而偏白。只影响观察,不影响数值结果。</td></tr>
<tr><td><span class="badge b-legacy">中</span></td><td><code>figure(2)</code> 指向初始化图窗</td><td><code>demo.m</code> 先建 Figure 1(roipoly),再建 <code>fig1=figure</code>(Figure 2),再建 <code>fig2=figure</code>(Figure 3);循环里的 <code>figure(2)</code> 实际激活的是初始化图窗,不是 <code>fig2</code>。</td></tr>
<tr><td><span class="badge b-legacy">中</span></td><td><code>covar</code> / <code>m</code> 可能未赋值</td><td>若首轮 <code>rcond</code> 检查一直不通过,<code>covar</code> 仍是 <code>[]</code>;若某区域概率和为 0,<code>m(:,j)</code> 不会更新。多数 demo 有 <code>if sum(pr)&gt;0</code> 保护,但 legacy SSCE 版本没有。</td></tr>
<tr><td><span class="badge b-legacy">中</span></td><td><code>deterministicKmeans</code> 的 <code>allDataMean(d)</code></td><td>第 39 行用维度数 <code>d</code> 作为索引,疑似应为 <code>allDataMean(i)</code>;该路径未充分验证。</td></tr>
<tr><td><span class="badge b-legacy">中</span></td><td><code>biadjacency_matrix</code> 的 3 参数版本不可用</td><td>gateway 的 <code>nrhs != 2</code> 检查会让第三个参数直接报错;请只传 <code>(A,M)</code>。</td></tr>
<tr><td><span class="badge b-legacy">中</span></td><td><code>initial_mask</code> 与 <code>demo_v2</code> 的坐标顺序不同</td><td><code>initial_mask</code> 的 <code>X</code> 沿行、<code>Y</code> 沿列;<code>demo_v2</code> 的 <code>xx</code> 沿列。两套初始化混用时核对 <code>center=[x,y]</code> 的顺序。</td></tr>
<tr><td><span class="badge b-helper">低</span></td><td>MEX 版本兼容</td><td>当前 <code>.mexw64</code> 是 2016 年编译的,已在本机 R2025a 验证可用;换机器/版本若报 “Invalid MEX-file”,运行 <code>compile_mex_functions</code> 重编。</td></tr>
<tr><td><span class="badge b-helper">低</span></td><td><code>gaussianValue</code> 的数值稳定性</td><td>用 <code>inv</code> 和 <code>det</code> 直接算高斯;协方差加 <code>eps*eye(d)</code> 并做 <code>rcond</code> 检查,但极端情况下仍可能下溢/病态。</td></tr>
<tr><td><span class="badge b-helper">低</span></td><td>内存与速度</td><td><code>B</code> 是 <code>N×(M²K)</code> 稀疏矩阵;<code>M=15</code> 时列数很大,但每行非零元固定为 <code>M²</code>。大图建议先用较小 <code>patch_size</code> 或减少 <code>nToClust</code> 试跑。</td></tr>
</tbody>
</table>
</div>
</section>

<section id="run">
<h2>十一、运行速查</h2>
<div class="card">
<h3>交互式运行(推荐)</h3>
<pre><span style="color:#94a3b8"># PowerShell:启动 MATLAB(不在 PATH,需用完整路径)</span>
&amp; 'E:\xbb\ProgramFiles\MATLAB\R2025a\bin\matlab.exe'</pre>
<pre><span style="color:#94a3b8">% MATLAB 命令窗口</span>
cd('E:/xbb/Workspace/MATLAB/ZiDianXueXi/code')
addpath texture_functions level_set_functions auxiliary_functions
demo        <span style="color:#64748b">% 手动 roipoly 圈选初始区域(需要交互)</span>
demo_v2     <span style="color:#64748b">% 自动圆初始化;等价于 demo1</span>
demo3       <span style="color:#64748b">% test_case=13,21 步,显示概率图</span></pre>

<h3>无界面运行</h3>
<pre>&amp; 'E:\xbb\ProgramFiles\MATLAB\R2025a\bin\matlab.exe' -batch "cd('E:/xbb/Workspace/MATLAB/ZiDianXueXi/code'); demo_v2"</pre>
<p class="small muted"><code>demo.m</code> 不能用 <code>-batch</code>,因为 <code>roipoly</code> 需要鼠标交互。</p>

<h3>重编 MEX(换 MATLAB 版本或报 Invalid MEX-file 时)</h3>
<pre>cd('E:/xbb/Workspace/MATLAB/ZiDianXueXi/code/texture_functions')
compile_mex_functions</pre>

<h3>换测试用例</h3>
<p>打开 <code>demo.m</code> / <code>demo_v2.m</code> / <code>demo1.m</code> / <code>demo3.m</code>,修改开头的 <code>test_case</code>;参数表在 <code>load_test_case.m</code>。注意:自动初始化建议先用 case 1(<code>134052.jpg</code>)或 case 13(<code>0145.jpg</code>),多圆用例改用 <code>demo.m</code> 手动圈选。</p>

<h3>验证环境是否就绪</h3>
<pre><span style="color:#94a3b8">% 三个 MEX 应返回 3(MEX 文件)并能实际调用</span>
exist('build_km_tree','file')
exist('search_km_tree','file')
exist('biadjacency_matrix','file')
<span style="color:#94a3b8">% 图像处理工具箱函数(roipoly/bwdist)应存在</span>
exist('roipoly','file')</pre>
</div>

<div class="callout ok"><b>建议的阅读顺序</b>先看 <code>demo.m</code> 的“设置 → 预处理 → 初始化 → 循环”四段,再看 <code>transition_matrix.m</code> 和 <code>gaussianValue.m</code>,最后对照 <code>demo.m逐行分析.md</code> 与本文的图 1/图 2。C++ 只在需要理解字典树或改 MEX 时展开。</div>
</section>

<h2 style="margin-top:60px">附:一句话总结每个函数</h2>
<div class="card tight">
<table>
<thead><tr><th style="width:26%">函数</th><th>一句话</th></tr></thead>
<tbody>
<tr><td><code>load_test_case</code></td><td>按编号返回图片路径和全部算法参数。</td></tr>
<tr><td><code>convertJxN</code></td><td>把彩色图变成“每列一个像素”的 <code>d×N</code> 矩阵。</td></tr>
<tr><td><code>build_km_tree</code></td><td>用随机小块逐层 k-means,建一棵纹理字典树。</td></tr>
<tr><td><code>search_km_tree</code></td><td>每个像素沿树贪心下降,得到它的字典词号。</td></tr>
<tr><td><code>biadjacency_matrix</code></td><td>把“像素 ↔ 词 + 块内位置”写成稀疏双邻接矩阵。</td></tr>
<tr><td><code>transition_matrix</code></td><td>由 B 得到标签→字典、字典→像素的两个归一化算子。</td></tr>
<tr><td><code>curvature</code></td><td>算水平集曲率,提供边界平滑力。</td></tr>
<tr><td><code>gaussianValue</code></td><td>算每个像素在给定高斯下的概率密度。</td></tr>
<tr><td><code>show_contour</code></td><td>初始化时显示底图和初始零等值线(会清屏)。</td></tr>
<tr><td><code>showContour</code></td><td>每轮显示当前零等值线(不清屏,会叠加)。</td></tr>
<tr><td><code>showLevleSetFunction</code></td><td>把最终 φ 画成 3D 曲面并标出零等值线。</td></tr>
<tr><td><code>initial_mask</code></td><td>生成一个或多个圆盘的并集掩膜(legacy)。</td></tr>
<tr><td><code>mask2sdf</code></td><td>把掩膜转成符号距离场(注意符号与主流程相反)。</td></tr>
<tr><td><code>relative_probability</code></td><td>多相水平集的相对概率(当前未使用)。</td></tr>
<tr><td><code>show_phi</code></td><td>φ 的 2D 蓝白红热图(legacy)。</td></tr>
<tr><td><code>blue_white_red</code></td><td>生成蓝-白-红色图。</td></tr>
<tr><td><code>deterministicKmeans</code></td><td>确定性 k-means + 高斯拟合,供自动初始化使用(依赖 kmeans 工具箱)。</td></tr>
<tr><td><code>compile_mex_functions</code></td><td>一键编译三个 C++ MEX。</td></tr>
</tbody>
</table>
</div>

<p class="muted small" style="margin-top:34px">本文档由源码静态分析生成;图 1、图 2 为内联 SVG,不依赖任何外部 CDN 或 JavaScript。代码块来自工程原文件,部分 GBK 文件已在生成时按 GBK 解码。行号以当前文件为准。</p>
</main>
</div>
<a class="top" href="#overview">↑ 顶部</a>
</body>
</html>

posted @ 2026-10-20 22:21  小小冰&碎冰冰  阅读(5)  评论(0)    收藏  举报