企业官网建设流程全解析

简介：这份zip压缩包是密歇根图像重建工具箱MIRT的Matlab版本，由美国密歇根大学团队开源，面向医学成像、遥感、半导体检测等领域的图像重建研究者和工程师，可解决CT、MRI、PET等模态下的重建与仿真问题。资源共1569个文件、约2.07MB，其中m脚本1187个、c源文件87个、h头文件50个，并辅以readme说明、tex文档、mat数据等，涵盖滤波反投影（FBP）、代数重建（ART/SART）及期望最大化（EM）等可调用算法，便于用户从源码层面理解实现细节并扩展自定义功能。目前已有219人学习浏览，适合具备Matlab基础、希望快速上手MIRT的研究者。包内除主程序外，还收录元素周期表各元素的质量衰减数据、不同组织/材料参数、测试脚本和平台编译文件，可直接配合MIRT的数据模拟模块构建数字体模、模拟不同成像系统的探测器响应、噪声特性及射线传递过程，验证重建效果，也可用于科研教学中的图像质量对比与算法参数优化。目录结构完整，便于按需检索和二次开发。

1. 图像重建绕不开的 MIRT：这个 Matlab 工具箱到底解决了什么

如果你手里有一组 CT 或 PET 投影数据，用 Matlab 自带的 iradon 直接反投影，出来往往糊成一片，边缘全是放射状伪影——不是 Matlab 不行，而是你没把成像物理和噪声统计放进重建过程。Michigan Image Reconstruction Toolbox (MIRT) 就是干这个的：它是密歇根大学 Fessler 组维护的 Matlab 图像重建工具箱，把从滤波反投影（FBP）到统计迭代重建（OS-SQS、CG、ADMM 类算法）这一整套方法，打包成了能直接调用的类与函数。你拿到的这个 Matlab 版压缩包，解压后就是完整源码目录，不是编译好的黑匣子。它适合三类人：被重建论文公式绕晕、想复现结果的人；手头有投影数据但不知道怎么建模系统矩阵的人；被 Matlab 图像处理大作业逼到墙角、想走正规路径的人。下面拆解基于我实际跑过的流程，照着走，一晚上能跑通。

2. 拆包与安装：把 MIRT 源码装进你的 Matlab 环境

2.1 目录结构：先分清哪些目录进路径、哪些目录先别碰

拿到压缩包先别急着全盘解压到桌面，MIRT 的目录结构是有讲究的。解压后你会看到mirt、data、doc、examples、mex这几个主要目录。mirt是算法主体，按ct、mri、pet、reg子模块分门别类，重头戏都在这里；data里是测试数据和生成测试数据的脚本，体积不大但非常管用；doc里有 Fessler 那套图像重建讲义的 PDF，我最常翻的是正则化那一部分，公式和代码能对上；examples是官方 demo 脚本，安装验证全靠它。

至于mex目录，第一次接触时先别碰。MIRT 较早版本里有一部分系统矩阵算子是用 C 写的 mex，需要本地编译；但近几年版本里，绝大多数算子已经改成纯 Matlab 实现了。你先把主路径配好，如果运行时报缺少某个Gtomo开头的 mex 函数，再回头编译它，省得一开始就卡在编译器配置上。另一个建议是解压路径保持纯英文，不要带空格。Matlab 的类加载对路径字符非常敏感，我见过好几次image_geom类加载失败，最后发现是路径里夹了个中文文件夹名。

2.2 三条命令完成安装与路径固化

安装 MIRT 不需要安装程序，本质就是把源码目录加进 Matlab 搜索路径。打开 Matlab，在命令窗口执行：

% 递归添加 MIRT 所有子目录到 Matlab 搜索路径 addpath(genpath('D:/tools/mirt')); % 保存路径配置，下次启动 Matlab 不用重复添加 savepath

