一个做图像重建的朋友问我phantom函数里那些椭圆参数到底是怎么控制图像的他说自己在网上找了一圈大部分帖子都在教怎么调用phantom(Shepp-Logan)但一涉及到自定义参数矩阵全都语焉不详。这问题我太熟悉了因为我刚开始用Matlab做CT成像仿真时也被这个函数绕得晕头转向。后来花了两个晚上把它的参数矩阵挨个调了一遍才算真正搞明白这个函数在干什么。今天这篇就专门讲透Matlab的phantom函数——从它最基础的原理到Shepp-Logan模型的来龙去脉再到如何自己写参数矩阵构造一个包含颅骨、脑组织、病灶等结构的自定义头部仿真。本文面向所有做图像重建、CT仿真、医学图像处理以及想用这个函数做算法验证的读者我会把每一步的操作细节和踩过的坑都写清楚。1. phantom函数解决的是什么问题为什么图像重建离不开标准模体1.1 没有标准答案算法好坏就无从谈起先想一个问题你开发了一个新的滤波反投影算法想验证它重建出来的图像是不是准确拿什么图来测如果你直接拿一张真实的CT图像做投影再重建看起来好像很真实但是问题在于——你并不知道这张图在理论上的完美重建应该长什么样。真实图像本身就存在噪声、伪影和各种不确定性拿它当测试基准最后算出来的PSNR、SSIM指标其实没有一个可对照的真值。这就是标准模体存在的意义。模体phantom是一张我们预先知道精确灰度分布的人造图像我们可以对这张图做正向投影比如Radon变换得到正弦图然后让重建算法从正弦图里反推图像再和原始模体比较。因为真值就在手里算法的误差就能被精确测量。phantom函数生成的就是这种标准模体。它不是照片不是真实扫描结果而是用若干个几何形状标准情况下是椭圆按不同灰度和位置叠加起来的合成图像。最经典的用途是模拟人体头部X射线CT成像时的横断面所以也常被叫做头部模体。1.2 Shepp-Logan头模型的诞生与改进phantom函数最原始的版本对应的是Shepp和Logan在1974年提出的头模型。他们在研究计算机断层扫描重建算法时构造了一个由10个椭圆组成的图像用来模拟人类头部的横截面。这些椭圆分别代表颅骨、脑组织、肿瘤、血肿等不同密度区域。由于椭圆在数学上很容易通过解析表达式描述它的投影线积分也能用解析方法快速计算因此在验证重建算法时非常高效。后来人们发现原始的Shepp-Logan模体对比度太低灰度范围很窄直接显示出来几乎看不清细节。所以在1980年前后出现了Modified Shepp-Logan版本主要调整了各个椭圆的灰度值把相对差异放大同时依然保持解剖结构的相对关系让视觉效果更适合观察和量化分析。现在Matlab里的phantom函数如果你直接用phantom(n)或者phantom(Modified Shepp-Logan, n)默认输出的其实是后者这个改良版。1.3 打开Matlab第一次看到phantom的感觉我第一次跑phantom(256)这行代码时输出一张灰底上有个椭圆形头部轮廓的图像里面有几个稍亮稍暗的斑块。当时觉得这东西不就是几个椭圆的堆叠吗后来才知道正是这几个椭圆的堆叠因为它的每一个像素灰度值都能被精确复现才让它成为图像重建领域横跨几十年的金标准。你甚至可以把phantom的输出当成一个理想的数字人体断面所有的高精度重建论文里几乎都少不了这张图。所以phantom函数不只是一个画几个椭圆的玩具它是连接解析几何、成像物理和重建算法的桥梁。理解了它的参数矩阵你就等于掌握了一套自己构造仿真模型的工具这也是本文要深入下去的核心。2. 读懂phantom的核心参数矩阵才是控制图像的唯一入口2.1 五种调用方式与默认参数Matlab的phantom函数有好几种调用方式先把它们列全P phantom生成默认的256×256的Modified Shepp-Logan模体。P phantom(n)生成n×n的默认模体。P phantom(256)等价于上面。P phantom(E, n)通过自定义参数矩阵E配合n×n的图像尺寸生成模体。[P, E] phantom不仅返回图像还返回默认的参数矩阵E。这个非常有用因为你可以拿到E之后改一改再传给phantom实现自定义。需要特别注意的是phantom函数输出的是一个double类型的二维矩阵而不是通常意义上0到255的灰度图像。也就是说它的每个像素点都可能是0.2、0.9这样的小数。默认背景是0各个椭圆区域的灰度值不同从低到高有以下几种0、0.2、0.3、0.4、1等。所以你直接调用imshow(P)看到的效果可能会是一片暗沉沉的图——因为imshow默认把double矩阵当作0到1的数来显示如果矩阵里有负值或大于1的值显示还会出问题。2.2 矩阵里每一列到底代表什么现在来看phantom函数的精髓——参数矩阵E。默认的E是一个10行6列的矩阵每一行代表一个椭圆。官方文档里对每一列的定义如下列号含义第1列椭圆中心的x坐标相对于图像中心范围约为-1到1第2列椭圆中心的y坐标第3列椭圆长轴半径长度单位为图像宽度的一半第4列椭圆短轴半径第5列椭圆长轴与x轴正方向的旋转角度单位度第6列椭圆区域的灰度值这个坐标系是以图像中心为原点的归一化坐标整个图像的横纵坐标范围大致在[-1, 1]之间。图像横向对应x方向纵向对应y方向。所有椭圆的数值都是相对单位比如长轴半径为0.69就表示它的半长轴约为图像宽度一半的0.69倍。如果你用[P,E] phantom然后disp(E)会得到这样一个矩阵数值可能随版本略有差异但结构一致0 0 0.6900 0.9200 0 0 0 -0.0184 0.6624 0.8740 0 -0.0200 0.22 0.1100 0.1100 0.3100 -0.1800 -0.0200 0.22 0.1100 0.1600 0.4100 -0.1800 -0.0200 0 0 0.2100 0.2500 0 0.0100 0 0 0.0460 0.0460 0 0.0100 0 -0.3500 0.0460 0.0230 0 0.0100 -0.0800 0 0.0460 0.0230 0 0.0100 0 -0.0800 0.0460 0.0230 0 0.0100 -0.0800 0.6500 0.0460 0.0230 -0.1800 0.0100第一行是最外层的大椭圆长轴0.69、短轴0.92注意短轴数值大于长轴这说明这个椭圆其实是竖着放的y方向更长。第二行是颅骨轮廓灰度为-0.02负值代表比背景更暗的区域。后面几行分别模拟脑组织、肿瘤和血肿等结构。2.3 用一个小实验验证矩阵与图像的对应关系为了验证理解是否正确你可以做一个小实验。把默认矩阵E复制出来只保留第一行其余行全部删掉然后调用phantom(E,256)。猜一下结果会是什么应该是一个只有一个椭圆没有内部结构的模体。但如果你直接运行会发现——图像变成了一个里面什么都没有的椭圆灰度值统一等于第一行的灰度0。背景为0椭圆内部为0那图像看起来就是全黑的这里有个容易混淆的点如果椭圆内部灰度设为0和背景0相同那你看不出任何形状。这正好说明参数矩阵里的灰度值可以任意取并不强制要求图像前景比背景亮。在默认的Shepp-Logan模体中背景是0颅骨对应的灰度是-0.02比背景暗脑组织等区域是0.01或0.02比背景亮通过正负灰度差来区别不同组织。把矩阵多保留几行你就能看到形状一层层叠加出来。这正是phantom函数的工作方式从背景开始依次绘制每个椭圆用对应灰度值填充椭圆内部。绘制顺序就是矩阵行的顺序后出现的椭圆会覆盖先出现的椭圆——这一点在做自定义时非常重要因为你要是把颅骨写在脑组织后面脑组织就会把颅骨区域覆盖掉。3. 自定义头部仿真亲手搭建颅骨、组织与病灶3.1 确定头部仿真的解剖结构清单理解了参数矩阵后我们完全可以自己设计一个更符合需求的头部模体。假设你现在想模拟一个经过简化的头部横断面包含以下结构颅骨最外层的闭合椭圆环用两个椭圆相减来模拟但phantom只支持填充不支持打孔——所以更常见的是用一个较暗的椭圆做颅骨背景再在里面填充一个较大的亮椭圆做脑组织这样两者之间形成环形边界。脑组织位于颅骨内侧的大椭圆。脑脊液区域在脑组织和颅骨之间或内部亮度略低的结构。肿瘤/病灶一个或多个小的亮斑或暗斑。血肿一般用一个小的高灰度椭圆模拟。如果想让图像更丰富还可以加入一些组织如白质、灰质的分界用两个相邻的不同灰度的椭圆或者模拟钙化点一个小而亮的区域。3.2 从解剖坐标到椭圆参数的计算与转换如果你熟悉解剖坐标想设计一个更精确的头部模体就需要把实际毫米尺寸转换成phantom的归一化坐标。假设整个头部横断面左右宽度约160mm前后长度约200mm图像中心为原点那么x方向归一化坐标 实际x(mm) / 80y方向归一化坐标 实际y(mm) / 100。椭圆的半径也按同样的比例换算。举个例子一个肿瘤位于左侧x-20mm中心略偏前y10mm直径约10mmx方向和8mmy方向那么它的参数就是x-0.25y0.1长轴半径0.0625短轴半径0.05角度根据形状方向设定如果轴对齐则角度为0灰度根据你要模拟的CT值设定。但说实话设计仿真模体时不需要太纠结于绝对毫米精度更关键的是结构之间的空间关系、对比度和尺寸比例。你需要让这些椭圆有足够的几何区分度保证后续做投影和重建时各个结构的边界是清晰的。3.3 写一个自己的phantom生成函数我们可以直接基于默认矩阵E来修改而不是从零写一个复杂的参数矩阵。这样既保留了Shepp-Logan的整体解剖结构又能加入自定义病变。思路如下% 获取默认参数矩阵 [~, E] phantom; % 新增一个肿瘤中心(-0.3, 0.1)半径0.05 x 0.04角度30度灰度0.8 newEllipse [-0.3, 0.1, 0.05, 0.04, 30, 0.8]; E2 [E; newEllipse]; % 生成自定义模体 P_custom phantom(E2, 512); figure; imagesc(P_custom); axis image; colormap gray; colorbar;运行之后你会在原本的头部图像中看到一个额外的亮斑。这个操作看起来很简单但如果你不把这个椭圆的灰度设得比周边组织高很多视觉上可能不明显。所以做自定义时要清楚每个结构你想让它亮还是暗以及它应该出现在哪一层。如果不想改动默认矩阵而是完全自己构造矩阵你可以从最简单的三行开始颅骨外环、脑组织、一个小病灶。比如E3 [ 0, 0, 0.9, 1.2, 0, 0; % 外层颅骨灰度0暗 0, 0, 0.8, 1.0, 0, 0.9; % 脑组织灰度0.9亮 0.2, -0.3, 0.1, 0.1, 45, 1.5; % 肿瘤亮斑 ]; P3 phantom(E3, 256);这种写法直观容易控制。我建议新手从这种全自定义开始每加一行就显示一次观察变化慢慢叠加出自己的头部模型。4. 把静态模体变成临床场景噪声、对比度与CT投影4.1 调整灰度分布模拟不同组织对比度默认的Modified Shepp-Logan虽然看起来很美但灰度值范围在-0.02到1之间各组织之间的对比度是固定的。做重建实验时你可能希望调整组织的CT值差。比如颅骨在CT图像中应该是高亮高衰减脑组织是中等而空气是低值。但Shepp-Logan原始模型为了显示方便把颅骨设为负值这就和真实CT值有出入。matlab的phantom函数本身不做任何物理意义层面的保证你可以随手把灰度改为任意值。比如把颅骨改成2.0脑组织改为1.0肿瘤改为3.0背景改为0这样模体呈现出的就是骨骼亮、组织暗的临床风格。需要注意的是phantom输出矩阵的值都是线性比例你后续如果要用它模拟CT的HU值得自己再加一个线性映射。比如把你想要的HU值除以1000再赋给对应椭圆。4.2 添加噪声与分辨率评估真实的CT图像总会有噪声用纯净的phantom测试算法虽然理想但和实际差距较大。通常的做法是往phantom里添加高斯噪声或泊松噪声。比如P_noise imnoise(P, gaussian, 0, 0.01);但要注意imnoise默认期望输入范围是[0,1]的双精度图像如果你的p值有负数建议先归一化再加噪否则噪声强度不一致。还有一种做法是手动加高斯噪声sigma 0.05; P_noise P sigma * randn(size(P));模拟低剂量CT时泊松噪声更符合物理过程可以先用radon得到投影再在投影域加泊松噪声最后重建。这样更接近真实成像链。phantom还能用来评价重建算法的分辨率。因为你确切地知道模体中小椭圆的尺寸和位置重建后如果这些小肿瘤被模糊成一片或者完全消失说明算法的空间分辨率不够。用不同类型和尺寸的椭圆添加多个小目标就能在模体里做分辨率测试卡。4.3 与radon配合生成正弦图并做滤波反投影重建phantom和radon是Matlab里做CT仿真最常用的搭档。先对phantom做Radon变换模拟各个角度下的X射线投影再用不同算法从正弦图重建图像最后和原始phantom比较。P phantom(256); theta 0:179; % 投影角度0到179度 [R, xp] radon(P, theta); % 用滤波反投影重建 I iradon(R, theta, linear, Ram-Lak); figure; subplot(1,3,1); imshow(P, []); title(原始模体); subplot(1,3,2); imshow(R, [], XData, theta, YData, xp); axis normal; xlabel(角度); ylabel(探测器位置); title(正弦图); subplot(1,3,3); imshow(I, []); title(滤波反投影重建);这个流程是图像重建实验的经典套路。你还可以改变投影角度数量比如只给90个角度、加噪声、换滤波器如Shepp-Logan、Hamming然后用PSNR、SSIM指标来定量评估重建质量。这时候phantom的价值体现得淋漓尽致因为原始图像已知你可以把误差精确到每个像素。用自定义参数矩阵生成的模体同样可以直接喂给radon。radon函数不关心你的图像是从哪里来的它只做数学投影。所以只要你构造出任何想要的图像都能用它做仿真。5. 实战避坑笔记phantom使用中经常翻车的五个细节5.1 显示问题imshow和imagesc的灰度尺度这是最容易被坑的地方。phantom返回的是double矩阵值范围通常不是标准的0-255。如果你直接用imshow(P)由于imshow对double类型的默认显示范围是[0,1]而Shepp-Logan模体里存在负值这会导致负值显示为黑色所有低灰度细节全部丢失。更稳妥的方法是用imagesc(P)加colorbar或者imshow(P, [])——方括号表示让Matlab自动把最小值映射为黑色、最大值映射为白色。但注意如果你后面要定量比较不同模体的视觉效果用imshow(P, [])会因为每张图的最大最小值不同而产生不同的映射导致视觉对比不一致。这时建议统一显示范围比如imshow(P, [-0.1, 1.2])。5.2 翻转与旋转坐标轴方向对结果的影响phantom函数内部使用的坐标是x水平、y垂直角度是逆时针方向为正。但Matlab的图像矩阵索引是先行后列也就是第一维是y方向、第二维是x方向。当你自己设计椭圆参数时很容易把x、y搞反。比如你想把椭圆中心放在图像的上方y为正但在矩阵里你写成了第2列为负——结果椭圆跑到了下方。另外国家标准图像坐标通常是y轴向下而phantom的数学坐标是y轴向上。这意味着如果你在matlab里用imagesc(P)显示图像的上半部分在数组中是靠前的行也就是实际y坐标为正的部分对应显示在图像上方。这点在旋转椭圆时特别容易混乱。我的建议是在纸上先画出坐标轴标出每个椭圆的中心和轴方向再转换成矩阵行多做几次后自然就熟练了。5.3 自定义椭圆越界时图像可能完全变样phantom里椭圆的坐标和半径是相对值但不强制约束椭圆必须在[-1,1]范围内。如果你把椭圆中心放到(0.8, 0.8)长轴0.6这超出了图像边界Matlab会照常渲染——超出部分自然被截断图像里只出现一部分椭圆。这本身不是什么错误但如果你没意识到这一点可能会被怎么多了一块灰白色的新月形区域搞懵。在设计自定义模体时建议先用简单的检查脚本判断每个椭圆是否都在有效范围内for i 1:size(E,1) center E(i,1:2); radii E(i,3:4); xlim center(1) radii(1) * [-1 1]; ylim center(2) radii(2) * [-1 1]; if xlim(1) -1 || xlim(2) 1 || ylim(1) -1 || ylim(2) 1 warning(第%d个椭圆越界, i); end end5.4 别把phantom的输出当真图像它的灰度范围是抽象的phantom函数的输出仅仅是一个数学表达式的离散化结果它不代表真实的CT值也不代表某种物理衰减系数。很多初学者把phantom直接当作模拟CT图像来展示但实际它更像一个理想的测试图案。如果你需要模拟16位CT图像那种HU值范围比如-1000到3000需要自己做一个线性映射并转换为uint16。一个常见操作是P phantom(256); P_hu P * 3000 - 1000; % 映射到-1000到2000范围粗略示意 P_uint16 uint16(P_hu 1024); % 加上偏移并转无符号整数但要注意这种映射完全是你自己定义的phantom不负责保证物理正确性。所以不要在你的论文里说使用phantom生成了临床CT数据而应该说使用数值模体仿真了CT成像过程。5.5 随机数种子与噪声叠加的稳定性如果你在仿真实验中需要保证结果可重复叠加随机噪声之前一定要设置随机数种子。Matlab不同版本用的函数不一样现在推荐用rngrng(2024); % 设置种子 P_noise P 0.05 * randn(size(P));如果不设置种子每次运行生成的噪声都不同实验结果就无法复现。这在写论文或者做多人协作时是灾难——你辛苦调出来的参数别人那边跑出来完全是另一组数据。6. 从phantom出发的进阶玩法3D扩展与算法评测6.1 扩展到3Dphantom3m的用法基础phantom函数只能生成2D断面但很多研究如锥形束CT重建、三维图像配准需要三维模体。好在Matlab还提供了phantom3m这个函数在某些版本中可能位于自定义工具箱或Image Processing Toolbox相关文件里它能生成3D的Shepp-Logan头部模型返回一个三维数组。P3 phantom3m(128); % 生成128x128x128的三维头模 % 查看中间切面 figure; imagesc(squeeze(P3(64,:,:))); axis image; colormap gray;如果你手头没有phantom3m也可以自己写一个生成3D椭球的模体函数。3D椭球的参数比2D多一组中心x,y,z三个半径以及绕三个轴的旋转角度。通过在不同切片上画不同的2D椭圆可以近似构建一个3D头部模型。不过这种做法效率不高推荐研究一下phantom3m内部的实现逻辑体会一下作者是怎么把多个椭球叠加成三维图像的。6.2 用phantom做重建算法的压力测试phantom真正的大用途是作为算法评测的基准。你可以通过修改参数矩阵在模体里放置一系列从大到小的椭圆用来测试重建算法对细节的分辨能力也可以把椭圆的对比度设得非常低用来测试算法对低对比度物体的检测灵敏度。比如你设计一个包含5个不同灰度梯度的椭圆组重建后再测量每个区域的灰度均值看看算法是否引入了伪影或偏移。一个标准的评测流程应该是用自定义矩阵生成高分辨率模体比如1024×1024确保结构清晰。对其进行radon变换得到投影数据。对投影数据加一定强度的高斯噪声。用测试算法重建得到1024×1024的重建图像。将重建图像降采样到256×256与256×256的原始模体对比计算RMSE、PSNR、SSIM等指标。改变角度采样数、噪声水平、重建滤波器画出精度曲线。这样的一套流程能非常客观地反映一个重建算法的性能。你能确切地知道每个误差来源对最终结果的贡献。这些都依赖phantom提供一个可精确复现的标准图像。6.3 结合图像处理和深度学习的仿真数据增强如果你正在做深度学习图像重建相关的研究phantom同样是一个非常好用的数据生成器。你可以在基础模体上做大量变换合成大量的真值图-投影数据对用来训练网络。比如对自定义模体做不同角度的旋转、缩放、平移、添加不同形状的病变快速生成上千个带标签的训练样本。尽管phantom的结构比较简单但它能提供严格对齐的真值这在监督训练里非常宝贵。更进一步你还可以把phantom的输出作为基础用插值和形变场让它更接近真实解剖数据的形态。我之前做过一个实验在Shepp-Logan模体中加入10个随机的椭圆每个椭圆的位置、大小、角度和灰度都随机然后生成批量训练数据。本来以为这个做法太简单但实际测试下来模型在真实CT数据上的泛化能力确实有提升。这说明phantom虽然老但在现代研究中的价值依然不可替代。最后再分享一个我自己的小习惯每次拿到一个新的重建算法我不急着用真实数据跑而是先在phantom上做一遍全流程调试参数、评估效果。因为phantom能将变量控制到最小一旦发现异常我会先检查算法本身而不是怀疑模体数据有问题。这个习惯帮我避开了很多假象伪影干扰。如果你也想深入学习和使用Matlab图像处理仿真从phantom开始绝对是一个性价比极高的选择。