《数学建模与数学实验》教学教学资源（案例）CT图像的代数重建问题

CT图像的代数重建问题 X射线透视可以得到3维对象在2维平面上的投影，CT则通过不同角度的X 射线得到3维对象的多个2维投影，并以此重建对象内部的3维图像.代数重建 方法就是从这些2维投影出发，通过求解超定线性方程组，获得对象内部3维图 像的方法 图11双层螺旋CT 图12CT图像 这里我们考虑一个更简单的模型，从2维图像的1维投影重建原先的2维图 像.一个长方形图像可以用一个横竖均匀划分的离散网格来覆盖，每个网格对应 一个像素，它是该网格上各点像素的均值.这样一个图像就可以用一个矩阵表示， 其元素就是图像在一点的灰度值（黑白图像），.下面我们以3×3图像为例来说明. 表4消耗与产出情况 3×3图像 水平方向上 各点的灰度值 的登加值 x1=12=0r3=0 r1++r3=1 x4=0r=05x=05 7=0.5x8=0x9=1 7+8+9=1.5 竖直方向上 1+4十72+5十83+6十和 的叠加值 =1.5 =0.5 =1.5 每个网格中的数字x代表其灰度值，范围在0,1]内.0表示白色，1表示黑色，0.5 表示灰色加如果我们不知首网格中的数值只知首沿竖直方向和水平方向的叠加 值，为了确定网格中的灰度值，可以建立线性方程组（含有6个方程，9个未知数） [x+x2+x3=1 X+xg+x。=1 x,+X+x。=1 显然该方程组的解是不唯一的，为了重建图像，必须增加叠加值.如我们增加从 右上方到左下方的叠加值，则方程组将增加5个方程 X1=1, x2+X4=0, X+x5+7=1 x6+x8=0.5, X9= 和上面的6个方程放在一起构成一个含有11个方程，9个未知数的线性方程组

CT 图像的代数重建问题 X 射线透视可以得到 3 维对象在 2 维平面上的投影, CT则通过不同角度的 X 射线得到 3 维对象的多个 2 维投影, 并以此重建对象内部的 3 维图像. 代数重建 方法就是从这些 2 维投影出发, 通过求解超定线性方程组, 获得对象内部 3 维图 像的方法. 图 11 双层螺旋 CT 图 12 CT 图像 这里我们考虑一个更简单的模型, 从 2 维图像的 1 维投影重建原先的 2 维图 像. 一个长方形图像可以用一个横竖均匀划分的离散网格来覆盖, 每个网格对应 一个像素, 它是该网格上各点像素的均值. 这样一个图像就可以用一个矩阵表示, 其元素就是图像在一点的灰度值(黑白图像). 下面我们以 33 图像为例来说明. 表 4 消耗与产出情况 33 图像 各点的灰度值 水平方向上 的叠加值 x1 = 1 x2 = 0 x3 = 0 x1 + x2 + x3 = 1 x4 = 0 x5 = 0.5 x6 = 0.5 x4 + x5 + x6 = 1 x7 = 0.5 x8 = 0 x9 = 1 x7 + x8 + x9 = 1.5 竖直方向上 的叠加值 x1 + x4 + x7 = 1.5 x2 + x5 + x8 = 0.5 x3 + x6 + x9 = 1.5 每个网格中的数字 xi 代表其灰度值, 范围在[0, 1]内. 0 表示白色, 1 表示黑色, 0.5 表示灰色. 如果我们不知道网格中的数值, 只知道沿竖直方向和水平方向的叠加 值, 为了确定网格中的灰度值, 可以建立线性方程组(含有 6 个方程, 9 个未知数) 1 2 3 4 5 6 3 6 9 1 1 1 x x x x x x x x x  + + =  + + =    + + = 显然该方程组的解是不唯一的, 为了重建图像, 必须增加叠加值. 如我们增加从 右上方到左下方的叠加值, 则方程组将增加 5 个方程 x1 = 1, x2 + x4 = 0, x3 + x5 + x7 = 1, x6 + x8 = 0.5, x9 = 1, 和上面的 6 个方程放在一起构成一个含有 11 个方程, 9 个未知数的线性方程组

【模型准备】设3x3图像中第一行3个点的灰度值依次为x1,第二行3个 点的灰度值依次为x4,x5,稀，第三行3个点的灰度值依次为7，x8,9.沿竖直方向 的叠加值依次为1.5,0.5,1.5，沿水平方向的叠加值依次为1,1,1.5，沿右上方到 左下方的叠加值依次为1,0,1,0.5,1.确定x1,x2, ,的值 【模型建立】由已知条件可得（含有11个方程，9个未知数的）线性方程组 +名+书= x4++x6=1 (-1 【模型求解】在Matlab命令窗口输入以下命令 >A=[1,1,1,0,00,0,0,0:0,0,0,1,1,1,0,0,0:0,0,0,0,0,0,1,1,1: 1,0,0,1,0,0,1,0,0:0,1,0,0,1,0,0,1,0:0,0,1,0,0,1,0,0,1 1,0,0,0,0,0,0.0,0:0,1,0,1,0,0,0,0,0:0,0,1,0,1,0,1,0,0 0,0,0,0,0,1,0,1,0:0,0,0,0,0,0,0,0,1]: >b=[1;1:1.5;1.5:0.5;1.5;1:0:1;0.5:1 >x A\b;x' MatIab执行后得 Warning:Rank deficient,rank =8 tol =4.2305e-015 ans 1.00000.00000-0.00000.50000.50000.5000-0.0000 1.0000 可见上述方程组的解不唯一，其中的一个特解为 x1=1,2=03=0,4=0,x=0.5,6=0.5,x7=0.5,x8=0,9=1. 【模型分析】上述结果表明，仅有三个方向上的叠加值还不够.可以再增加从左上 方到右下方的叠加值.在实际情况下 由于测量误差， 述线性方程组可能是超 定的.这时可以将超定方程组的近似解作为重建的图像数 Matlab实验题 给定一个3x3图像的2个方向上的灰度香加值，沿左上方到右下方的灰度叠 加值依次为0.8,12,17,02,03，沿右上方到左下方的灰度叠加值依次为0.6,02 16,12,0 ()建立可以确定网格数据的线性方程组，并用Matlab求解 (2)将网格数据乘以256，再取整，用Matlab绘制该灰度图像

【模型准备】设 33 图像中第一行 3 个点的灰度值依次为 x1, x2, x3, 第二行 3 个 点的灰度值依次为 x4, x5, x6, 第三行 3 个点的灰度值依次为 x7, x8, x9. 沿竖直方向 的叠加值依次为 1.5, 0.5, 1.5, 沿水平方向的叠加值依次为 1, 1, 1.5, 沿右上方到 左下方的叠加值依次为 1, 0, 1, 0.5, 1. 确定 x1, x2, ., x9 的值. 【模型建立】由已知条件可得(含有 11 个方程, 9 个未知数的)线性方程组 1 2 3 4 5 6 9 1 1 1 x x x x x x x  + + =  + + =    = 【模型求解】在 Matlab 命令窗口输入以下命令 >> A = [1,1,1,0,0,0,0,0,0;0,0,0,1,1,1,0,0,0;0,0,0,0,0,0,1,1,1; 1,0,0,1,0,0,1,0,0;0,1,0,0,1,0,0,1,0;0,0,1,0,0,1,0,0,1; 1,0,0,0,0,0,0,0,0;0,1,0,1,0,0,0,0,0;0,0,1,0,1,0,1,0,0; 0,0,0,0,0,1,0,1,0;0,0,0,0,0,0,0,0,1]; >> b = [1;1;1.5;1.5;0.5;1.5;1;0;1;0.5;1]; >> x = A\b; x’ Matlab 执行后得 Warning: Rank deficient, rank = 8 tol = 4.2305e-015. ans = 1.0000 0.0000 0 -0.0000 0.5000 0.5000 0.5000 -0.0000 1.0000 可见上述方程组的解不唯一. 其中的一个特解为 x1 = 1, x2 = 0, x3 = 0, x4 = 0, x5 = 0.5, x6 = 0.5, x7 = 0.5, x8 = 0, x9 = 1. 【模型分析】上述结果表明, 仅有三个方向上的叠加值还不够.可以再增加从左上 方到右下方的叠加值. 在实际情况下, 由于测量误差, 上述线性方程组可能是超 定的. 这时可以将超定方程组的近似解作为重建的图像数据. Matlab 实验题 给定一个 33 图像的 2 个方向上的灰度叠加值: 沿左上方到右下方的灰度叠 加值依次为 0.8, 1.2, 1.7, 0.2, 0.3; 沿右上方到左下方的灰度叠加值依次为 0.6, 0.2, 1.6, 1.2, 0.6. (1) 建立可以确定网格数据的线性方程组, 并用 Matlab 求解. (2) 将网格数据乘以 256, 再取整, 用 Matlab 绘制该灰度图像