CT断层成像系列01——Shepp-Logan头模型(附Matlab|C++代码实现)
时间地点_20260122上海同步辐射光源SSRF
序言
开门先介绍一本神级好书，我淘宝京东扫货了中文互联网所有的CT相关的书。医疗CT工业CT小动物CT等等甚至探矿地震CT，花了四五千，这本书是教的最细最深最好的，最关键的事，他每个章节每个算法都给了Matlab代码（有注释），同步辐射平行束CT，医疗扇束CT，工业锥束CT都有，解析法FBP迭代法ART甚至还有压缩感知CS. 这给了我很大的启发, 我发现世界上的知识但凡学不会就是因为课本不给例题不给答案不给公式不给源码不给老师亲手边实现边录的项目实战视频, 都形成行业了这些知识怎么会有门槛, 纯纯的防自学机制. 我在学习路上随时记录学习心得,并附上Matlab和C++两个版本的代码, 以方便后来的人可以轻松学会CT的知识并复现. 为西安电子科技大学的师生们狠狠点赞。桃李不言，下自成蹊。
Shepp-Logan头模型
下面开始正题
在图像重建算法的研究中，为了比较客观的评价各种算法的有效性，人们往往选择一种公认的模型作为研究对象。Shepp-Logan（S-L）头模型是目前CT领域仿真计算最常用的模型。
Shepp-Logan头是由L.A.Shepp和B.F.Logan于1974年首次提出的，它由10个位置、大小、方向、密度各异的椭圆叠加而成，来模拟一个医疗CT中脑部的断层。参见下图2.1(a)，其中英文字母代表10个椭圆的编号，数字则代表各个椭圆的密度 \(\rho\) 。图2.1(b)是断层模型的原始灰度图像。要想绘制椭圆，需要六个参数，他们分别是中心位置横坐标\(x_0\) ，中心位置纵坐标\(y_0\) ，椭圆长轴\(a\), 椭圆短轴\(b\), 旋转角度\(\alpha\), 椭圆密度\(\rho\). 值得注意的是, 最终的头模型的灰度值是由几个互相重叠的椭圆密度值累加而成.
我们可以直接调函数phantom在Matlab实现SL头:
clc;clear all;close all;
%自定义头模型的参数矩阵
shep=[
1 .69 .92 0 0 0
-.8 .6624 .8740 0 -.0184 0
-.2 .1100 .3100 .22 0 -18
-.2 .1600 .4100 -.22 0 18
.1 .2100 .2500 0 .35 0
.1 .0460 .0460 0 .1 0
.1 .0460 .0460 0 -.1 0
.1 .0460 .0230 -.08 -.605 0
.1 .0230 .0230 0 -.606 0
.1 .0230 .0460 .06 -.605 0
];
N=256;
I=phantom(shep,N); %生成自定义的Shepp-Logan头模型
%I=phantom(N); %生成默认的Shepp-Logan头模型
figure;imshow(I,[]);
OpenCV（C++ 版本）本身并没有像 MATLAB 那样内置的 phantom()函数 来直接生成 Shepp–Logan 头模型。因此没有可以直接调包的代码, 更下面有我手搓的C++代码, 写C++特征就是手搓轮子
值得注意的是, 参数力得长轴、短轴、中心点横纵坐标都是百分比的，这样你设置多大的图像，生成的头都是一个样。
运行结果如下:
不调包手搓Shepp-Logan头模型的数学实现
先附上两张插图方便对照理解, 尤其是图2.6和2.7, 后面就是描述2.7的三个坐标系变换:
画SL头就是画椭圆，对于一个规规整整坐落于中心的椭圆，绘图可用如下方程表示：
对于处在任意位置任意旋转的椭圆，我们要从公式（1）用三个坐标系变换来得到它。任意一个椭圆，以它的椭圆中心为坐标系\(X_1O'Y_1\)(这个坐标系相对于你是旋转的,歪的),其椭圆方程为：
坐标系（\(X_2O'Y_2\)）也过椭圆心，但是它旋转后与中央坐标系平行，所以椭圆要作旋转坐标系表达处理。有：
画个图就很好理解旋转坐标系的映射关系了:
将公式3带如公式2. 所以在\(X_2O'Y_2\)(这个坐标系相对于你是平直的,但是中心在椭圆心处不在原点), 椭圆方程为:
现在要从\(X_2O'Y_2\)回到中央坐标系\(XOY\), 需要坐标系平移, 有:
将(5)代入(4), 我们就得到了任意椭圆在中央坐标系\(XOY\)下的最终方程:
此公式代入(1), 我们就能画出Shepp_Logan模型头了:
公式推导完毕, 下面上C++的代码实现, 记得先配好OpenCV:
#include <iostream>
#include <opencv2/opencv.hpp>
using namespace cv;
using namespace std;
Mat phantomSheppLogan(int size = 256)
{
Mat img(size, size, CV_32F, Scalar(0));
//直角坐标系中心center(c),配合图像尺寸scale,将从左上角计数的坐标轴转换为中心百分比坐标轴
const int c = size / 2;
const float scale = size / 2.0f;
// 标准 Shepp-Logan 参数
//椭圆Ellipse绘图的六个参数 ρ密度|振幅 A长轴 B短轴 x0椭圆心横坐标 y0椭圆形纵坐标 α旋转角度
struct Ellipse { float rho, A, B, x0, y0, alpha; };
const Ellipse E[] = {
{1.0f, .69f, .92f, 0.f, 0.f, 0.f}, // 大椭圆
{-.8f, .6624f, .874f, 0.f, -.0184f, 0.f}, // 内椭圆
{-.2f, .11f, .31f, .22f, 0.f, -18.f * CV_PI / 180},
{-.2f, .16f, .41f, -.22f, 0.f, 18.f * CV_PI / 180},
{.1f, .21f, .25f, 0.f, .35f, 0.f},
{.1f, .046f, .046f, 0.f, .1f, 0.f},
{.1f, .046f, .046f, 0.f, -.1f, 0.f},
{.1f, .046f, .023f, -.08f, -.605f, 0.f},
{.1f, .023f, .023f, 0.f, -.605f, 0.f},
{.1f, .023f, .046f, .06f, -.605f, 0.f}
};
for (auto& e : E)
{
float ct = cos(e.alpha), st = sin(e.alpha);
for (int v = 0; v < size; ++v)
{
for (int u = 0; u < size; ++u)
{
//变换坐标系居中,再缩放比例至0-1,以使得所有尺寸的图绘制的一样
float x = (u - c) / scale;
float y = (v - c) / scale;
float x2 = x - e.x0;
float y2 = y - e.y0;
//公式(6),算这个像素与椭圆心的距离
float d2 = (x2 * ct + y2 * st) * (x2 * ct + y2 * st)
/ (e.A * e.A) + (-x2 * st + y2 * ct) * (-x2 * st + y2 * ct) / (e.B * e.B);
//如果距离小于1,本像素在椭圆内, 这个像素画上椭圆的密度,见公式(7)
if (d2 <= 1.0f)
img.at<float>(v, u) += e.rho;
}
}
}
normalize(img, img, 0, 1, NORM_MINMAX);
Mat img_updownflip;
flip(img, img_updownflip, 0);
return img_updownflip;
}
int main()
{
Mat phantom = phantomSheppLogan(1024);
imshow("Shepp-Logan Phantom", phantom);
waitKey(0);
}
代码运行结果如下:
下面附上Matlab的代码实现:
P = phantomSheppLogan(1024);
imshow(P, []); title('Shepp-Logan Phantom');
function P = phantomSheppLogan(N)
% P = phantomSheppLogan(N) 生成 N×N 的 Shepp-Logan 头模型
% 返回单精度图像，灰度范围 0~1
if nargin < 1, N = 256; end
% ---------- 1. 预定义 10 个椭圆 ----------
% 每行：[rho, A, B, x0, y0, alpha(度)]
ellipseParam = [ ...
1.0 0.69 0.92 0.0 0.0 0;
-0.8 0.6624 0.874 0.0 -0.0184 0;
-0.2 0.11 0.31 0.22 0.0 -18;
-0.2 0.16 0.41 -0.22 0.0 18;
0.1 0.21 0.25 0.0 0.35 0;
0.1 0.046 0.046 0.0 0.1 0;
0.1 0.046 0.046 0.0 -0.1 0;
0.1 0.046 0.023 -0.08 -0.605 0;
0.1 0.023 0.023 0.0 -0.605 0;
0.1 0.023 0.046 0.06 -0.605 0];
% ---------- 2. 坐标网格 ----------
% 将像素坐标归一化到 [-1,1]，原点位于图像中心
[y, x] = ndgrid( (0:N-1)-N/2+0.5 , (0:N-1)-N/2+0.5 );
x = x / (N/2); % 归一化
y = y / (N/2);
P = zeros(N, 'single'); % 累加密度
% ---------- 3. 逐椭圆叠加 ----------
for k = 1:size(ellipseParam,1)
rho = ellipseParam(k,1);
A = ellipseParam(k,2);
B = ellipseParam(k,3);
x0 = ellipseParam(k,4);
y0 = ellipseParam(k,5);
alpha = deg2rad( ellipseParam(k,6) ); % 转弧度
% 旋转矩阵
ca = cos(alpha); sa = sin(alpha);
% 相对椭圆中心的坐标
xx = x - x0;
yy = y - y0;
% 旋转后坐标
xr = xx*ca + yy*sa;
yr = -xx*sa + yy*ca;
% 归一化距离平方
d2 = (xr/A).^2 + (yr/B).^2;
% 在椭圆内的像素累加密度
mask = d2 <= 1;
P(mask) = P(mask) + rho;
end
% ---------- 4. 0-1 归一化 + 上下翻转 ----------
P = (P - min(P(:))) / (max(P(:)) - min(P(:)));
P = flipud(P); % 与 C++ 版保持一致
end
运行结果如下:
结束.
浙公网安备 33010602011771号