线性代数实验题03-CT图像的代数重建问题
ct重建算法
CT重建算法1. 介绍计算机断层扫描（Computed Tomography, CT）是一种通过旋转式X射线扫描来获取物体内部详细结构的成像技术。
CT重建算法是将获得的一系列投影数据转化为图像的过程。
本文将介绍CT重建算法的原理、常见算法以及应用。
2. 原理CT重建算法的原理基于X射线的相对吸收特性。
当X射线通过物体时，被吸收的程度与物体的密度有关。
通过在不同角度上获得物体的吸收投影数据，可以得到物体的密度分布。
CT重建算法将这些投影数据转换为物体的二维或三维图像。
3. 常见算法3.1 过滤回投影算法（Filtered Backprojection）过滤回投影算法是最常用的CT重建算法之一。
它在重建过程中使用反投影和滤波两个步骤。
反投影（Backprojection）是将投影数据沿着投影路径反向投射到图像平面上。
滤波（Filtering）是为了抵消投影数据中带来的伪影，通常使用高通滤波器来增强边缘。
过滤回投影算法的优点是简单、快速，适用于大部分CT重建应用。
然而，它对数据质量要求较高，容易受到噪声的影响。
3.2 代数重建技术（Algebraic Reconstruction Technique，ART）代数重建技术是一种迭代重建算法。
它通过假设一个初始图像，然后通过反复调整该图像，使其产生的投影数据与实际投影数据越来越接近。
ART算法的优点是对噪声更加稳健，并且可以提供更好的图像质量。
然而，它的计算量较大，需要更长的重建时间。
3.3 迭代重建算法除了ART算法，还有其他一些迭代重建算法，如最小二乘迭代算法、最小均方偏差迭代算法等。
这些算法的思想都是通过迭代过程逐步调整图像，使其产生的投影数据与实际投影数据更接近。
迭代重建算法的优点是能够处理高噪声情况下的重建问题，并且可以提供更好的图像质量。
然而，它的计算量较大，需要更长的重建时间。
4. 应用CT重建算法在医学领域有着广泛的应用。
它可以用于诊断与鉴别诊断，如放射影像学、肿瘤检测和血管成像等。
CT图像重建
L257
L256Байду номын сангаас
L75
L74
L46
找纵坐标
L379 L378 L257 y=6.2149 L256 h i L89 …………
3.确定接受值和所穿越长度之间的比值
• 已知中心点（-9.3133， 6.2149 ）
Radon变换
• 从前面的分析知道：旋转角度大概为-60至120度
3、射线的180个方向角度
题目分析
从这几张图和数据能获得哪些信息？ 1. 我们假设水平方向为X轴，射线与x正方向的夹角代表射线角度。能否在附件2中， 找出第几列对应90度射线？ [a b]=find(AS1==max(AS1(:))) 2. 能否在附件2中，找出第几列对应着0度射线？ 3. 以90度方向为准，第1条射线到第512条射线是从左到右，还是从右到左？ 问：你知道机器大概是从多少度旋转到多少度吗？
• 刚才我们根据猜测，认为没旋转1次为1度，但是现实中，机器也 会旋转也会有差异。
习题：
• 请编程实现求解： f(x) = x^2 - 2x – 3
开始编程：10分钟
问题：
开始编程：35分钟
求极值的函数
• fminunc Find minimum of unconstrained multivariable function 求如下函数的最小值：
fun = @(x)3*x(1)^2 + 2*x(1)*x(2) + x(2)^2 - 4*x(1) + 5*x(2);
上午总结
• 1、非线性极值问题 • 2、符号运算编程 • 3、字符操作 • 4、文档阅读
2、CT系统旋转中心位置
[[a b]=find(AS1==max(AS1(:)))
线性代数在医学图像处理中的应用 案例解析
线性代数在医学图像处理中的应用案例解析线性代数在医学图像处理中的应用近年来，随着科技的不断发展，医学图像处理技术在医疗领域中扮演着越来越重要的角色。
而线性代数作为一门重要的数学分支，也被广泛应用于医学图像处理中。
本文将以多个案例来解析线性代数在医学图像处理中的应用，展示其在提高医学诊断准确性、减少操作风险等方面的积极作用。
案例一：三维重建技术在医学图像处理中，三维重建技术是一项常用的技术。
通过将多幅二维医学图像进行重建，可以得到一个三维的结构模型，从而更准确地了解患者的病情。
在这个过程中，线性代数起到了至关重要的作用。
首先，我们可以将每一幅二维医学图像视为一个二维向量，然后将这些向量构成一个矩阵。
通过对这个矩阵进行分解和运算，可以得到一个近似原始三维结构的矩阵。
然后，通过对这个近似矩阵进行优化和逼近，最终可以得到一个高精度的三维结构模型。
其次，线性代数的矩阵运算还可以用于解决三维重建中的一些实际问题。
例如，在重建过程中，可能会遇到数据缺失或者不完整的情况。
通过利用线性代数中的矩阵填补方法，可以将缺失的数据进行估计，从而得到一个更完整的三维结构模型。
案例二：图像增强和恢复在医学图像处理中，图像增强和恢复技术被广泛应用于提高图像质量和清晰度。
而线性代数提供了一种有效的数学工具来实现图像的增强和恢复。
一种常用的图像增强技术是滤波操作。
通过对图像进行滤波，可以去除图像中的噪声，并提高图像的清晰度。
在这个过程中，线性代数中的卷积运算被广泛应用。
通过将图像视为矩阵，可以利用线性代数中的卷积定理和矩阵运算，对图像进行滤波操作，从而实现图像的增强。
此外，在医学图像处理中，还常常需要对低质量的图像进行恢复。
这种情况下，线性代数中基于最小二乘法的技术被广泛应用。
通过对图像进行建模，利用线性代数中的最小二乘法，可以对低质量的图像进行修复，从而恢复其细节和清晰度。
案例三：图像分割和分类在医学图像处理中，图像的分割和分类是非常关键的步骤。
线性代数与图像处理
线性代数与图像处理线性代数是一门研究向量空间和线性映射的代数学科，它在众多领域中具有重要的应用价值。
其中之一就是在图像处理方面的应用。
本文将探讨线性代数如何在图像处理中发挥作用。
一、像素和向量在图像处理中，图像可以被看作是一个由像素组成的网格。
每个像素都有自己的属性，如颜色、亮度等。
而这些属性可以被看作是一个向量。
例如，在RGB颜色空间中，一个像素的颜色可以由一个三维向量表示。
二、图像的线性操作线性代数的一个重要概念是线性操作。
在图像处理中，我们可以对图像进行各种线性操作，如平移、旋转、缩放等。
这些操作可以通过矩阵乘法来实现。
例如，我们可以通过矩阵乘法来对图像进行平移变换。
三、图像的滤波图像滤波是图像处理中常用的一种操作。
它可以通过线性代数中的矩阵运算来实现。
滤波可以增强或者抑制图像中的某些频率分量。
常见的滤波器包括均值滤波器、高斯滤波器等。
这些滤波器可以通过卷积操作来实现，其中卷积核可以被看作是一个滤波器的参数。
四、图像的压缩图像压缩是一种将图像数据用更小的空间存储或传输的技术。
在图像处理中，压缩可以通过线性代数中的奇异值分解（SVD）来实现。
SVD可以将一个矩阵分解为三个矩阵的乘积，其中一个矩阵包含了图像的主要信息，而其他两个矩阵包含了图像的噪声或细节信息。
通过保留主要信息，我们可以实现对图像的压缩。
五、图像的重建图像重建是图像处理中的另一个重要任务。
它可以通过线性代数中的逆运算来实现。
例如，当我们对图像进行压缩后，我们可以使用逆SVD来重建原始图像。
此外，我们还可以使用其他线性代数中的技术，如线性插值、多项式拟合等来进行图像的重建。
六、图像的特征提取图像特征提取是一种将图像中的信息抽象为数值或向量的技术。
线性代数可以帮助我们对图像进行特征提取。
例如，我们可以使用主成分分析（PCA）来提取图像的主要特征。
PCA可以通过线性变换来找到一组新的变量，这些变量具有最大的方差，并且可以用来描述图像的形状、纹理等特征。
sart 代数重建算法
sart 代数重建算法代数重建算法（Algebraic Reconstruction Technique，简称ART）是一种用于三维图像重建的算法。
该算法基于代数的原理，通过迭代运算将平面投影数据转化为三维体数据。
ART算法在医学影像领域得到广泛应用，特别是在计算机断层扫描（CT）的重建过程中。
在CT扫描中，X线束通过患者的身体产生一系列平面投影图像。
这些图像包含了关于内部结构的信息。
传统的重建算法，如滤波反投影（Filtered Back Projection，简称FBP）算法，通过将投影数据反向投影到体素空间中进行重建。
然而，FBP算法在噪声和散射等方面存在一些问题，它使用直线模型来重建图像，无法涵盖复杂的非线性成像系统。
相比之下，ART算法采用了一种不同的方法。
它使用一个重建矩阵，通过迭代更新来逼近真实的三维体数据。
ART算法的核心思想是假设体素内部的吸收系数是均匀的，投影数据之间存在线性关系。
重建过程中，根据投影数据与当前体像素之间的残差，调整每个迭代步的体像素值。
ART算法的重建过程可以简化为以下几个步骤：1.初始化重建矩阵：将重建矩阵的每个像素值初始化为某个初始值，例如零或均匀分布的值。
2.计算投影矩阵预测：使用当前体像素猜测值计算投影数据的预测值。
3.计算残差：将实际投影数据与预测的投影数据之间的差异作为残差。
4.更新体像素值：根据残差和投影矩阵来更新每个迭代步的体像素值。
更新规则可以根据实际应用的需求来设计，例如最小二乘法。
5.重复迭代：重复步骤2至4，直到达到预设的收敛条件。
ART算法的优点之一是可以在重建过程中考虑非线性估计和噪声模型。
它可以通过增加迭代次数来增加重建图像的准确性，但也需要注意迭代次数过多可能导致图像模糊或产生伪影。
然而，ART算法也存在一些挑战。
首先，它对初始模型非常敏感，初始模型的选择可能会影响最终的重建结果。
其次，ART算法在处理大规模数据时，计算复杂度较高，需要较长的运算时间。
数学建模案例分析--线性代数建模案例20例
线性代数建模案例汇编目录案例一. 交通网络流量分析问题1案例二. 配方问题4案例三. 投入产出问题6案例四. 平板的稳态温度分布问题7案例五. CT图像的代数重建问题11案例六. 平衡结构的梁受力计算13案例七. 化学方程式配平问题16案例八. 互付工资问题17案例九. 平衡价格问题19案例十. 电路设计问题20案例十一. 平面图形的几何变换22案例十二. 太空探测器轨道数据问题24案例十三. 应用矩阵编制Hill密码25案例十四. 显示器色彩制式转换问题27案例十五. 人员流动问题29案例十六. 金融公司支付基金的流动31案例十七. 选举问题33案例十八. 简单的种群增长问题34案例十九. 一阶常系数线性齐次微分方程组的求解36 案例二十. 最值问题38附录数学实验报告模板错误!未定义书签。
案例一. 交通网络流量分析问题城市道路网中每条道路、每个交叉路口的车流量调查，是分析、评价及改善城市交通状况的基础。
根据实际车流量信息可以设计流量控制方案，必要时设置单行线，以免大量车辆长时间拥堵。
【模型准备】 某城市单行线如下图所示, 其中的数字表示该路段每小时按箭头方向行驶的车流量(单位: 辆).图3 某城市单行线车流量(1) 建立确定每条道路流量的线性方程组.(2) 为了唯一确定未知流量, 还需要增添哪几条道路的流量统计? (3) 当x 4 = 350时, 确定x 1, x 2, x 3的值.(4) 若x 4 = 200, 则单行线应该如何改动才合理?【模型假设】 (1) 每条道路都是单行线. (2) 每个交叉路口进入和离开的车辆数目相等.【模型建立】 根据图3和上述假设, 在①, ②, ③, ④四个路口进出车辆数目分别满足500 = x 1 + x 2① 400 + x 1 = x 4 + 300 ② x 2 + x 3 = 100 + 200 ③ x 4 = x 3 + 300 ④ 【模型求解】根据上述等式可得如下线性方程组12142334500100300300x x x x x x x x +=⎧⎪-=-⎪⎨+=⎪⎪-+=⎩其增广矩阵(A , b ) =1100500100110001103000011300⎛⎫ ⎪--⎪ ⎪ ⎪-⎝⎭−−−−→初等行变换10011000101600001130000000--⎛⎫ ⎪⎪-- ⎪⎪⎝⎭由此可得142434100600300x x x x x x -=-⎧⎪+=⎨⎪-=-⎩ 即142434100600300x x x x x x =-⎧⎪=-+⎨⎪=-⎩. 为了唯一确定未知流量, 只要增添x 4统计的值即可. 当x 4 = 350时, 确定x 1 = 250, x 2 = 250, x 3 = 50.若x 4 = 200, 则x 1 = 100, x 2 = 400, x 3 = -100 < 0. 这表明单行线“③←④”应该改为“③→④”才合理.【模型分析】(1) 由(A , b )的行最简形可见, 上述方程组中的最后一个方程是多余的. 这意味着最后一个方程中的数据“300”可以不用统计.(2) 由142434100600300x x x x x x =-⎧⎪=-+⎨⎪=-⎩可得213141500200100x x x x x x =-+⎧⎪=-⎨⎪=+⎩, 123242500300600x x x x x x =-+⎧⎪=-+⎨⎪=-+⎩, 132343200300300x x x x x x =+⎧⎪=-+⎨⎪=+⎩, 这就是说x 1, x 2, x 3, x 4这四个未知量中, 任意一个未知量的值统计出来之后都可以确定出其他三个未知量的值.Matlab 实验题某城市有下图所示的交通图, 每条道路都是单行线, 需要调查每条道路每小时的车流量. 图中的数字表示该条路段的车流数. 如果每个交叉路口进入和离开的车数相等, 整个图中进入和离开的车数相等.图4 某城市单行线车流量(1)建立确定每条道路流量的线性方程组.(2)分析哪些流量数据是多余的.(3)为了唯一确定未知流量, 需要增添哪几条道路的流量统计.案例二. 配方问题在化工、医药、日常膳食等方面都经常涉及到配方问题. 在不考虑各种成分之间可能发生某些化学反应时, 配方问题可以用向量和线性方程组来建模. 【模型准备】一种佐料由四种原料A 、B 、C 、D 混合而成. 这种佐料现有两种规格, 这两种规格的佐料中, 四种原料的比例分别为2:3:1:1和1:2:1:2. 现在需要四种原料的比例为4:7:3:5的第三种规格的佐料. 问: 第三种规格的佐料能否由前两种规格的佐料按一定比例配制而成?【模型假设】 (1) 假设四种原料混合在一起时不发生化学变化. (2) 假设四种原料的比例是按重量计算的. (3) 假设前两种规格的佐料分装成袋, 比如说第一种规格的佐料每袋净重7克(其中A 、B 、C 、D 四种原料分别为2克, 3克, 1克, 1克), 第二种规格的佐料每袋净重6克(其中A 、B 、C 、D 四种原料分别为1克, 2克, 1克, 2克). 【模型建立】 根据已知数据和上述假设, 可以进一步假设将x 袋第一种规格的佐料与y 袋第二种规格的佐料混合在一起, 得到的混合物中A 、B 、C 、D 四种原料分别为4克, 7克, 3克, 5克, 则有以下线性方程组24,327,3,2 5.x y x y x y x y +=⎧⎪+=⎨+=⎪+=⎩ 【模型求解】上述线性方程组的增广矩阵(A , b ) =214327113125⎛⎫ ⎪ ⎪ ⎪ ⎪⎝⎭−−−−→初等行变换101012000000⎛⎫ ⎪⎪ ⎪ ⎪⎝⎭,可见{1,2.x y == 又因为第一种规格的佐料每袋净重7克, 第二种规格的佐料每袋净重6克, 所以第三种规格的佐料能由前两种规格的佐料按7:12的比例配制而成. 【模型分析】(1) 若令α1 = (2, 3, 1, 1)T , α2 = (1, 2, 1, 1)T , β = (4, 7, 5, 3)T , 则原问题等价于“线性方程组Ax = b 是否有解”, 也等价于“β能否由α1, α2线性表示”.(2) 若四种原料的比例是按体积计算的, 则还要考虑混合前后体积的关系(未必是简单的叠加), 因而最好还是先根据具体情况将体积比转换为重量比, 然后再按上述方法处理.(3) 上面的模型假设中的第三个假设只是起到简化运算的作用. 如果直接设x 克第一种规格的佐料与y 克第二种规格的佐料混合得第三种规格的佐料, 则有下表因而有如下线性方程组214(),7619327(),7619113(),7619125().7619x y x y x y x y x y x y x y x y ⎧+=+⎪⎪⎪+=+⎪⎨⎪+=+⎪⎪⎪+=+⎪⎩(*) 【模型检验】把x = 7, y = 12代入上述方程组(*), 则各等式都成立. 可见模型假设中的第三个假设不影响解的正确性.Matlab 实验题蛋白质、碳水化合物和脂肪是人体每日必须的三种营养, 但过量的脂肪摄入不利于健康.人们可以通过适量的运动来消耗多余的脂肪. 设三种食物(脱脂牛奶、大豆面粉、乳清)每100克中蛋白质、碳水化合物和脂肪的含量以及慢跑5分钟消耗蛋白质、碳水化合物和脂肪的量如下表.问怎样安排饮食和运动才能实现每日的营养需求?案例三. 投入产出问题在研究多个经济部门之间的投入产出关系时, W. Leontief 提出了投入产出模型. 这为经济学研究提供了强有力的手段. W. Leontief 因此获得了1973年的Nobel 经济学奖.【模型准备】某地有一座煤矿, 一个发电厂和一条铁路. 经成本核算, 每生产价值1元钱的煤需消耗0.3元的电; 为了把这1元钱的煤运出去需花费0.2元的运费; 每生产1元的电需0.6元的煤作燃料; 为了运行电厂的辅助设备需消耗本身0.1元的电, 还需要花费0.1元的运费; 作为铁路局, 每提供1元运费的运输需消耗0.5元的煤, 辅助设备要消耗0.1元的电. 现煤矿接到外地6万元煤的订货, 电厂有10万元电的外地需求, 问: 煤矿和电厂各生产多少才能满足需求? 【模型假设】假设不考虑价格变动等其他因素.【模型建立】设煤矿, 电厂, 铁路分别产出x 元, y 元, z 元刚好满足需求. 则有下表根据需求, 应该有(0.60.5)60000(0.30.10.1)100000(0.20.1)0x y z y x y z z x y -+=⎧⎪-++=⎨⎪-+=⎩, 即0.60.5600000.30.90.11000000.20.10x y z x y z x y z --=⎧⎪-+-=⎨⎪--+=⎩ 【模型求解】在Matlab 命令窗口输入以下命令>> A = [1,-0.6,-0.5;-0.3,0.9,-0.1;-0.2,-0.1,1]; b = [60000;100000;0]; >> x = A\bMatlab 执行后得 x =1.0e+005 *1.99661.84150.5835可见煤矿要生产1.9966⨯105元的煤, 电厂要生产1.8415⨯105元的电恰好满足需求.【模型分析】令x =xyz⎛⎫⎪⎪⎝⎭, A =00.60.50.30.10.10.20.10⎛⎫⎪⎪⎝⎭, b =60000100000⎛⎫⎪⎪⎝⎭, 其中x称为总产值列向量,A称为消耗系数矩阵, b称为最终产品向量, 则Ax =00.60.50.30.10.10.20.10⎛⎫⎪⎪⎝⎭xyz⎛⎫⎪⎪⎝⎭=0.60.50.30.10.10.20.1y zx y zx y+⎛⎫⎪++⎪+⎝⎭根据需求, 应该有x-Ax = b, 即(E-A)x = b. 故x = (E-A)-1b.Matlab实验题某乡镇有甲、乙、丙三个企业. 甲企业每生产1元的产品要消耗0.25元乙企业的产品和0.25元丙企业的产品. 乙企业每生产1元的产品要消耗0.65元甲企业的产品, 0.05元自产的产品和0.05元丙企业的产品. 丙企业每生产1元的产品要消耗0.5元甲企业的产品和0.1元乙企业的产品. 在一个生产周期内, 甲、乙、丙三个企业生产的产品价值分别为100万元, 120万元, 60万元, 同时各自的固定资产折旧分别为20万元, 5万元和5万元.(1) 求一个生产周期内这三个企业扣除消耗和折旧后的新创价值.(2) 如果这三个企业接到外来订单分别为50万元, 60万元, 40万元, 那么他们各生产多少才能满足需求?案例四. 平板的稳态温度分布问题在热传导的研究中, 一个重要的问题是确定一块平板的稳态温度分布. 根据…定律, 只要测定一块矩形平板四周的温度就可以确定平板上各点的温度.图8 一块平板的温度分布图【模型准备】如图9所示的平板代表一条金属梁的截面. 已知四周8个节点处的温度(单位°C), 求中间4个点处的温度T 1, T 2, T 3, T 4.图9 一块平板的温度分布图【模型假设】假设忽略垂直于该截面方向上的热传导, 并且每个节点的温度等于与它相邻的四个节点温度的平均值.【模型建立】根据已知条件和上述假设, 有如下线性方程组1232143144231(90100)41(8060)41(8060)41(5050)4T T T T T T T T T T T T ⎧=+++⎪⎪⎪=+++⎪⎨⎪=+++⎪⎪=+++⎪⎩ 【模型求解】将上述线性方程组整理得1231241342344190414041404100T T T T T T T T T T T T --=⎧⎪-+-=⎪⎨-+-=⎪--+=⎪⎩. 在Matlab 命令窗口输入以下命令T 1T 2 T 3 T 4 10080908060506050>> A = [4,-1,-1,0;-1,4,0,-1;-1,0,4,-1;0,-1,-1,4]; b = [190;140;140;100];>> x = A\b; x’Matlab执行后得ans =82.9167 70.8333 70.8333 60.4167可见T1 = 82.9167, T2 = 70.8333, T3 = 70.8333, T4 = 60.4167.参考文献陈怀琛, 高淑萍, 杨威, 工程线性代数,: 电子工业, 2007. 页码: 15-16.Matlab实验题假定下图中的平板代表一条金属梁的截面, 并忽略垂直于该截面方向上的热传导. 已知平板内部有30个节点, 每个节点的温度近似等于与它相邻的四个节点温度的平均值. 设4条边界上的温度分别等于每位同学学号的后四位的5倍, 例如学号为16308209的同学计算本题时, 选择T l = 40, T u = 10, T r = 0, T d = 45.图10 一块平板的温度分布图(1) 建立可以确定平板内节点温度的线性方程组.(2) 用Matlab软件求解该线性方程组.(3) 用Matlab中的函数mesh绘制三维平板温度分布图.案例五. CT图像的代数重建问题X射线透视可以得到3维对象在2维平面上的投影, CT则通过不同角度的X射线得到3维对象的多个2维投影, 并以此重建对象内部的3维图像. 代数重建方法就是从这些2维投影出发, 通过求解超定线性方程组, 获得对象内部3维图像的方法.图11双层螺旋CT 图12 CT图像这里我们考虑一个更简单的模型, 从2维图像的1维投影重建原先的2维图像. 一个长方形图像可以用一个横竖均匀划分的离散网格来覆盖, 每个网格对应一个像素, 它是该网格上各点像素的均值. 这样一个图像就可以用一个矩阵表示,其元素就是图像在一点的灰度值(黑白图像). 下面我们以3⨯3图像为例来说明.3⨯3图像各点的灰度值水平方向上的叠加值x1 = 1 x2 = 0 x3 = 0 x1 + x2 + x3 = 1x4 = 0 x5 = 0.5 x6 = 0.5 x4 + x5 + x6 = 1x7 = 0.5 x8 = 0 x9 = 1 x7 + x8 + x9 = 1.5 竖直方向上的叠加值x1 + x4 + x7= 1.5x2 + x5 + x8= 0.5x3 + x6 + x9= 1.5i色. 如果我们不知道网格中的数值, 只知道沿竖直方向和水平方向的叠加值, 为了确定网格中的灰度值, 可以建立线性方程组(含有6个方程, 9个未知数)123456369111x x xx x xx x x++=⎧⎪++=⎪⎨⎪++=⎪⎩显然该方程组的解是不唯一的, 为了重建图像, 必须增加叠加值. 如我们增加从右上方到左下方的叠加值, 则方程组将增加5个方程x1 = 1,x2 + x4 = 0,x3 + x5 + x7 = 1,x 6 + x 8 = 0.5, x 9 = 1,和上面的6个方程放在一起构成一个含有11个方程, 9个未知数的线性方程组. 【模型准备】设3⨯3图像中第一行3个点的灰度值依次为x 1, x 2, x 3, 第二行3个点的灰度值依次为x 4, x 5,x 6, 第三行3个点的灰度值依次为x 7, x 8, x 9. 沿竖直方向的叠加值依次为1.5, 0.5, 1.5, 沿水平方向的叠加值依次为1, 1, 1.5, 沿右上方到左下方的叠加值依次为1, 0, 1, 0.5, 1. 确定x 1, x 2, …, x 9的值.【模型建立】由已知条件可得(含有11个方程, 9个未知数的)线性方程组1234569111x 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 可见上述方程组的解不唯一. 其中的一个特解为x 1 = 1, x 2 = 0, x 3 = 0, x 4 = 0, x 5 = 0.5, x 6 = 0.5, x 7 = 0.5, x 8 = 0, x 9 = 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 绘制该灰度图像.案例六. 平衡结构的梁受力计算在桥梁、房顶、铁塔等建筑结构中, 涉及到各种各样的梁. 对这些梁进行受力分析是设计师、工程师经常做的事情.图14 埃菲尔铁塔局部下面以双杆系统的受力分析为例, 说明如何研究梁上各铰接点处的受力情况. 【模型准备】在图15所示的双杆系统中, 已知杆1重G1 = 200牛顿, 长L1 = 2米, 与水平方向的夹角为θ1 = π/6, 杆2重G2 = 100牛顿, 长L2 = 2米, 与水平方向的夹角为θ2 = π/4. 三个铰接点A, B, C所在平面垂直于水平面. 求杆1, 杆2在铰接点处所受到的力.图15双杆系统【模型假设】假设两杆都是均匀的. 在铰接点处的受力情况如图16所示.【模型建立】对于杆1:水平方向受到的合力为零, 故N1 = N3,竖直方向受到的合力为零, 故N2 + N4 = G1,以点A为支点的合力矩为零, 故(L1sinθ1)N3 + (L1cosθ1)N4 = (12L1cosθ1)G1.图16 两杆受力情况对于杆2类似地有AC杆1杆2CN1N2N3N5N6G1G2A B杆1杆2π/6π/4N 5 = N 7, N 6 = N 8 + G 2, (L 2sin θ2)N 7 = (L 2cos θ2)N 8 + (12L 2cos θ2)G 2.此外还有N 3 = N 7, N 4 = N 8. 于是将上述8个等式联立起来得到关于N 1, N 2, …, N 8的线性方程组:132414800N N N N G N N -=⎧⎪+=⎪⎨⎪⎪-=⎩ 【模型求解】在Matlab 命令窗口输入以下命令>> G1=200; L1=2; theta1=pi/6; G2=100; L2=sqrt(2); theta2=pi/4; >> A = [1,0,-1,0,0,0,0,0;0,1,0,1,0,0,0,0;0,0,L1*sin(theta1),L1*cos(theta1),0,0,0,0;0,0,0,0,1,0,-1,0; 0,0,0,0,0,1,0,-1;0,0,0,0,0,0,L2*sin(theta2),-L2*cos(theta2); 0,0,1,0,0,0,-1,0;0,0,0,1,0,0,0,-1];>> b = [0;G1;0.5*L1*cos(theta1)*G1;0;G2;0.5*L2*cos(theta2)*G2;0;0]; >> x = A\b; x ’ Matlab 执行后得 ans =95.0962 154.9038 95.0962 45.0962 95.0962 145.0962 95.0962 45.0962【模型分析】最后的结果没有出现负值, 说明图16中假设的各个力的方向与事实一致. 如果结果中出现负值, 则说明该力的方向与假设的方向相反. 参考文献陈怀琛, 高淑萍, 杨威, 工程线性代数,: 电子工业, 2007. 页码: 157- 158.Matlab 实验题有一个平面结构如下所示, 有13条梁(图中标号的线段)和8个铰接点(图中标号的圈)联结在一起. 其中1号铰接点完全固定, 8号铰接点竖直方向固定, 并在2号, 5号和6号铰接点上, 分别有图示的10吨, 15吨和20吨的负载. 在静平衡的条件下,任何一个铰接点上水平和竖直方向受力都是平衡的. 已知每条斜梁的角度都是45º.(1) 列出由各铰接点处受力平衡方程构成的线性方程组. (2) 用Matlab 软件求解该线性方程组, 确定每条梁受力情况.图17 一个平面结构的梁案例七. 化学方程式配平问题在用化学方法处理污水过程中, 有时会涉及到复杂的化学反应. 这些反应的化学方程式是分析计算和工艺设计的重要依据. 在定性地检测出反应物和生成物之后,可以通过求解线性方程组配平化学方程式.【模型准备】某厂废水中含K, 其浓度为650mg/L. 现用氯氧化法处理, 发生如下反应:K + 2KOH + Cl 2 = KO+ 2KCl + H 2O.投入过量液氯, 可将氰酸盐进一步氧化为氮气. 请配平下列化学方程式:KO +KOH +Cl 2 ===CO 2+N 2+KCl +H 2O.(注: 题目摘自XX 省XX 外国语学校2008-2009学年高三第三次月考化学试卷) 【模型建立】设x 1KO ＋x 2KOH ＋x 3Cl 2 === x 4CO 2 ＋x 5N 2 ＋x 6KCl ＋x 7H 2O,则1261247141527362222x x x x x x xx x x x x x x x +=⎧⎪+=+⎪⎪=⎪⎨=⎪⎪=⎪=⎪⎩, 即1261247141527360200202020x x x x x x x x x x x x x x x +-=⎧⎪+--=⎪⎪-=⎪⎨-=⎪⎪-=⎪-=⎪⎩ 【模型求解】在Matlab 命令窗口输入以下命令>> A = [1,1,0,0,0,-1,0;1,1,0,-2,0,0,-1;1,0,0,-1,0,0,0;1,0,0,0,-2,0,0;0,1,0,0,0,0,-2;0,0,2,0,0,-1,0];>> x = null(A,’r ’); format rat, x ’Matlab 执行后得 ans =1 2 3/2 1 1/2 3 1 可见上述齐次线性方程组的通解为x = k (1, 2, 3/2, 1, 1/2, 3, 1)T .取k = 2得x = (2, 4, 3, 2, 1, 6, 2)T . 可见配平后的化学方程式如下2KO + 4KOH + 3Cl 2 ===2CO 2+ N 2+ 6KCl + 2H 2O.【模型分析】利用线性方程组配平化学方程式是一种待定系数法. 关键是根据化学方程式两边所涉及到的各种元素的量相等的原则列出方程. 所得到的齐次线性方程组Ax = θ中所含方程的个数等于化学方程式中元素的种数s , 未知数的个数就是化学方程式中的项数n .当r(A ) = n -1时, Ax = θ的基础解系中含有1个(线性无关的)解向量. 这时在通解中取常数k 为各分量分母的最小公倍数即可. 例如本例中1, 2, 3/2, 1, 1/2, 3, 1分母的最小公倍数为2, 故取k = 2.当r(A ) ≤n -2时, Ax = θ的基础解系中含有2个以上的线性无关的解向量. 这时可以根据化学方程式中元素的化合价的上升与下降的情况, 在原线性方程组中添加新的方程. Matlab 实验题配平下列反应式(1) FeS + KMnO 4 + H 2SO 4—— K 2SO 4 + MnSO 4 + Fe 2(SO 4)3 + H 2O + S ↓ (2) Al 2(SO 4)3 + Na 2CO 3 + H 2O —— Al(OH)3↓+ CO 2↑+ Na 2SO 4案例八. 互付工资问题互付工资问题是多方合作相互提供劳动过程中产生的. 比如农忙季节, 多户农民组成互助组, 共同完成各户的耕、种、收等农活. 又如木工, 电工, 油漆工等组成互助组, 共同完成各家的装潢工作. 由于不同工种的劳动量有所不同, 为了均衡各方的利益, 就要计算互付工资的标准.【模型准备】现有一个木工, 电工, 油漆工. 相互装修他们的房子, 他们有如下协议:(1) 每人工作10天(包括在自己家的日子), (2) 每人的日工资一般的市价在60~80元之间, (3) 日工资数应使每人的总收入和总支出相等.求每人的日工资. 【模型假设】假设每人每天工作时间长度相同. 无论谁在谁家干活都按正常情况工作, 既不偷懒, 也不加班.【模型建立】设木工, 电工, 油漆工的日工资分别为x , y , z 元, 则由下表可得2610451044310x y z xx y z y x y z z++=⎧⎪++=⎨⎪++=⎩, 即8604504470x y z x y z x y z -++=⎧⎪-+=⎨⎪+-=⎩【模型求解】在Matlab 命令窗口输入以下命令>> A = [-8,1,6;4,-5,1;4,4,-7];>> x = null(A,’r ’); format rat, x ’ Matlab 执行后得ans =31/36 8/9 1可见上述齐次线性方程组的通解为x = k (31/36, 8/9, 1)T . 因而根据“每人的日工资一般的市价在60~80元之间”可知60 ≤3631k <98k < k ≤ 80, 即 312160≤k ≤ 80.也就是说, 木工, 电工, 油漆工的日工资分别为3631k 元, 98k 元, k 元, 其中312160≤k ≤ 80. 为了简便起见, 可取k = 72, 于是木工, 电工, 油漆工的日工资分别为62元, 64元, 72元.【模型分析】事实上各人都不必付自己工资, 这时各家应付工资和各人应得收入如下6845447y z x x z y x y z +=⎧⎪+=⎨⎪+=⎩, 即8604504470x y z x y z x y z -++=⎧⎪-+=⎨⎪+-=⎩ 可见这样得到的方程组与前面得到的方程组是一样的.Matlab 实验题甲, 乙, 丙三个农民组成互助组, 每人工作6天(包括为自己家干活的天数), 刚好完成他们三人家的农活, 其中甲在甲, 乙, 丙三家干活的天数依次为: 2, 2.5, 1.5; 乙在甲, 乙, 丙三家各干2天活, 丙在甲, 乙, 丙三家干活的天数依次为: 1.5, 2, 2.5. 根据三人干活的种类, 速度和时间, 他们确定三人不必相互支付工资刚好公平. 随后三人又合作到邻村帮忙干了2天(各人干活的种类和强度不变), 共获得工资500元.问他们应该怎样分配这500元工资才合理?案例九. 平衡价格问题为了协调多个相互依存的行业的平衡发展, 有关部门需要根据每个行业的产出在各个行业中的分配情况确定每个行业产品的指导价格, 使得每个行业的投入与产出都大致相等.【模型准备】假设一个经济系统由煤炭、电力、钢铁行业组成, 每个行业的产出在各个行业中的分配如下表所示:等的平衡价格.【模型假设】假设不考虑这个系统与外界的联系.【模型建立】把煤炭、电力、钢铁行业每年总产出的价格分别用x 1,x 2, x 3表示, 则123212331230.40.60.60.10.20.40.50.2x x x x x x x x x x x =+⎧⎪=++⎨⎪=++⎩, 即1231231230.40.600.60.90.200.40.50.80x x x x x x x x x --=⎧⎪-+-=⎨⎪--+=⎩. 【模型求解】在Matlab 命令窗口输入以下命令>> A = [1,-0.4,-0.6;-0.6,0.9,-0.2;-0.4,-0.5,0.8]; >> x = null(A,’r ’); format short, x ’ Matlab 执行后得ans =0.9394 0.8485 1.0000 可见上述齐次线性方程组的通解为x = k(0.9394, 0.8485, 1)T.这就是说, 如果煤炭、电力、钢铁行业每年总产出的价格分别0.9394亿元, 0.8485亿元, 1亿元, 那么每个行业的投入与产出都相等.【模型分析】实际上, 一个比较完整的经济系统不可能只涉及三个行业, 因此需要统计更多的行业间的分配数据.Matlab实验题假设一个经济系统由煤炭、石油、电力、钢铁、机械制造、运输行业组成, 每个行业的产出在各个行业中的分配如下表所示:产出分配购买者煤炭石油电力钢铁制造运输0 0 0.2 0.1 0.2 0.2 煤炭0 0 0.1 0.1 0.2 0.1 石油0.5 0.1 0.1 0.2 0.1 0.1 电力0.4 0.1 0.2 0 0.1 0.4 钢铁0 0.1 0.3 0.6 0 0.2 制造0.1 0.7 0.1 0 0.4 0 运输等的平衡价格.案例十. 电路设计问题电路是电子元件的神经系统. 参数的计算是电路设计的重要环节. 其依据来自两个方面: 一是客观需要, 二是物理学定律.图22 USB扩展板【模型准备】假设图23中的方框代表某类具有输入和输出终端的电路. 用11vi⎛⎫⎪⎝⎭记录输入电压和输入电流(电压v以伏特为单位, 电流i以安培为单位), 用22vi⎛⎫⎪⎝⎭记录输出电压和输入电流. 若22vi⎛⎫⎪⎝⎭= A11vi⎛⎫⎪⎝⎭,则称矩阵A为转移矩阵.图23 具有输入和输出终端的电子电路图图24给出了一个梯形网络, 左边的电路称为串联电路, 电阻为R 1(单位: 欧姆). 右边的电路是并联电路, 电路R 2. 利用欧姆定理和楚列斯基定律, 我们可以得到串联电路和并联电路的转移矩阵分别是1101R -⎛⎫ ⎪⎝⎭和2101/1R ⎛⎫ ⎪-⎝⎭串联电路 并联电路图24 梯形网络设计一个梯形网络, 其转移矩阵是180.55-⎛⎫⎪-⎝⎭. 【模型假设】假设导线的电阻为零.【模型建立】设A 1和A 2分别是串联电路和并联电路的转移矩阵, 则输入向量x 先变换成A 1x , 再变换到A 2(A 1x ). 其中A 2A 1 =2101/1R ⎛⎫ ⎪-⎝⎭1101R -⎛⎫ ⎪⎝⎭=121211/1/R R R R -⎛⎫ ⎪-+⎝⎭就是图22中梯形网络的转移矩阵.于是, 原问题转化为求R 1, R 2的值使得121211/1/R R R R -⎛⎫ ⎪-+⎝⎭=180.55-⎛⎫ ⎪-⎝⎭. 【模型求解】由121211/1/R R R R -⎛⎫ ⎪-+⎝⎭=180.55-⎛⎫ ⎪-⎝⎭可得121281/0.51/5R R R R -=-⎧⎪-=-⎨⎪+=⎩. 根据其中的前两个方程可得R 1 = 8, R 2 = 2. 把R 1 = 8, R 2 = 2代入上面的第三个方程确实能使等式成立. 这就是说在图22中梯形网络中取R 1 = 8, R 2 = 2即为所求.【模型分析】若要求的转移矩阵改为180.54-⎛⎫⎪-⎝⎭, 则上面的梯形网络无法实现. 因为v 2这时对应的方程组是121281/0.51/4R R R R -=-⎧⎪-=-⎨⎪+=⎩. 根据前两个方程依然得到R 1 = 8, R 2 = 2, 但把R 1= 8, R 2 = 2代入上第三个方程却不能使等式成立.练习题根据基尔霍夫回路电路定律(各节点处流入和流出的电流强度的代数和为零, 各回路中各支路的电压降之和为零), 列出下图所示电路中电流i 1, i 2, i 3所满足的线性方程组, 并用矩阵形式表示:图25简单的回路案例十一. 平面图形的几何变换随着计算机科学技术的发展, 计算机图形学的应用领域越来越广, 如仿真设计、效果图制作、动画片制作、电子游戏开发等.图形的几何变换, 包括图形的平移、旋转、放缩等, 是计算机图形学中经常遇到的问题. 这里暂时只讨论平面图形的几何变换.【模型准备】平面图形的旋转和放缩都很容易用矩阵乘法实现, 但是图形的平移并不是线性运算, 不能直接用矩阵乘法表示. 现在要求用一种方法使平移、旋转、放缩能统一用矩阵乘法来实现. 【模型假设】设平移变换为(x , y ) → (x +a , y +b )旋转变换(绕原点逆时针旋转θ角度)为(x , y ) → (x cos θ-y sin θ, x sin θ + y cos θ)放缩变换(沿x 轴方向放大s 倍, 沿y 轴方向放大t 倍)为(x , y ) → (sx , ty )【模型求解】R 2中的每个点(x , y )可以对应于R 3中的(x , y , 1). 它在xOy 平面上方1单E 12位的平面上. 我们称(x , y , 1)是(x , y )的齐次坐标. 在齐次坐标下, 平移变换(x , y ) → (x +a , y +b )可以用齐次坐标写成(x , y , 1) → (x +a , y +b , 1).于是可以用矩阵乘积1001001a b ⎛⎫ ⎪ ⎪⎝⎭1x y ⎛⎫ ⎪ ⎪⎝⎭=1x a y b +⎛⎫⎪+ ⎪⎝⎭实现.旋转变换(x , y ) → (x cos θ-y sin θ, x sin θ + y cos θ)可以用齐次坐标写成(x , y , 1) → (x cos θ-y sin θ, x sin θ + y cos θ, 1). 于是可以用矩阵乘积cos sin 0sin cos 0001θθθθ-⎛⎫ ⎪ ⎪⎝⎭1x y ⎛⎫ ⎪ ⎪⎝⎭=cos sin sin cos 1x y x y θθθθ-⎛⎫⎪+ ⎪⎝⎭实现.放缩变换(x , y ) → (sx , ty )可以用齐次坐标写成(x , y , 1) → (sx , ty , 1).于是可以用矩阵乘积0000001s t ⎛⎫ ⎪ ⎪⎝⎭1x y ⎛⎫ ⎪ ⎪⎝⎭=1sx ty ⎛⎫⎪ ⎪⎝⎭实现.【模型分析】由上述求解可以看出, R 2中的任何线性变换都可以用分块矩阵1⎛⎫⎪⎝⎭A O O 乘以齐次坐标实现, 其中A 是2阶方阵. 这样, 只要把平面图形上点的齐次坐标写成列向量, 平面图形的每一次几何变换, 都可通过左乘一个3阶变换矩阵来实现.参考文献David C. Lay, 线性代数及其应用, 沈复兴, 傅莺莺等译,: 人民邮电, 2009. 页码: 139-141.Matlab 实验题在Matlab 命令窗口输入以下命令 >>clear all , clc,>>t=[1,3,5,11,13,15]*pi/8; >>x=sin(t); y=cos(t); >>fill(x,y,'r'); >>grid on ;>>axis([-2.4, 2.4, -2, 2])运行后得图25.图26Matlab绘制的图形(1) 写出该图形每个顶点的齐次坐标;; 最后进行横(2) 编写Matlab程序, 先将上面图形放大0.9倍; 再逆时针旋转3坐标加0.8, 纵坐标减1的图形平移. 分别绘制上述变换后的图形.案例十二. 太空探测器轨道数据问题太空航天探测器发射以后, 可能需要调整以使探测器处在精确计算的轨道里. 雷达监测到一组列向量x1, …, x k,它们给出了不同时刻探测器的实际位置与预定轨道之间的偏差的信息.图28 火星探测器【模型准备】令X k = [x1, …, x k]. 在雷达进行数据分析时需要计算出矩阵G k = X k X k T. 一旦接收到数据向量x k+1,必须计算出新矩阵G k+1. 因为数据向量到达的速度非常快, 随着k的增加, 直接计算的负担会越来越重. 现需要给出一个算法, 使得计算G k的负担不会因为k的增加而加重.【模型求解】因为G k = X k X k T=[x 1, …, x k ]T 1T k⎡⎤⎢⎥⎢⎥⎢⎥⎣⎦x x =T 1k i i i =∑x x ,G k +1 = X k +1T1k +X =[X k , x k +1]T T 1k k +⎡⎤⎢⎥⎣⎦X x = X k X k T +x k +1T 1k +x =G k +x k +1T 1k +x ,所以一旦接收到数据向量x k +1, 只要计算x k +1T1k +x , 然后把它与上一步计算得到的G k相加即可. 这样计算G k 的负担不会因为k 的增加而加重.【模型分析】计算机计算加法的时间与计算乘法的时间相比可以忽略不计. 因此在考虑计算矩阵乘积的负担时, 只要考察乘法的次数就可以了. 设x k 的维数是n , 则X k = [x 1, …, x k ]是n ⨯k 的矩阵, G k = X k X k T 是n ⨯n 的矩阵. 直接计算G k = X k X k T 需要做n 2k 次乘法. 因而计算的负担会随着k 的增加而增加. 但是对于每一个k , 计算x k Tk x 始终只要做n 2次乘法.Matlab 实验题用Matlab 编写一个程序用于处理这个问题.案例十三. 应用矩阵编制Hill 密码密码学在经济和军事方面起着极其重要的作用. 现代密码学涉及很多高深的数学知识. 这里无法展开介绍.图29 XX 通信的基本模型密码学中将信息代码称为密码, 尚未转换成密码的文字信息称为明文, 由密码表示的信息称为密文. 从明文到密文的过程称为加密, 反之为解密. 1929年, 希尔(Hill)通过线性变换对待传输信息进行加密处理, 提出了在密码史上有重要地位的希尔加密算法. 下面我们略去一些实际应用中的细节, 只介绍最基本的思想.【模型准备】若要发出信息action, 现需要利用矩阵乘法给出加密方法和加密后得到的密文, 并给出相应的解密方法.。
CT图像重建技术
CT图像重建技术由于csdn贴图不⽅便，并且不能上传附件，我把原⽂上传到了资源空间1.引⾔计算机层析成像(Computed Tomography，CT)是通过对物体进⾏不同⾓度的射线投影测量⽽获取物体横截⾯信息的成像技术，涉及到放射物理学、数学、计算机学、图形图像学和机械学等多个学科领域。
CT技术不但给诊断医学带来⾰命性的影响．还成功地应⽤于⽆损检测、产品反求和材料组织分析等⼯业领域。
CT技术的核⼼是由投影重建图像的理论，其实质是由扫描所得到的投影数据反求出成像平⾯上每个点的衰减系数值。
图像重建的算法有很多，本⽂根据CT扫描机的发展对不同时期CT所采⽤重建算法分别进⾏介绍。
2. 平⾏束投影重建算法第⼀代和第⼆代CT机获取⼀个单独投影的采样数据是从⼀组平⾏射线获取的，这种采样类型叫平⾏投影。
平⾏投影重建算法⼀般分为直接法与间接法两⼤类。
直接法是直接计算线性⽅程系数的⽅法，如矩阵法、迭代法等。
间接法是先计算投影的傅⽴叶变换，再导出吸收系数的⽅法，如反投影法、⼆维傅⽴叶重建[1]。
法和滤波反投影法等2.1 直接法2.1.1 矩阵法设⼀个物体的内部吸收系数矩阵为：(1)为了求得该矩阵中的元素值，我们可以先计算该矩阵在T个⾓度下的T组投影值,如设⽔平⽅向时,则：(2)同样其它⾓度下也有类似⽅程，把所有⽅程联⽴得到求解，即可求得所有u值。
通常情况下，由于联⽴⽅程组的数⽬往往不同于未知数个数，且可能有不少重复的⽅程，这样形成的不是⽅阵，所以⼀般不满秩，此时需要利⽤⼴义逆矩阵法进⾏求解。
2.1.2 迭代法实际应⽤中，由于图像尺⼨较⼤，联⽴的⽅程个数较多，采⽤直接采⽤解析法难度较⼤，因此提出了迭代重建⽅法。
迭代法的主要思想是：从⼀个假设的初[2]。
始图像出发，采⽤迭代的⽅法，将根据⼈为设定并经理论计算得到的投影值同实验测得的投影值⽐较，不断进⾏逼近，按照某种最优化准则寻找最优解[2]：通常有两种迭代公式，⼀种是加法迭代公式（3）[2]：另⼀种是乘法迭代公式(4)两式中是相邻两次迭代的结果；是某⼀⾓度的实测投影值，是计算过程的计算投影值,是投影的某⼀射线穿过点的点数，即计算投影值的射线所经过的像素的数⽬，是松弛因⼦。
第3节X-CT图像重建方法
反投影法 --------即把各个方向得到的投影值沿反方向反 馈到相应方向的体素中 ,并把各个投影值迭 加, 然后把每个体素的μ值减去底数,再除以 它们的最大公约数,使各体素μ值间为简单的 比例.最后得到各体素的μ值. 例: 2х2矩阵 反投影法的缺点: 星状伪影 滤波反投影法: 将投影信号进行滤波处理,使信号两侧各 加一个正\负脉冲,然后再进行反投影计算.
三.反投影法
第四节 X-CT扫描机简介
一.X-CT扫描机的主要组成部分 扫描机架: 内有X射线管和检测器阵(扫描系统) 病人床及控制系统: 步进电机等 电源 监视器 二.CT扫描方式 1. 单束扫描 (第一代) 2.窄角扇束扫描 (第二代) 3. 广角扇束扫描 (第三代)-----多层螺旋扫描* 4. 反扇束扫描 (第四代) 5.动态空间重建技术(第五代) 6.电子束扫描 (第六代)
第三节 X-CT图像重建方法
一.联立方程法 二.逐步近似法(迭代法) 三.反投影法 四.二维傅立叶变换法
一.联立方程法
•
• • • •
即根据Σμi=lnI。/I ,列方程组,求解各体 素的μ值. 其中μi:物质的线性吸收系数 lnI。/I---- X射线沿某方向的投影值 ** 独立方程的个数大于体素的个数 例: 一个2х2 矩阵 ,水平方向两个投影值分 别为5和15,垂直方向的两个投影值分别为9 和11,左对角线方向三个投影值分别为4,8,8, 右对角线方向三个投影值分别为 1,12,7, 求 每个像素的线性吸收系数.
三.X-CT机的基本结构
• 一.数据采集系统 • 1.扫描装置 • 2.X射线发生器 • 3. 检测系统 • 二. 计算机图象处理系统 • 见下节
CT原理部分 CT图像重建 CT图像重建
640×640矩阵。
（三）投影
• 投照受检体后出射的X线束强度 I 称投影 (projection，P) ，投影值的分布称为投影 函数。
• 1．近似单能X线束获取 • 使发射的X线束中主要是标识辐射的X线；
•基本原理：是将所测得的投影值按其原路 径平均的分配到每一点上，各个方向上投影 值反投影后，在影像处进行叠加，推断出原 图像。
•缺点：影像边缘处不清晰。
CT机装置
• 探测器：碘化钠、锗酸铋—高压氙气电离 室—稀土陶瓷探测器—平板探测器(非晶 硅、非晶硒）
（二）傅里叶变换法
• 是基于使图像矩阵的求解与图像投影的 傅里叶变换间建立确定的关系；或为修 正反投影法中模糊因子，从频域上校正 图像模糊部分的图像重建方法。
★(一) 反投影法
• 缺点：影像边缘处不清晰。 • 如果在一均匀的组织密度内，存在吸收系
数极不均匀的部分时，反投影图像与原图 像会出现伪影(image artifact)。 • 如圆柱形单密度体，利用反投影法所重建 图像的结果呈现出星形伪影。 • 反投影数量愈多，重建图像愈接近于原图 像，但由于存在星形伪影，而使得重建图 像的边缘部分模糊不清。
地滤波，达到满意的重建图像效果。 卷积函数h（t）选取是卷积计算的关键， h（t）称为卷积核(convo1ution kernet)。
（四）卷积反投影法
• 与滤波反投影法相比，卷积反投影法避 免了FT运算。
（四）卷积反投影法
卷积的滤波作用
（四）卷积反投影法
x cos y sin R
CT图象的重建
ai1x1 ai2x2 aiNxN bi
aij 选取取决于Ｘ射线和体素的关系
❖ （１）中心法
1 i束X射线穿过j个体素的中心
aij 0
i束X射线不穿过j个体素的中心
❖ （２）中心线法
aij
i束X射线的中心线位于 j个体素内的长度 j个体素内的长度
方程组：Y AX
其中： Y ( y1, y2 yJ )T 为投影数据；
X (x1, x2 xJ )T 为待求图像向量。
级数展开法就是由给定的投影数据Y，估计 出图像向量X； 要求所估X与一组Ｊ个基图象 ｛b1,b2,…,bJ} 的组合能足够逼近所希望重建的 图象 (x, y) 。
射线投影可由Radon变换:
X min 2 X 2 X T X
求使得
min X T X , st AX b
l
❖ （３）面积法
Ｍ束Ｘ射线，通过Ｎ个体素 Ｎ个未知数，Ｍ个方程线性方程组
N
aijxj bi i 1, , M
j 1
１、Ｍ＝Ｎ的解法有解
（１）雅可比迭代 （２）高斯－塞德尔迭代 （３）超松弛迭代
２、Ｍ和Ｎ不相等无解的解法
最小二乘法 ＡＸ＝Ｂ
＝＝》 ＡＴＡＸ＝Ａｂ
３、有无穷多解的解法
定理：存在 X min ， 使得
CT的基本原理
X射线穿过物质时其强度按指数衰减
穿过强度 Iout
， 入射强度 Iin
Iout Iin el
l 传播的距离，
衰减系数
X射线穿过一组物质时其强度变化
Iout Iin e1l enl 1 n 各物质的衰减强度
CT的一般描述
被测物体对X射线的线性衰减系数为： (x, y)
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
实验题：CT 图像的代数重建问题（线性方程组的应用）
X 射线透视可以得到3维对象在2维平面上的投影，CT 则通过不同角度的X 射线得到3维对象的多个2维投影，并以此重建对象内部的3维图像。
代数重建方法就是从这些2维投影出发，通过求解超定线性方程组，获得对象内部3维图像的方法。
这里我们考虑一个更简单的模型，从2维图像的1维投影重建原先的2维图像。
一个长方形图像可以用一个横竖均匀划分的离散网格来覆盖，每个网格对应一个像素，它是该网格上各点像素的均值。
这样一个图像就可以用一个矩阵表示，其元素就是图像在一点的灰度值（黑白图像）。
下面我们以33⨯图像为例来说明。
每个网格中的数字代表其灰度值，范围在]1,0[内，记0表示白，1表示黑，0.5为中间的灰色，沿某个方向的投影就将该方向上的灰度值相加（见上图所示）。
如果我们不知道网格中的数值，只知道沿竖直方向和水平方向的投影。
设网格按第1列、第2列、第3列的顺序排列，为了确定网格中的灰度值，可以建立线性方程组：
⎪⎪⎪⎪⎪⎪⎪⎪⎭
⎫ ⎝⎛=⎪⎪⎪⎪⎪⎪⎪⎪⎭⎫ ⎝⎛5.1115.15.05.1100100100010010010001001001111000000000111000
000000111X 显然该方程组的解是不唯一的，为了重建图像，必须增加投影数量。
如我们增加从右上到左下的投影，则方程组将增加5个方程，成为超定方程组。
考虑到测量误差，可以将超定方程组的近似解作为重建的图像数据。
1 1 1.5 1
1
问题：给定一个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绘制该灰度图像。