第一行是你的实际解压路径，Windows 下正斜杠反斜杠都能识别。genpath的作用是把这个根目录下所有嵌套子目录都递归加进去，这样examples里的 demo 脚本引用mirt子模块时不会断。第二行savepath把路径配置写进 Matlab 的 pathdef 文件，下次启动自动生效。如果savepath报错说没有写权限，常见原因是 Matlab 装在了C:\Program Files下，启动时用管理员权限运行一次就行了。

我不建议把整个mirt解开以后一股脑扔进某个已有的大型工具箱目录里。MIRT 的函数虽然短名多，但前缀基本都很独特，一般不会和其他库冲突。不过如果你机器上同时装了 ASTRA 或 TIGRE 这类重建库，还是把addpath(genpath(...))改成只加mirt和examples两个目录更稳妥，避免同名函数互相覆盖。

2.3 验证装没装对：which 探测与官方 demo

路径加完后先做一次“探测”，确认关键函数真的被 Matlab 找到了：

% 验证三个核心入口是否可见 which fbp2 which os_sqs which poisson

如果每一行都返回了完整路径，说明路径配置成功；如果返回“未找到”，那就是genpath没包含到对应子目录，或者路径没有保存。接着跑官方 CT 重建示例，这一步是最快的端到端验证：

% 完整跑一遍官方 CT demo：模拟投影、FBP、迭代重建 demo_ct

demo_ct会依次做几件事：生成 Shepp-Logan 椭圆幻影、模拟 CT 扫描投影、加泊松噪声、做 FBP 重建和迭代重建，最后弹出几张图像窗口。我第一次跑的时候，图像弹出来的一瞬间就被迭代重建的效果震住了：同样的输入数据，FBP 图上布满星芒状伪影，迭代重建图干净得多。如果你跑完不报错、图像能正常弹出，环境就闭环了。

关于版本兼容性，我测试过的 R2021b 和 R2023b 都能完整跑通 MIRT 主流程。差异主要在两处：一是正则化对象从Robject这种传统结构体风格迁移到了Reg1类风格；二是部分系统矩阵构建函数改了名字或增加了类封装。如果你还在用 R2020a 之前的老版本，建议找 MIRT 历史 release 标签里的旧版本源码，直接用最新源码反而会遇到接口不匹配的问题。

3. 跑通第一个 CT 重建：从椭圆幻影到 FBP 与 OS-SQS

3.1 为什么用椭圆幻影而不是真实临床数据

MIRT 的官方 demo 几乎全部基于模拟数据，原因很简单：只有当你手里有“真值”时，重建算法的误差才是一个可以计算的数。临床 CT 数据没有标准答案，你怎么知道某条伪影是算法问题还是物理扫描问题？椭圆幻影（Shepp-Logan Phantom）是一组数学定义的椭圆叠加，投影可以精确计算，重建结果可以直接和真值做 RMSE、PSNR，这种定量能力是调参的基础。

模拟数据还有个好处：你可以精确控制噪声水平。真实 CT 的泊松噪声强度取决于 X 射线光子数，而模拟数据里光子数是自己设的，想测算法鲁棒性就调低光子数，想测收敛行为就调高光子数。这个“可控性”在初学阶段比什么都重要。我建议新手不要急着找真实数据，先用默认模拟流程把 MIRT 的几何定义、系统矩阵、重建算法三个环节全部打通，再换真实数据不迟。

3.2 先定义几何：image_geom 与 ct_geom 的参数在说什么

MIRT 的核心抽象是“几何对象”，它把重建网格和扫描几何分开描述。重建网格用image_geom，CT 扫描几何用ct_geom，先定义几何再生成系统矩阵，这是 MIRT 整个 Matlab 架构里最值得学的一层设计。看代码：

% 重建网格：256x256，像素间距 1mm ig = image_geom('nx', 256, 'ny', 256, 'dx', 1, 'dy', 1); % CT 扫描几何：扇束扫描，探测器方向采样 128，径向采样 256，120 个视角，全周 360 度 cg = ct_geom('nf', 128, 'nb', 256, 'na', 120, 'orbit', 360); % 生成 Shepp-Logan 椭圆幻影，作为重建真值 xtrue = ellipse_im(ig, 'shepp-logan'); % 显示真值图像 im(ig, xtrue);

image_geom的nx、ny是重建矩阵的宽和高，dx、dy是每个像素的物理尺寸（毫米）。这四个参数直接决定重建分辨率：nx越大图像越细，但系统矩阵的规模也越大，迭代重建耗时成倍上涨。ct_geom里nf和nb联合描述探测器方向的采样密度，na是投影视角数，orbit是扫描角度范围，CT 通常设 360 或 180。视角数na越密，数据量越大，重建角度采样越充分；欠采样时图像会出现角度方向的重影伪影。

ellipse_im第三个参数传'shepp-logan'是 MIRT 预置的经典幻影，不传也能生成默认椭圆，但 Shepp-Logan 的灰度层次更接近真实人体组织分布，对比算法效果更明显。im(ig, xtrue)是 MIRT 自带的显示函数，它会自动读ig的尺寸信息正确显示图像，比直接imagesc省事得多。

3.3 系统矩阵是隐式算子：正投影与噪声模拟

几何定义好之后，下一步是生成系统矩阵 A。这个环节是 MIRT 与普通 Matlab 图像处理代码最大的分水岭：

% 生成系统矩阵（积分算子），返回的是隐式算子对象而非普通矩阵 A = cg.make_G(ig); % 正投影：得到无噪声正弦图 sinogram ybar = A * xtrue; % 模拟泊松计数噪声，1e5 是光子计数标度，越小噪声越大 rng(0); yi = poisson(ybar, 1e5);

这里最关键的一点：A不是一坨普通 double 矩阵。如果你在 Matlab 工作区里看它，会发现它是一个封装好的算子对象，内部可能存了扫描几何、插值表、甚至 mex 函数指针。你只需要记住两件事：A * x做正投影（从图像到正弦图），A' * y做反投影（从正弦图回到图像域）。它对外表现就是“乘一个矩阵”，但内部永远不会把系统矩阵实体化。千万别手贱去执行full(A)或double(A)，256 乘 256 的重建网格加上几万条射线，实体化以后内存直接爆炸。

poisson(ybar, 1e5)是 MIRT 里的噪声模拟函数：第一个参数是把无噪声投影值当作泊松分布的期望计数，第二个参数是光子数标度。这个标度越大，得到的投影噪声越小，重建图越接近无噪声理想结果；调成 1e4 甚至 1e3，信噪比会肉眼可见地下降。这是后面所有调参实验的基准变量，建议记下来。rng(0)是为了让噪声可复现，跑实验时保证对比公平。

3.4 FBP 基线与 OS-SQS 迭代重建：第一次对比

有了噪声正弦图，接下来就是重建的主菜。先用 FBP 做一个基线，再用 MIRT 的 OS-SQS 做迭代重建，两者直接对比：

% FBP 基线重建：滤波反投影，一条命令出图 xfbp = fbp2(yi, cg, ig); % 构造正则化对象：二次惩罚，beta 控制平滑强度 R = Robject(ig.mask, 'beta', 2^8, 'potential', 'quad'); % OS-SQS 迭代重建：20 轮迭代，12 个子集 [xos, info] = os_sqs(yi, A, R, 'niter', 20, 'nblock', 12); % 定量比较：计算掩膜区域内 RMSE err_fbp = xtrue - xfbp; err_os = xtrue - xos; rmse_fbp = sqrt(mean(err_fbp(ig.mask(:)).^2)); rmse_os = sqrt(mean(err_os(ig.mask(:)).^2)); fprintf('FBP RMSE = %g\n', rmse_fbp); fprintf('OS-SQS RMSE = %g\n', rmse_os);

fbp2内部做了斜坡滤波加反投影，一条命令出结果，速度极快，但它对含噪投影几乎无计可施——斜坡滤波器本质是高通滤波，会把高频噪声放大，所以 FBP 图上满屏的放射状条纹。Robject是 MIRT 传统的正则化构造器，ig.mask限定了计算区域，'beta', 2^8是正则强度，'potential', 'quad'表示二次势能函数，也就是经典的吉洪诺夫平滑。如果你的 MIRT 版本较新提示找不到Robject，用Reg1(ig.mask, 'beta', 2^8)等效替代。

os_sqs是 ordered subsets + separable quadratic surrogates 的缩写：SQS 把目标函数在每个像素上分解成可分代理函数，保证每次迭代目标值单调下降；OS 把 120 个投影视角分成 12 个子集，每个子集做一次子迭代，加速约等于子集数倍。niter是外层迭代轮数，nblock是子集数。运行之后打印的 RMSE 会告诉你一个典型结果：FBP 的 RMSE 可能是迭代重建的两到三倍。第一眼看到这个数字，就明白迭代重建为什么是今天临床 CT 的主流了。

4. 算法与参数选型：目标函数、噪声模型与 beta 的调节门道

4.1 目标函数长什么样：数据保真项与正则项怎么配对

MIRT 里所有迭代重建算法，本质上都在最小化同一个家族的目标函数：数据保真项加正则项。用数学语言写就是：

x̂ = argmin_x D(Ax, y) + β R(x)

D(Ax, y)是数据保真项，衡量“当前图像正投影后与实测投影的差距”；R(x)是正则项，惩罚图像的不合理波动；β（代码里的beta）负责平衡两者。理解这个结构比记住任何函数名都重要。数据保真项的选择取决于噪声模型：如果观测噪声是高斯型，用加权最小二乘；如果是泊松计数，应该用泊松对数似然。正则项的选择取决于你对图像的先验：二次惩罚平滑均匀但不保边缘，TV 项保边缘但容易出块状伪影，Huber 势能在两者之间折中。

MIRT 把这些选择全部参数化了。你可以不写一行求导代码，只通过组合参数就切换出十几种不同算法。我见过不少新手直接把os_sqs当成黑盒用，换了个数据就报错，原因就是没意识到这个函数内部默认了一套泊松保真项加二次惩罚的组合，换了 MRI 的复数数据自然不匹配。

4.2 噪声模型决定算法入口：Poisson 走 OS-SQS，Gaussian 走 CG

MIRT 的算法入口选择，第一判断依据不是“哪个函数新”，而是“你的数据噪声是什么类型”。CT 和 PET 的投影数据本质是光子计数，服从泊松分布，首选os_sqs或qpsqs，它们对泊松似然的处理是经过代理函数优化的，收敛稳定。MRI 的 k 空间数据噪声接近高斯分布，加权最小二乘模型就够，用共轭梯度类求解器更高效，比如qpwls_pcg。下表是我平时选择的快速参照：

这里要泼一盆冷水：不要以为os_sqs是万能药。它对泊松模型很擅长，但你如果非要用它解高斯噪声问题，本身也不会出错，只是绕了弯路。反过来，用高斯假设处理低计数 CT 投影，重建图会在暗区出现明显伪影，这是模型错误不是算法 bug。我自己的习惯是收到一批新数据先画直方图，看噪声分布形态，再决定走哪条求解路线。

4.3 beta 怎么调：从量纲到指数网格扫描

beta是 MIRT 里最影响结果、也最容易翻车的参数。β 设太大，图像被过度平滑，细小结构全没了；β 设太小，噪声完全没被压住，迭代重建退化成带噪声的最小二乘解。理想情况是取在“误差曲线拐点”附近。我从来不用手搓的方式猜 β，而是固定几何和迭代参数，在指数网格上扫一遍：

% 候选 beta：从 2^6 到 2^14，每个数量级都覆盖 betas = 2.^(6:2:14); rmse = zeros(size(betas)); for ii = 1:numel(betas) R = Robject(ig.mask, 'beta', betas(ii), 'potential', 'quad'); [x, info] = os_sqs(yi, A, R, 'niter', 20, 'nblock', 12); rmse(ii) = sqrt(mean((xtrue - x).^2, 'all')); end % 画出来，找误差最低点 plot(betas, rmse, 'o-'); set(gca, 'XScale', 'log');

为什么用2^n网格？因为 β 在目标函数里的作用接近指数尺度，线性网格扫 0 到 1000 基本是浪费计算量。扫完之后看 RMSE 曲线，通常在某个数量级区间内最低，两边都会升高。这个最低点对应的 β 就是用在这个数据上的推荐值。如果曲线整体平缓，说明你的数据噪声很低，正则项影响有限；如果曲线尖锐，说明问题高度病态，β 选择必须谨慎。

值得注意的是，β 的最优值跟投影光子数I0、视角数na、像素尺寸dx都有关系。换了一组扫描参数，原来调好的 β 就要重扫。我的习惯是每次做新的模拟实验，哪怕只改了噪声水平，也会把 β 扫描脚本重跑一遍，而不是凭记忆沿用旧值。这一步在避坑章节里还会细说。

5. 避坑：MIRT 使用里最常见的五个翻车现场

5.1 坑一：路径配置不全，类加载失败

现象：启动后执行image_geom(...)直接报错“Undefined class 'image_geom'”，或者提示找不到函数，但明明已经把解压目录加进路径了。

原因：第一，addpath只加了根目录，没有用genpath递归包含子目录，而image_geom的类定义文件在mirt下一层子模块里，根本不在搜索路径上。第二，路径里有中文、空格或特殊符号，Matlab 类加载机制解析失败。第三，路径虽然加了，没执行savepath，重启后配置丢失。

解决：统一用addpath(genpath('D:/tools/mirt'))，把 MIRT 移到纯英文路径下，执行后跑一下which image_geom，返回路径就对了。如果还不行，执行rehash toolboxcache强制刷新工具箱缓存。我自己现在的习惯是安装后立刻把这三条命令写成一行存进startup.m，以后换机器也不怕。

5.2 坑二：mex 编译失败卡住安装流程

现象：运行某个 demo 时提示找不到Gtomo2_dsc之类函数，或者直接报“Invalid MEX file”。上网一查说是要编译 mex，但mex命令报错 C 编译器没配好。

原因：MIRT 老版本里一部分系统矩阵算子依赖 C 加速，新版虽然大部分是纯 Matlab，但某些特定几何仍会调用 mex。编译失败通常是 Matlab 没找到匹配的编译器，或者装了 MinGW 后没有在 Matlab 里选中它。

解决：先执行mex -setup C，选好编译器后进到 MIRT 的mex目录，看有没有一键编译脚本，没有就找目录下哪些.c文件对应的函数缺失，逐个mex编译。其实在较新版本里，绝大多数算子已经不需要 mex，我遇到这种情况会先检查是不是版本太老，再决定是否编译，而不是一上来就装编译器。顺带一提，Matlab 对 MinGW 版本有严格匹配要求，缺编译器时先查官方支持列表，别装最新版 MinGW，装了也白装。

5.3 坑三：系统矩阵实体化，内存直接爆掉

现象：代码一执行到A = cg.make_G(ig)之后，Matlab 内存直线上升，几秒钟后卡死，甚至直接“内存不足”崩溃。有人还会手滑执行full(A)或double(A)，瞬间把内存吃光。

原因：认为 A 是普通矩阵，想把它转成 double 看看内容。MIRT 的系统矩阵是隐式算子，内部封装了投影和反投影逻辑，但它被设计成不存储完整矩阵。256 乘 256 重建网格配上几万条测量的系统矩阵，实体化以后是上百 GB 的稠密矩阵，什么机器都扛不住。

解决：把 A 当作不透明对象，永远只通过A * x和A' * y使用它。真要看某个局部内容，用A(:, idx)取特定列或A(idx, :)取特定行，这样才会按需计算那一部分。这条规矩我是在连续爆了两次内存之后才彻底记住的。如果你内存有限，系统矩阵相关的算子会主动把数据分成块，不需要你手动处理，你需要做的就是别去拆它的封装。

5.4 坑四：poisson 函数的尺度参数理解反了

现象：用poisson(ybar, 1e5)生成噪声投影，结果图像比无噪声还干净；调成 1e3 又噪声大得离谱，完全找不到规律。

原因：poisson的第二个参数是“光子计数标度”而不是噪声方差。它表示的是期望计数的整体缩放：第二个参数越小，相当于探测器接收到的光子越少，泊松噪声相对幅度越大。参数方向理解反了，就会觉得这个函数行为诡异。

解决：每次用之前看help poisson，确认参数语义。实际实验时以 1e4 到 1e6 作为常见范围去试。如果你要对比不同噪声水平，应该固定种子（rng(0)），这样投影噪声模式一致，只是强度不同，定量对比才有意义。我在自己的实验流程里会把光子数当成一个显式变量写进脚本头部，而不是埋在代码中间。

5.5 坑五：迭代收敛判断错，跑了一堆无效实验

现象：迭代重建跑完几十轮，图像肉眼看着还行，但 RMSE 没有下降，或者目标函数值还在明显波动。继续增加迭代次数也没有改善，白白浪费了 GPU 和一下午时间。

原因：第一，niter或nblock设置不匹配，比如视角数只有 120，却设了 60 个子集，子集近似误差太大，算法一直在振荡。第二，β 设置过小，目标函数几乎没有正则约束，迭代过程在噪声拟合阶段反复横跳。第三，只看了图像没看info结构体，误以为“图像没变=收敛了”。

解决：MIRT 的迭代函数返回的info结构体里有cost（目标函数值）和rmse等字段。判敛的标准不是图像不动，而是目标函数值曲线进入平台期。把plot(info.cost)画出来看，末尾应该是平坦的；如果还在单调下降或波动，说明没收敛。子集数nblock一般设 8 到 20，总子迭代数要覆盖全部视角至少一遍，也就是说niter * nblock不要小于na。这条是新手最容易忽略的细节，我看到过太多人拿着未收敛的结果下结论。

6. 进阶：用 info 结构监控迭代，把模拟流程迁到真实数据

6.1 用 info 结构盯住每一轮迭代

os_sqs之类的迭代函数返回的第二个参数info不是摆设，它是判断算法是否健康运行的唯一可信指标。我每次跑完重建第一件事不是看图，而是画 target 曲线：

% 跑 30 轮完整迭代，返回 info 结构 [xos, info] = os_sqs(yi, A, R, 'niter', 30, 'nblock', 12); % 画目标函数值与 RMSE，两者要同时看 figure; subplot(2, 1, 1); plot(info.cost); title('Objective'); subplot(2, 1, 2); plot(info.rmse); title('RMSE');

info.cost曲线如果先快速下降、后进入平台，说明算法收敛路径正常；如果中段有上翘甚至振荡，基本可以断定子集数nblock太大，或者 β 小到接近零。info.rmse反映真实误差变化，它比cost晚几步进入平台，但趋势应当一致。这条监控管线我现在每次都跑，它帮我节省过太多调试时间。

6.2 从模拟数据切到真实采集数据，要动哪些地方

换到真实数据时，模拟流程里需要替换的是“投影数据来源”，而不是重建框架。真实 CT 给你的通常是经过系统校准的 sinogram，可能不是原始计数，需要先做对数变换或归一化。ct_geom和image_geom要按真实扫描几何配置，探测器单元尺寸、视角数、像素大小都必须和扫描协议对齐。系统矩阵A可以复用同一套构建代码，但要确认几何参数正确，否则重建结果会有系统性形变。

真实数据没有xtrue，所以 RMSE 这条评估线只能退化成视觉评估和物理校验，比如检查已知密度区域的重建值是否落在理论范围内。这也是为什么我坚持先用模拟数据调好 β 和迭代参数再上真实数据——在可控环境里把参数学明白，真实数据上只需要微调，而不是从零摸黑。从那以后，我每接到一套新重建任务，都强制走一遍“定义几何 → 模拟数据 → FBP 基线 → 迭代重建 → 看 info 曲线”这五步，确定参数和算法行为正常之后，才把真实数据接进来。希望帮到你。

本文还有配套的精品资源，点击获取