CN102782701A - 从x射线图像中的标志位置提取患者运动矢量 - Google Patents

Abstract

一种用于确定图像中的患者运动的方法和系统。该方法包括基于在扫描期间由扫描器产生的图像数据而获得图像。所述图像包含假定为具有刚性或半刚性构造的至少三个标志。所述至少三个标志中的各个具有在第一维度和第二维度上在扫描器的检测器面板上的实际测量的位置。该方法进一步包括为所述至少三个标志中的各个确定基准三维位置，并定义用于描述在基准三维位置和实际测量的位置之间的关系、扫描器的几何参数、以及患者运动的方程式。该方法最后包括：在数字上对这些方程式进行求解，以为对所述至少三个标志中的各个的基准三维位置和所述至少三个标志中的相应各个的实际测量的位置之间之差加以考虑的图像，得出用于描述患者运动的六分量运动矢量。

Description

从X射线图像中的标志位置提取患者运动矢量

相关申请

本专利申请是于2009年11月25日提交的、标题为“X射线图像中的标志识别和处理（MARKER IDENTIFICATION ANDPROCESSING IN X-RAY IMAGES）”的共同未决美国专利申请No.12/626,197的部分继续申请，以上申请的全部内容通过援引都加入到本申请中。本专利申请还涉及标题为“使用从X射线图像中的标志位置所提取的患者运动矢量来修正并重建X射线图像（CORRECTING ANDRECONSTRUCTING X-RAY IMAGES USING PATIENT MOTIONVECTORS EXTRACTED FROM MARKER PO SITIONS IN X-RAYIMAGES）”的申请No.12/700,032、标题为“用于X射线图像上的标志的准确子像素定位的方法（METHOD FOR ACCURATE SUB-PIXELLOCALIZATION OF MARKERS ON X-RAY IMAGES）”的申请No.12/626,218、标题为“在存在移动时用于在3D空间中的X射线标志定位的方法（METHOD FOR X-RAY MARKER LOCALIZATION IN 3DSPACE IN THE PRESENCE OF MOTION）”的申请No.12/626,268、以及标题为“用于在CT投影图像序列中跟踪X射线标志的方法（METHOD FOR TRACKING X-RAY MARKERS IN SERIAL CTPROJECTION IMAGES）”的申请No.12/626,310。

版权通告

本专利文献的一部分公开内容包含受到版权保护的材料。在其出现在专利商标局专利文件或记录时，版权拥有者不反对专利文献或专利公开内容中的任何一者被复制再现，但在其他无论任何方面都保留所有版权。

技术领域

本发明的实施例涉及图像处理。具体而言，本发明的实施例提供用于从X射线图像中的标志位置提取患者运动矢量，并使用这些向量针对患者的和其他的计划之外的运动来修正重建后的图像的方法和系统。

背景技术

现代成像系统在单个过程或扫描期间获取许多图像。使用已知为重建的过程来组合图像，以制作表示穿过患者的有限宽度“切片”的一个或多个图像。因为在为医疗或牙科过程做准备时或者为了适当地校准成像设备，经常要对所得图像加以分析，所以重要的是图像尽可能地清楚。然而，扫描期间患者移动会导致图像中的误差或模糊。

通常存在两种控制患者运动的影响的方法。一种方法是物理地限制患者。然而，物理限制往往不足以防止患者运动，这是因为患者不能保持静止（例如，年幼的孩子或由于物理极限）或是因为操作人员不适当地使用患者限制系统。另外，即使仔细地进行患者定位和限制，也会发生一些患者运动。

另一种用于控制患者运动的影响的方法是在患者运动发生时对其进行测量，并且针对其进行修正。现有的用于检测患者移动的机制包括基于激光（光学的）和超声（声音）的范围查找器（finder）或一组立体视镜（stereoscopic）相机以检测患者位置变化。然而，这些系统中所使用的检测器、范围查找器、或相机增加了系统的复杂性和成本。例如，大多数扫描环境中的旋转台架使得难以将检测器、范围查找器、或相机放置到在台架旋转期间不会在某些时刻发生阻碍的适当位置上。另外，将硬件安装在台架身上增加了检测过程的复杂性。

在计算机断层术（“CT”）成像时，将原始（raw）图像数据重建成最终图像依靠特定物理结构以数学上可描述的方式逐个投影图像地改变其定位这一事实。用于描述贯穿投影图像序列的、结构的U和V位置的图像到图像（image-to-image）轨迹的方程式具有用于描述台架的几何形状的参数、以及表示结构在三维（“3D”）空间中的物理定位的参数。如果参数被固定，则方程式准确地描述轨迹。然而，如果参数在扫描过程内变化（例如，由于患者运动），则实际轨迹不同于预期轨迹。结构的所测量的位置和结构的预期位置之间的数学关系是提取表示患者运动的运动矢量的基础。

在检测运动之后，存在多种方式以针对所检测的运动修正图像。一种方式是移动X射线源检测器以抵消运动的影响。这是相机中用于图像稳定化的主要方法。另一种方法是使用关于所检测的运动的信息在获取了图像数据之后对其进行修正。这可以包括删除并重新获取由于检测到的移动而被破坏的图像。然而，这种方法无法处理扫描期间连续的患者运动（即，破坏了图像的所有部分或者绝大部分的运动），并且增加了获取图像所需要的时间量。另外，图像数据的重新获取在一些形式的单个旋转扫描程序中是不可行的。

另一种在图像中针对所检测的运动进行修正的方法是以试图补偿所检测的运动并使图像显现为看似没有运动发生的方式平移、缩放和旋转图像。这种方法可以有效地处理患者的平移运动（translationalmotion）但其不能完全处理患者的旋转运动。例如，患者上的两个点可以是相异的，但在患者旋转之后，这两点在所得图像上会显现为一个点。图像平移或缩放不能针对这种情形进行修正，至少通常不能。

又一种修正所检测的运动的方法是如美国专利申请No.2004/0136590中所描述的那样使得重建格栅（grid）翘曲（warp）。使得重建格栅翘曲可以处理可变形的运动，例如在心脏运动期间发生的运动。因此，这种方法使用关于所检测的运动的信息修改重建格栅，以使其将假定的变形纳入考虑。尽管这种方法可以处理宽广范围的运动，但是形成格栅并将图像数据匹配到格栅是复杂的。因此，这种方法会是耗时的，并且受到各种类型的伪影（artifact）的影响。

发明内容

本发明的实施例提供用于确定用于每个投影图像的六分量患者位置（运动）矢量的系统和方法。一种方法包括，如果没有患者运动发生，则获取每个投影图像上的三个或更多个X射线基准（fiducial）标志中各个的U和V位置，并确定用以解释标志的实际定位和标志的预期定位之间的变化所需要的患者位置偏移。基于所述变化，该方法输出描述投影图像中患者的平移和旋转运动的六分量矢量。该方法允许以大约0.2毫米或更好地精确度在整个扫描期间连续跟踪患者运动。

本发明的实施例还提供用于将关于所检测的患者运动的信息结合到重建过程中的方法和系统，其使用该信息来补偿运动。一种方法包括获取关于各个投影图像的六分量位置误差向量（即，患者运动矢量）并将每个矢量与静态台架校准信息组合，以生成关于各个投影图像的投影矩阵。所述投影矩阵描述每个X射线是如何从X射线源穿过三维对象而行进到二维X射线检测器的。所述重建过程使用投影图像和所述投影矩阵来重建图像并对该图像针对所检测的患者运动进行修正。

在一个特定实施例中，提供一种用于确定图像中的患者运动的方法。该方法包括，基于在扫描期间由扫描器所产生的图像数据获得图像。所述图像包括被假定为具有刚性或半刚性构造的至少三个标志，其中，所述至少三个标志中的各个具有在第一维度和第二维度上在扫描器的检测器面板上的实际测量的位置。然后由电子处理单元确定关于所述至少三个标志中各个的基准三维位置。电子处理单元然后定义或设置用于描述基准三维位置和实际测量的位置之间的关系、扫描器的几何参数、以及患者运动的方程式。最后，电子处理单元在数字上对这些方程式进行求解，从而为所述图像得出描述患者运动的六分量运动矢量，所述六分量运动矢量对所述至少三个标志中的各个的基准三维位置和所述至少三个标志中的相应各个的实际测量的位置之间的差加以考虑。

在另一特定实施例中，基于所检测的患者运动修正并重建多个投影图像。该方法包括，在计算机处获得由扫描器产生的多个投影图像。所述多个投影图像中的各个包括至少三个标志，并且所述至少三个标志中的各个具有在第一维度和第二维度上的在扫描器的检测器面板上的所测量的位置、以及所测量的三维位置。基于所述多个投影图像中的每个中的至少三个标志，由电子处理单元确定关于所述多个投影图像中的各个的位置误差向量。所述位置误差向量定义投影图像中的患者运动。将关于所述多个投影图像中的各个的各个位置误差向量与和扫描器相关联的几何参数组合，以得出用于所述多个投影图像中的各个的投影矩阵。将所述多个投影图像以及用于所述多个投影图像中的各个的投影矩阵提供或供给到图像重建过程，该过程生成针对患者运动修正后的重建的图像。

本发明的其他方面将通过对具体说明书以及附图加以考虑而变得明显。

附图说明

图1是示出了根据一个实施例的锥形束计算机断层摄影系统。

图2示意性示出了图1的锥形束计算机断层摄影系统的部件及其几何关系。

图3示意性示出了根据一个实施例的图1的计算机。

图4示出了包含标志的图像。

图5是示出了根据一个实施例的向重建后的图像提供运动修正的过程的流程图。

图6用图表示出了相对于检测器面板的水平平面的患者的三类旋转运动、以及患者的两类平移运动。

图7用图表示出了相对于检测器面板的垂直平面的患者的三类旋转运动、以及患者的两类平移运动。

图8是示出了根据一个实施例的用于确定患者运动矢量的三部分（three-part）过程的流程图。

图9是示出了根据一个实施例的图像修正方法的流程图。

具体实施方式

在详细解释本发明的任何实施例之前，要理解的是，本发明在其应用时不限于在下面的描述中所阐述的或者在附图中示出的构造的细节以及部件的布置。本发明能够具有其他实施例，并且能够以各种方式实践或实施。

还应予以注意的是，可利用多个基于硬件和软件的装置、以及多个不同结构的部件来实现本发明。此外，如随后的篇章中所描述的，图中所示具体配置旨在例示本发明的实施例。替代配置也是可能的。

现代X射线过程使用一个或多个检测器而非膜。检测器以电子方式检测并量化到达检测器的X射线的量（即，X射线的密度）。使用检测器，建立X射线计算机断层术（“CT”）系统，其使得X射线源绕着患者旋转，并且利用在患者的与X射线源相对的一侧的单个宽条状（single-wide strip）检测器元件，以电子方式检测所得X射线。X射线源、检测器、以及允许X射线源旋转的机械结构的组合已知为“CT台架”。针对所有不同的X射线位置收集来自检测器的数据，然后在已知为重建的过程中对其进行组合。组合后的图像表示穿过患者的单个有限宽度“切片”，其中每个点处的图像密度表示特定物理定位处的组织的X射线密度。重建过程使用CT台架的固定几何形状和相对于患者的X射线源-检测器面板组合的可变化角度（即，θ）来处理所收集的投影图像。

通过在图像采集或扫描之间移动患者、X射线源、或者两者的同时重复上述过程，可以获取多个切片。例如，移动支承患者的平台且同时也旋转X射线源，则代替数据切片而产生“螺线（helix）”。另外，将检测器条或环的大小或宽度从单行检测器增加到多行（例如，达到256行），允许对更多数据更加快速的采集。此外，用较大的二维检测器替换检测器条，则获取在每个X射线源位置处的整个检测器面板图像，而非仅单条数据。对在数量上达到600或者更多的这些的图像的收集已知为投影图像。每个投影图像表示随着X射线源和检测器绕着患者同步地旋转而从不同的透视或角度对患者的X射线快照。因为锥形X射线束需要覆盖二维检测器，所以这种CT成像已知为锥形束（“CB”）CT成像。图1示意性示出了根据本发明的一个实施例的CB CT成像系统10，且图2示出了CB CT成像系统10的部件及其几何关系和参数。

CB CT成像系统10包括扫描器12和计算机14。扫描器12包括CT台架13，其包括X射线源16和检测器18。X射线源16和检测器18在转动载台20上在彼此对面对准，该转动载台20使X射线源16和检测器18绕着患者22移动。患者22由座位24所支承。图1中所示的成像系统10为牙科成像系统。因此，患者22坐在座位24中并将他或她的下颚放置到托架（rest）26中。当台架13转动以完成对患者的头部的扫描时，托架26保持患者的头部相对静止。

扫描器12将图像数据输出至计算机14。图像数据表示在扫描期间由检测器18所检测到的X射线的强度水平。计算机14连接至控制台30，其包括显示器32和一个或多个输入和/或输出装置，例如键盘34。使用者使用控制台30以与计算机14交互。例如，用户可以使用控制台30以从计算机14请求图像或其他数据。计算机14将所请求的信息提供至控制台30，并且控制台30将该信息显示在显示器32上、将该信息打印至打印机（未示出），和/或将该信息保存至计算机可读存储器模块（未示出）。

图3示意性示出了根据本发明的一个实施例的图1的计算机14。计算机14包括输入/输出接口40、电子处理单元（“EPU”）42、以及一个或多个存储器模块，诸如可通过磁盘驱动器（未示出）访问的计算机可读磁盘、随机存取存储器（“RAM”）模块44、只读存储器（“ROM”）模块46、或其组合。输入/输出接口40从扫描器12接收图像数据，并将该图像数据提供给EPU 42。

在一些实施例中，完整的台架旋转耗时大约8~40秒。在此期间，患者可能会移动或CT台架可能会意外地移动，从而导致所得图像的模糊。典型的图像分辨率在0.25毫米的量级。因此，该相同量级的患者运动经常导致图像模糊，且大范围的患者移动可能使得所得图像对于其所期望的门诊的目的而言变得不可接受。另外，即使当模糊并不过度时，模糊也会导致图像质量普遍下降。例如，台架振动会导致模糊以及降低的图像质量。

计算机14通过跟踪图像中刚性或半刚性物体的移动来针对患者移动对图像进行修正。例如，在其中没有患者移动的理想条件下，被成像的刚性物体随着台架绕着患者旋转而以明确定义的方式改变投影图像的二维上的定位。患者移动（或预期之外的台架移动）导致图像中物体的预期定位之间的偏离。因此，通过测量明确定义的物体与其期望位置的偏离，可以测量并修正患者（以及预期之外的台架）移动量。具体而言，如下所述，如果图像序列中存在至少三个物体，则可以组合所测量的这些物体与其预期定位的偏离来确定患者运动矢量，其可以应用于图像以针对患者移动进行修正。

为了确保在图像中有所期望的数量的明确定义的刚性物体，可以在扫描之前将基准标志（例如，三个或多个）置于患者上。这些标志一般包括铅或钢BB，其较为致密并可防止或限制X射线透过。但是，这些标志也可以由其他材料制成，并且构造成在扫描期间产生的投影图像中以相对高比例可见的其他的形式或形状。每个标志或多个标志可位于粘合剂的层之间，并且该粘合剂可以涂覆至患者，以确保标志在程序期间不移动。另外，尽管更为复杂，但也可以使用内部的解剖标志物（landmark）作为标志物，而非外部应用的标志。

将标志置于患者上，使得由CB CT成像系统所生成的每个视野或图像均包括至少三个标志。例如，在患者上可放置七到九个标志，以确保至少三个标志在每个图像面内，以减少位置测量噪音，并允许结果的统计组合。在一些实施例中，在患者上标志均匀地分隔开，以关于患者具有最大的空间分布，并避免干扰图像判读（interpretation）。

在将标志置于患者上之后，对患者进行扫描，诸如患者的头部，并将表示投影图像序列的所得图像数据传输至计算机14。图4示出了包括八个标志102的示例性投影图像100。应予以注意的是，“标志102”是标志的表示而非实际的标志（即，标志102不是实际的bb或类似的装置）。使用术语标志来指代实际标志和标志的表示或图像与本领域技术人员所使用的普通语法学句法相一致。因此，在下面的书面描述中，并非总是对实际的标志和标志的表示之间做出明确区分。

EPU 42接收图像数据并通过执行一个或多个应用软件或模块来处理信息。应用软件或模块储存在存储器模块中，诸如ROM模块46或计算机可读磁盘、或者可由EPU 42通过驱动器或外部接口而访问的介质。如图3所示，存储器模块储存标志处理模块50、运动矢量提取模块52、投影矩阵计算模块54和重建模块56。图5示出了方法60，该方法60示出了经过这四个模块的数据流。如图5所示，EPU 42取得来自输入/输出接口40的投影图像，并执行标志处理模块50，以识别每个投影图像上的标志点、识别每个所识别的标志的垂直和水平位置、将每个位置及其维度分配给适当的物理标志、并使用图像位置数据来估计每个标志的基准三维（“3D”）物理定位（步骤62和64）。在所分配的代理案卷号为No.026212-9015-00的、标题为“X射线图像中的标志识别和处理（MARKER IDENTIFICATION AND PROCESSINGIN X-RAY IMAGES）”的共同未决申请中描述了此标志识别和处理过程的各种实施例。

如图5所示，在EPU 42获得用于每个投影图像的标志点列表之后，EPU 42执行患者运动矢量得出模块52，以得出用于各个投影图像的六分量患者运动矢量（步骤66和67）（参见图6-8）。EPU 42接着执行投影矩阵计算模块54，以基于六分量患者运动矢量以及台架几何参数而确定用于各个投影图像的运动修正后的投影矩阵（步骤68）（参加图9）。然后EPU 42执行重建模块54，以基于运动修正后的投影矩阵们而生成运动修正体积图像（步骤69）（参见图9）。下面详细地描述方法60的各部分。

患者运动矢量得出

如图5所示，方法60的第一步骤包括获得投影图像中的每个标志的标志位置信息（步骤62）。在所分配的代理案卷号为No.026212-9015-00的、标题为“X射线图像中的标志识别和处理（MARKER IDENTIFICATION AND PROCESSING IN X-RAYIMAGES）”共同未决申请中公开了用于确定投影图像中的每个标志的标志位置信息的各种实施例。所述标志位置信息包括每个投影图像上的每个标志的U（行）和V（列）定位。每个标志的U和V定位提供与扫描器的检测器面板有关的标志的2D位置（即，U_M1、V_M1、U_M2、V_M2、U_M3和V_M3）。这些测量值是在存在患者运动时所检测的标志的实际测量的位置。如图6和7所示，第一维度U与CT台架的旋转圆形相切。第二维度V是与CT台架的旋转轴线相平行的平面。

在步骤64，方法60基于标志的U和V位置信息而确定用于各个标志的基准3D物理定位（即，X、Y和Z）。由标志的基准位置判断各个标志的运动。基准位置是任意的，并且可以各种方式加以定义。在一个实施例中，其定义为扫描开始时的标志的位置。该位置可以各种方式得出，诸如发现用于生成最佳拟合所观察的标志的轨迹的理论轨迹的X、Y和Z坐标。当使用三个标志时，用于这些标志的基准定位由标注为X₁、Y₁、Z₁、X₂、Y₂、Z₂、X₃、Y₃和Z₃的九个坐标值（即，三个标志乘以各个的三个坐标）加以定义。

下列方程式用数学方式描述了定位X_i、Y_i和Z_i处的基准物理3D点是如何投影到特定台架旋转值θ_j处的2D检测器平面上的：

$U_{\mathrm{ij}}=\frac{\mathrm{DSD}*(Y_i\cos {\mathrm{\θ}}_j-X_i\sin {\mathrm{\θ}}_j)}{(\mathrm{DSO}+Y_i\sin {\mathrm{\θ}}_j+X_i\cos {\mathrm{\θ}}_j)}$ （方程式1）

$V_{\mathrm{ij}}=\frac{\mathrm{DSD}*Z_i}{(\mathrm{DSO}+Y_i\sin {\mathrm{\θ}}_j+X_i\cos {\mathrm{\θ}}_j)}$ （方程式2）

下标i表示标志编号，下标j表示投影编号。DSO是从X射线源到台架旋转中心的距离，且DSD是从X射线源到检测器面板的距离。变量θ_j是当前台架旋转角度。例如，当台架源位于对象的头部正后方时，θ_j为零。图2中大体上示出了参数DSO、DSD和θ。

接着在步骤66，方法60确定含有用于特定投影图像的六个运动变量（即，ΔX、ΔY、ΔZ、α、β和γ）的患者运动矢量。通常，当所扫描的对象被视为刚体时，该步骤确定描述该扫描对象（例如，患者的头部）的总体运动的用于投影图像的六分量运动矢量。每个运动矢量均考虑了投影图像中各个标志的所观察的位置（U和V）和相应各个标志的预测定位之差，所述预测定位假定理想的台架运动以及各个标志的已知的、固定（例如，非移动）的物理X、Y和Z基准定位。

例如，在该步骤，计算机接收用于每个标志的对象列表作为输入，所述对象列表列出了用于在每个投影图像中的标志的位置信息。如果每个投影图像中存在至少三个标志，则计算机基于用于该图像中所包含标志的标志对象列表，而输出用于各个投影图像的表示用于该图像的所检测到的运动的六分量运动矢量。六个分量中的三个表示患者在X、Y和Z维度上的平移移位（shift）或运动（即，变量ΔX、ΔY和ΔZ），另外三个分量表示患者关于X、Y和Z轴的旋转移位或运动（即，变量α、β和γ）。

如图5所示，不依赖于其他投影图像，为各个投影图像确定患者运动矢量（步骤67）。因此，步骤66和67的最后结果是包含等于投影图像的数量的大量向量的一组运动矢量。例如，对600个投影图像的研究将产生具有600个独立向量的运动矢量组。然而，如果一个或多个投影图像未包含三个或更多个标志，则所产生的向量的实际数量会较少，这会妨碍方法60得出运动信息。

应予以注意的是，在一些实施例中，运动矢量是投影图像的准确运动矢量的估计值。例如，在扫描期间，每个标志在其X、Y和Z位置处经历不同变化。方法60使用每个投影图像中每个标志的不同运动以生成六分量运动矢量。如果存在三个标志，则可以使用六个方程式（即，用于三个标志中的各个的U和V方程式）和六个未知量（即，运动矢量的分量）来确定准确的答案或运动矢量。然而，如果图像中存在多于三个的标志，则与未知量相比，存在更多方程式，且该问题是超定的（overspecified）。另外，标志测量由于噪音而不能总是准确测量。因此，由于问题是超定的且存在测量噪音，所以通过方法60所确定的运动矢量是用于特定投影图像的准确运动矢量的估计值。

回到步骤66，可以对各种方程式和公式加以定义，以确定患者运动矢量。例如，如果不存在患者运动，则用于特定标志的X_i、Y_i和Z_i对于所有投影图像都为常量，这暗示U和V仅由于变化的θ而改变。然而，当存在患者运动时，X_i、Y_i和Z_i由于运动的定位影响而改变。因此，在存在患者运动时，各个坐标均由下列方程式更加适当地表示：

X_i=X_iref+x_ij    Y_i=Y_iref+y_ij  Z_i=Z_iref+z_ij

其中，x_ij、y_ij和z_ij表示特定标志的i位置与其在与患者运动相关联的特定投影j处的基准定位之差。将这些值代入上面所示的方程式1和2中，得到下列修正后的方程式：

$U_{\mathrm{MEASij}}=\frac{\mathrm{DSD}*\left(\right(Y_{\mathrm{refi}}+y_{\mathrm{ij}})*\cos {\mathrm{\θ}}_j-(X_{\mathrm{refi}}+x_{\mathrm{ij}})*\sin {\mathrm{\θ}}_j)}{(\mathrm{DSO}+(X_{\mathrm{refi}}+x_{\mathrm{ij}})*\cos {\mathrm{\θ}}_j+(Y_{\mathrm{refi}}+y_{\mathrm{ij}})*\sin {\mathrm{\θ}}_{\text{j}})}$ （方程式3）

$V_{\mathrm{MEASij}}=\frac{\mathrm{DSD}*(Z_{\mathrm{refi}}+z_{\mathrm{ij}})}{(\mathrm{DSO}+(X_{\mathrm{refi}}+x_{\mathrm{ij}})*\cos {\mathrm{\θ}}_j+(Y_{\mathrm{refi}}+y_{\mathrm{ij}})*\sin {\mathrm{\θ}}_j)}$ （方程式4）

对于使用三个标志时的情况，在特定θ_j处，存在六个方程式（用于三个标志中的各个的U和V）和九个未知量（用于三个标志中的各个的x_ij、y_ij和z_ij），其在没有附加信息的情况下，是不可能求解的。

然而，如果增加假定所有标志都固定在刚体内或刚体上这一限制，则任何特定标志的运动都可以由其基准定位以及由六个变换参数（即，运动矢量分量）ΔX、ΔY、ΔZ、α、β和γ所提供的整个物体的运动来完整描述。变量ΔX、ΔY和ΔZ定义投影图像中患者的平移运动，且变量α、β和γ表示投影图像中患者的旋转运动，如图6和7所示。因此，通过增加刚体限制，至此仅存在六个未知量和六个方程式，这使问题可求解。

为了进一步简化方程式和公式，用于物体的基准框架可以从基准的固定（即，相对于地面）框架改变为基准的旋转的台架框架。可以将固定到旋转框架变换与患者运动变换组合，以提供总体运动变换。下列矩阵提供该组合：

|Γ_transformed|=|A_overall||Γ_ref|=|θ_rotation||A_{patient_motion}||Γ_ref|（方程式5）

其中，|Γ|是包含分量X、Y和Z的物体的位置向量，|А_{patient_motion}|是包含与ΔX、ΔY、ΔZ、α、β和γ平移和旋转相关联的分量的运动变换矩阵，|θ_rotation|是用于特定的投影图像的与台架旋转相关联的旋转矩阵。这暗示以上的方程式3和4的替代公式将是：

$U_{\mathrm{MEASij}}=\frac{\mathrm{DSD}*Y_{\mathrm{transformes\; j}}}{\mathrm{DSO}+X_{\mathrm{transformed}}}$ （方程式6）

$V_{\mathrm{MEASij}}=\frac{\mathrm{DSD}*Z_{\mathrm{transformed\; j}}}{\mathrm{DSO}+X_{\mathrm{transformed}}}$ （方程式7）

其中：

X_transformedj=f₁(X_refi,Y_refi,Z_refi,ΔX,ΔY,ΔZ,α,β，γ,θ)（方程式8）

Y_transformedj=f₂(X_refi,Y_refi,Z_refi,ΔX,ΔY,ΔZ,α,β，γ,θ)（方程式9）

Z_transformedj=f₃(X_refi,Y_refi,Z_refi,ΔX,ΔY,ΔZ,α,β，γ,θ)（方程式10）

并且f₁、f₂和f₃由上面的变换方程式（即，方程式5）定义。

因此，给出“固定”参数DSD、DSO、X_ref1、Y_ref1、Z_ref1、X_ref2、Y_ref2、Z_ref2、X_ref3、Y_ref3以及Z_ref3以及特定θ处的U_MEAS1、V_MEAS1、U_MEAS2、V_MEAS2、U_MEAS3、以及V_MEAS3的测量值，方法60使用方程式5、6和7得出在该θ处，以及用于投影图像序列中各个其他θ处的运动参数ΔX、ΔY、ΔZ、α、β和γ。

应予以理解的是，所示出的方程式6和7本身是非线性的。由方程式5、8、9和10所表示的运动变换具有与变换的旋转部分相关联的正弦（sine）和余弦（cosine）分量。而且，由于每个中所固有的除法，因而方程式6和7的同时求解表示非线性操作。所有这些因素都暗示无法得出准确解，并且迭代求解将是过分耗时的。然而，通过采用一些简化，可以确定可行的解。图8是示出了用于确定根据本发明一个实施例的患者运动的三部分方法70的流程图。当EPU 42执行患者运动矢量得出模块时，由计算机14的EPU 42执行方法70。方法70采用两个简化，都是基于实际的患者运动与台架旋转相比相对较小的前提而进行的。

第一个简化通过对旋转变换时的正弦函数（sine）和余弦函数（cosine）应用泰勒级数展开，并且只考虑每个中的第一项来线性化方程式。有效地，这意味着cos(x)用1代替且sin(x)用x代替。该假定对于小值的旋转分量而言是有效的。然而，如下所述，方法70可以采用允许此限制放宽的迭代方法。

第二个简化有利地利用V方程式对α、β和ΔZ的变化更为敏感而对ΔX、ΔY和γ的变化不太敏感这一点。相反地，U方程式对ΔX、ΔY和γ的变化更为敏感而对α、β和ΔZ的变化不太敏感。这意味着，对于小运动量，总体问题可以通过将问题分解成两部分来加以简化。具体而言，首先，可以对用于α、β和ΔZ的同时线性化后的V方程式进行求解，然后可以将先前使用V方程式得出的用于α、β和ΔZ的值作为常量，来对用于ΔX、ΔY和γ的同时线性化后的U方程式进行求解。

如上所述，使用迭代可以提高求解的精确度。迭代（步骤76）包含重叠地对U方程式（步骤72）和V方程式（步骤74）求解。在每个方程式的解内，存在关于迭代是如何增加精确度的两个分量。首先，V的每个求解均使用由U先前的求解而得出的用于ΔX、ΔY和γ的最近的、精确度逐渐提高的值。同样，U的每个求解均使用由V先前的求解而得出的关于α、β和ΔZ的最近的值（步骤74）。第二分量是残差（residual）的使用。为了做到这一点，将U和V方程式的较高阶（非线性分量）视为伪常量，以允许线性公式。由于相对于正常的台架运动而言患者运动可以认为是小的，所以可以做出该假定。在U或V方程式中任一个的各个求解（步骤71、72或74）之后，将伪常量更新至其精确度逐渐提高的值。对于V方程的第一个求解（步骤71），在执行求解之前没有信息存在以设定用于伪常量或ΔX、ΔY和γ的值。在执行此步骤之前，这些参数被设定为“零”。

下面提供求解V方程式和U方程式的进一步细节。

求解V方程式

为了求解V方程式，注意，分量α定义为绕X轴的旋转，β定义为绕Y轴的旋转，γ定义为绕Z轴的旋转。据此，用α表示的旋转的矩阵形式如下：

$\mathrm{rot\α}=\left[\begin{array}{ccc}1& 0& 0\\ 0& \cos \left(\mathrm{\α}\right)& -\sin \left(\mathrm{\α}\right)\\ 0& \sin \left(\mathrm{\α}\right)& \cos \left(\mathrm{\α}\right)\end{array}\right]$

对此矩阵应用泰勒展开，使该矩阵中的正弦函数和余弦函数线性化。

$\cos \left(\mathrm{\α}\right)=1-\frac{{\mathrm{\α}}^2}{2!}+\frac{{\mathrm{\α}}^4}{4!}+...$

$\sin \left(\mathrm{\α}\right)=\mathrm{\α}-\frac{{\mathrm{\α}}^3}{3!}+\frac{{\mathrm{\α}}^5}{5!}+...$

对于小角度，第一项1和α分别表示cos(α)和sin(α)的合理近似。替代地，可以通过用较高项的估计值修正第一项来提高精确度。然而，在没有现有了解的情况下，这些较高阶项的最佳估计值为零。但是，一旦将用于sin(α)和cos(α)的线性近似用于得出α的近似值，则该α值可以用于改进用于高阶项的近似。换言之，可以通过含有以下两项的表达式来给出cos(α)和sin(α)的各个的表示：线性（α中的）项和非线性残差，该非线性残差在任何迭代期间都可以视为是常量。

例如，下面提供关于cos(α)和sin(α)的数学表示，其中cos(α)已用c1代替，且sin(α)已用s1代替以更加简明。

cos(α)=c1=1-h1

h1=1-c1

sin(α)=s1=a1-g1

g1=a1-s1

其中，h1表示cos(α)残差且g1表示用于sin(α)的残差。后缀数字1用于指明旋转α。

使用后缀数字2表示旋转β，后缀数字3表示旋转γ，且后缀数字4表示台架旋转θ，下列矩阵定义所关注的每个旋转。

$\mathrm{rot\α}=\left[\begin{array}{ccc}1& 0& 0\\ 0& 1-h1& -a1+g1\\ 0& a1-g1& 1-h1\end{array}\right]$ （方程式11）

$\mathrm{rot\β}=\left[\begin{array}{ccc}1-h2& 0& -a2+g2\\ 0& 1& 0\\ a2-g2& 0& 1-h2\end{array}\right]$ （方程式12）

$\mathrm{rot\γ}=\left[\begin{array}{ccc}c3& -s3& 0\\ s3& c3& 0\\ 0& 0& 1\end{array}\right]$ （方程式13）

$\mathrm{rot\θ}=\left[\begin{array}{ccc}c4& -s4& 0\\ s4& c4& 0\\ 0& 0& 1\end{array}\right]$ （方程式14）

注意，a1＝α且a2=β，并且由于其是需要求解的两个变量而已被分解。相反，γ和θ由于这些变量中每个在V方程式求解期间均视为固定参数因而未被分解。

对于运动的平移分量，下列向量表示特定标志的3分量基准位置。

${\mathrm{position}}_{\mathrm{original}}=\left[\begin{array}{c}X0\\ Y0\\ Z0\end{array}\right]$

标志从一个位置到另一个位置的变换可以涉及多达3个位置移位（平移）和3个旋转移位。这些移位和平移的准确级别是任意的，但必须针对具体公式加以定义。例如，运动变换可以定义为首先包含三分量平移移位，接着是α旋转、β旋转和γ旋转，且最后是θ（台架）旋转。此定义用于下列推导中，但其他定义也将等同适合。

可以通过仅增加运动位移（displacement）（XD、YD和ZD）到标志的对应的基准位置，来定义用于特定标志的新的X、Y和Z位置。然而，在此V方程式求解内，X和Y位置假定为常量，并可以通过使用两个符号XS和YS来简化。使用这些符号，下列向量表示标志的平移后的位置（适用于平移移动而非旋转之后的位置）。

${\mathrm{position}}_{\mathrm{tramslated}}=\left[\begin{array}{c}X0+\mathrm{XD}\\ Y0+\mathrm{YD}\\ Z0+\mathrm{ZD}\end{array}\right]=\left[\begin{array}{c}\mathrm{XS}\\ \mathrm{YS}\\ Z0+\mathrm{ZD}\end{array}\right]$

在施加了旋转之后，至此标志处于其完全运动位移后的定位。这一新的位置被称为变换后的位置向量，并在下列方程式中给出。

$\left[\begin{array}{c}{X}_{\mathrm{tramsformed}}\\ Z_{\mathrm{tramsformed}}\\ Z_{\mathrm{tramsformed}}\end{array}\right]=\mathrm{rot\θ}\×\mathrm{rot\γ}\×\mathrm{rot\β}\×\cot \mathrm{\α}\×\left[\begin{array}{c}\mathrm{XS}\\ \mathrm{YS}\\ Z0+\mathrm{ZD}\end{array}\right]$ （方程式15）

此时，V方程式可以求解。例如，上面给出了V方程式作为方程式6，这里重复作为方程式16：

$V=\frac{\mathrm{DSD}*Z_{\mathrm{transformed}}}{\mathrm{DSO}-X_{\mathrm{transformed}}}$ （方程式16）

方程式16表示用于一个特定标志的V和该特定标志的运动变换后的标志定位之间的关系。如果将方程式16与方程式11、12、13、14和15组合，则所得方程式包含之后求解的三个变量：a1、a2和ZD（对应于α、β和ΔZ）。

0＝-DSD*YS*g1-DSD*YS*h2*a1+DSD*h2*h1*Z0+V1*c4*c3*g2*Z0+V1*XS*c4*c3-V1*YS*c4*s3+DSD*h2*h1*ZD-V1*YS*s4*c3+V1*s4*s3*a2*ZD-V1*DSO-DSD*h1*Z0-DSD*h2*ZD+V1*YS*s4*s3*a2*a1-V1*XS*s4*s3-DSD*XS*g2+DSD*YS*a1-DSD*h2*Z0+DSD*XS*a2-DSD*h1*ZD+DSD*YS*h2*g1+V1*c4*c3*g2*ZD+V1*YS*a2*c4*c3*g1-V1*s4*s3*g2*ZD-V1*s4*c3*g1*ZD+V1*YS*c4*c3*g2*a1-V1*YS*a2*s4*s3*g1-V1*YS*s4*s3*g2*a1-V1*YS*c4*c3*g2*g1+DSD*Z0+DSD*ZD-V1*c4*s3*g1*ZD-V1*s4*c3*g1*Z0+V1*YS*s4*s3*g2*g1-V1*s4*s3*g2*Z0-V1*c4*s3*g1*Z0+V1*c4*s3*a1*Z0+V1*c4*s3*a1*ZD+V1*s4*c3*a1*Z0+V1*s4*c3*a1*ZD-V1*c4*c3*g2*h1*Z0-V1*s4*s3*a2*h1*Z0-V1*s4*s3*a2*h1*ZD+V1*XS*s4*s3*h2-V1*XS*c4*c3*h2+V1*YS*c4*s3*h1-V1*c4*c3*a2*Z0-V1*c4*c3*a2*ZD+V1*s4*s3*a2*Z0-V1*YS*c4*c3*a2*a1-V1*c4*c3*g2*h1*ZD+V1*s4*s3*g2*h1*Z0+V1*s4*s3*g2*h1*ZD+V1*YS*s4*c3*h1+V1*c4*c3*a2*h1*Z0+V1*c4*c3*a2*h1*ZD

对于每个标志均存在一个这样的方程式。

为了求解上述方程式，其必须设置为下列形式：

a1*A+a2*B+ZD*C+D=0    （方程式17）

这通过将在不高于一次幂（power）状态下只含有三个变量之一a1、a2或ZD以及不含任何正弦函数或余弦函数残差g1、g2、h1或h2的这些子表达式分为一组来完成。这三组子表达式变为系数A、B和C。所有其他项都一起集合到“常数”项D中。

一旦进行此分类，就得到方程式17所示出的形式的m线性方程组（每个标志具有一个方程式）。如果存在3个标志，则通过简单的矩阵求逆而得出解。如果存在多于三个的标志，则可以使用奇异值分解（“SVD”）而得出解。SVD有效地提供令与每个方程式相关联的方差最小化的、用于a1、a2和ZD的值。

对于第一次迭代，构成包含a1、a2和/或ZD的较高阶组合或者含有正弦函数或余弦函数残差g1、g2或h1或h2的D的子表达式假定为零。然而，在随后的迭代中，可以提供用于这些较高阶项的改进的估计值。

求解U方程式

对于分别表示ΔX、ΔY和γ的变量XD、YD和a3的求解，可以使用提供给出特定的变换后X和Y位置的U位置的方程式，而以与求解V方程式完全类似的方式来实现：

$U=\frac{\mathrm{DSD}*Y_{\mathrm{transformed}}}{\mathrm{DSO}-X_{\mathrm{transformed}}}$

在此情况下，方程式4仍然相同，而方程式11、12、13和15变成：

$\mathrm{rot\α}=\left[\begin{array}{ccc}1& 0& 0\\ 0& c1& -s1\\ 0& s1& c1\end{array}\right]$ （方程式18）

$\mathrm{rot\β}=\left[\begin{array}{ccc}c2& 0& -s2\\ 0& 1& 0\\ s2& 0& c2\end{array}\right]$ （方程式19）

$\mathrm{rot\γ}=\left[\begin{array}{ccc}1-h3& -s3+g3& 0\\ s3-g3& 1-h3& 0\\ 0& 0& 1\end{array}\right]$ （方程式20）

$\left[\begin{array}{c}{X}_{\mathrm{tramsformed}}\\ Z_{\mathrm{tramsformed}}\\ Z_{\mathrm{tramsformed}}\end{array}\right]=\mathrm{rot\θ}\×\mathrm{rot\γ}\×\mathrm{rot\β}\×\cot \mathrm{\α}\×\left[\begin{array}{c}X0+\mathrm{XD}\\ Y0+\mathrm{YD}\\ \mathrm{ZS}\end{array}\right]$ （方程式21）

在替换之后，必须储存所得方程式，以提供下列形式的方程式：

$\mathrm{XS}*E+\mathrm{YD}*F+a3*G+H=0$ （方程式22）

再者，对于每个标志将存在这些方程式之一。可以利用矩阵求逆或利用奇异值分解来求解这些方程式。

尽管方程式17和22可以按任一次序进行求解，但优选的实施方式为先求解方程式17之后求解方程式22，然后通过这两个方程式迭代两次或更多次。已发现两次迭代即可提供用于运动的合理有用的估计值，但优选具有更多次迭代。利用仅五次迭代就可以获得精确度（比六个小数位数更好）。因为当N典型地小于10时，每次迭代均包含仅反转3×N个矩阵，所以可以快速地执行总体求解。

附录A中提供了根据本发明一个实施例的上述用于求解U和V方程式的方法的示例实施方式，其使用的是Matlab编程语言。在Intel（英特尔）Core 2U7600CPU上，用于执行对于300个投影中的每个处的24个标志和5次迭代的未经编译的Matlab代码的计算时间大约为0.2秒，这表示用于每个投影的近似求解时间为7毫秒。该计算时间中大部分与奇异值分解矩阵求解相关联。

回到图5，在为投影图像序列确定六分量运动矢量之后，EPU 42在步骤68和69处检索并执行投影矩阵计算模块54和重建模块56，以生成修正的和重建后的投影图像。这些步骤在下一节中详细描述（参见图9）。

应注意的是，在上述对方法70的描述中，假定了台架以完美的方式（即，除期望旋转以外，不存在台架旋转）绕着患者旋转。然而，可放宽此条件，且可修改方法70以考虑计划之外的台架运动。也可修改方法70，以适用于更加复杂的台架几何形状和旋转。

此外，在方法70中，标志无需在放置于患者上时彼此之间具有固定关系。标志之间的关系可以由标志轨迹行为的分段式（piecewise）分析而得出，只要患者运动不太严重。然而，在一些实施例中，一旦被放置，标志就需要相互维持固定关系。这意味着方法70更好地适用于刚体，例如人的头部。然而，如果采用更多标志，则方法70也可用于刚性较少的体。例如，只要三个（或更多个）标志对于特定区域具有固定关系，则可以对超过一个的区域进行评价。这一构思对于当颌正在移动时分离地评价头部中下颌和上颌的运动是有用的。另外，如果三个标志固定在刚性框架上，则也可以评价该框体的运动。此方法对于解剖学上的可变形运动，例如呼吸是有用的。

另外，方法70的标志放置就测量误差而言通常很重要，这是因为标志放置影响方法70的结果的精确度。例如，标志之间分得越开，测量误差将越小。在一个实施例中，围绕患者的头部180度放置标志，这提供大约0.2毫米或更好的运动矢量分量检测精确度。

因为方法70为每个投影图像产生独立的运动估计值，所以该方法提供附加的灵活性。首先，如果在特定投影图像中检测到少于三个的标志，则该方法可以由相邻的投影图像对该投影图像中的运动进行内插（interpolate）。其次，该方法可以使用与每个投影图像有关的数据来对相邻的投影结果进行滤波，并且进一步提供该过程的修正作用。

此外，因为方法70使用成像数据来提取运动信息，所以所检测的运动包括患者移动以及未经修正的或者计划之外的台架移动这两者。可以使用关于计划之外的台架移动的信息来动态地修正台架运动。例如，使用方法70，修正可以在患者扫描（例如，如果患者被用仪器装备）期间进行，或者作为使用多个标志人体模型（phantom）的患者扫描之间的校准过程的一部分。

图像修正和重建

回到图5，在投影图像中检测运动之后，计算机14基于所检测的运动重建并修正这些图像（步骤68和69）。图9示出了根据本发明一个实施例的执行这些步骤的方法80。当EPU 42执行投影矩阵计算模块54和重建模块56时，由计算机14的EPU 42执行图像修正方法80。

通常，方法80通过将关于所检测的运动的信息结合到重建过程中来针对运动修正投影图像。方法80处理在图像中所检测到的非可变形运动，其由六分量向量（即，患者运动矢量）针对每个投影图像加以规定。如上所述，六分量向量包括用于X、Y和Z平移运动的三个参数（即，ΔX、ΔY和ΔZ）和用于绕X、Y和Z轴的旋转运动的三个参数（即，α、β和γ）。各个向量表示患者运动（例如，在外部测量的运动的情况下）或患者相对于台架的运动（例如，如果运动信息是从图像推知的）。位置误差向量可表示患者运动相关误差、动态台架位置误差、或者两者。

回到图9，作为方法80的初始步骤，获取CB CT投影图像（步骤82）。接着，为各个投影图像确定位置误差向量（即，患者运动矢量）（步骤84）。可通过包括外部检测器和图像分析器的各类装置，来确定位置误差向量。在一个实施例中，方法80使用上面参照图5-8描述的方法60来完成该步骤。

接着，方法80将位置误差向量与描述扫描器几何形状的静态校准参数或信息（例如，在扫描过程中恒定不变的）参数组合，以得出用于各个投影图像的投影矩阵（步骤86）。用于扫描器几何形状的参数包括从X射线源到旋转中心的距离（DSO）、从X射线源到检测器的距离（DSD）、以及各个投影图像的θ值。投影矩阵描述各个X射线是如何从X射线源行进穿过患者而抵达检测器面板上的特定点（spot）的。投影矩阵等同于用于在3D计算机图形操作中在计算机屏幕上将虚拟3D体积转换成2D图像的矩阵。示例性投影矩阵如下所示：

ξ=A00*X+A01*Y+A02*Z+A03*1.0

ψ=A10*X+A11*Y+A12*Z+A13*1.0

ζ=A20*X+A21*Y+A22*Z+A23*1.0

U=ξ/ζ

V=ψ/ζ

坐标X、Y和Z表示患者上的点在3D空间中的位置，U和V表示相同点在2D检测器面板上的位置，且A00到A23表示用于给定投影图像的投影矩阵的分量。

在生成投影矩阵之后，方法80将投影图像以及对应的投影矩阵提供给重建过程（88）。所述重建过程使用这些投影图像和投影矩阵以确定在称为背投（back-projection）的期间如何将每个投影图像像素映射到投影体积上。替代地，重建过程可以使用迭代的重建来修正并重建图像，其使用投影矩阵以将像素数据映射到体积数据。重建过程输出重建后的图像，这些图像是基于所检测的运动（数据90）而修正后的。

即使在没有测量患者运动时，方法80也提供针对动态台架误差进行修正的机制，只要台架误差在扫描之间是可重复的即可。方法80也可以通过使用图形处理单元（“GPU”）或其他专门的并行处理硬件而进一步加速。另外，图像修正方法80可适于具有不完整的运动数据的情形。例如，如果标志的位置仅在扫描开始和结束处是已知的，则中间投影图像可通过利用其他投影图像的已知的位置误差向量而对与该投影图像相关联的位置误差向量进行线性内插来加以修正。

虽然已经结合CB CT成像对患者运动矢量得出方法60和图像修正方法80进行了描述，但这些方法也可以用于CT、MRI、超声、其他形式的医学成像以及非医学成像，例如摄影。例如，图像修正方法80可以用于其中将3D体积作为2D图像序列加以评价的成像技术，包括3D超声和3D光学一致性X线断层摄影术（coherence tomography）。

此外，虽然图3将ROM模块46示出为存储有单独的模块（例如，50、52、54、56等），但也可以组合或分离储存在ROM模块46或其他存储器中的模块。例如，在一些实施例中，相同的模块执行患者运动矢量得出方法60和图像修正方法80。另外，ROM模块46可以包含用于执行除上述的那些功能以外的功能的其他模块。此外，在一些实施例中，模块可以存储在可读磁盘上，并传送至RAM模块44和/或ROM模块46以进行执行。

因此，除了其他因素之外，本发明提供用于确定图像中的患者运动以及基于所检测的患者运动修正并重建图像的方法和系统。本发明的各种特征和优点将在随附权利要求书中阐明。

附录A

1％Load U and V data(n markers x m projections)into

2load′data\test.mat′；

3

4％Load Reference 3D X，Y，and Z positions of each marker(3xn)

5load data\RefLoc.mat；

6

7％Scanner Geometry Parameters

8DSD＝716.59；

9DSO=488.111；

10SCANANGLE=360.03；

11

12％Global Parameters

13numMarkers＝size(u_mm,2)；

14numProj=size(u_mm,1)；

15theta＝

(-30：-(SCANANGLE/numProj)：-30+(SCANANGLE/numProj)*numProj+1))′/180.0＊pi；

16

17％Initialize motion parameters to zero(first guess)

18CU＝zeros(size(u_mm,1)，3)；

19CV＝zeros(size(u_mm,1)，3)；

20

21％The following starts the main iterative loop

22for iter＝1：5

23％＊＊＊＊＊Solve for V associated parameters(Z，Alpha，and Beta)

24for j=1：numProj

25for i＝1：numMarkers

26％The following are″pseudo-constants″in the lineari zed equations

27V1＝v_mm(j，i)；

28XS=RefLoc(1，i)+CU(j，1)；

29YS＝RefLoc(2，i)+CU(j，2)；

30Z0＝RefLoc(3，i)；

31ZD＝CV(j，1)；

32a1＝CV(j，2)；

33a2＝CV(j，3)；

34h1＝1-cos(CV(j，2))；

35h2＝1-cos(CV(j，3))；

36g1＝CV(j，2)-sin(CV(j，2))；

37g2＝CV(j，3)-sin(CV(j，3))；

38s3=sin(CU(j，3))；

39c3＝cos(CU(j，3))；

40s4＝sin(theta(j))；

41c4＝cos(theta(j))；

42

43Alpha_vect(j，i)＝DSD＊YS+V1＊(-s4＊c3＊Z0+c4＊s3＊Z0)；

44Beta_vect(j，i)=DSD＊XS+V1＊(-s4＊s3＊Z0-c4＊c3＊Z0)；

45delZ_vect(j，i)=DSD；

46errorV_A(j，i)＝(V1＊DSO-DSD＊Z0+V1＊(XS＊(-s4＊s3-c4＊c3)+YS＊(c4＊s3-s4＊c3)))...

47...

48+(V1＊(XS＊(c4＊c3＊h2+s4＊s3＊h2)...

49+YS＊(c4＊c3＊g2＊g1-a2＊c4＊c3＊g1-s4＊s3＊g2＊a1-c4＊c3＊g2＊a1+s4＊c3＊h1...

50-a2＊s4＊s3＊g1+s4＊s3＊a2＊a1+c4＊c3＊a2＊a1+s4＊s3＊g2＊g1-c4＊s3＊h1)...

51+ZD＊(c4＊c3＊a2+s4＊c3＊a1+s4＊s3＊g2＊h1-s4＊c3＊g1+c4＊c3＊g2＊h1-c4＊c3＊g2...

52-s4＊s3＊g2+c4＊s3＊g1-s4＊s3＊a2＊h1-c4＊s3＊a1+s4＊s3＊a2-c4＊c3＊a2＊h1)...

53+Z0＊(-s4＊s3＊a2＊h1-c4＊c3＊a2＊h1-s4＊c3＊g1-s4＊s3＊g2...

54+c4＊c3＊g2＊h1-c4＊c3＊g2+c4＊s3＊g1+s4＊s3＊g2＊h1))...

55+DSD＊(XS＊g2+YS＊(h2＊a1-h2＊g1+g1)...

56+ZD＊(h1-h2＊h1+h2)+Z 0＊(h1-h2＊h1+h2)))；

57end

58

59BV=errorV_A(j，：)′；

60AV＝[delZ_vect(j，：)；Alpha_vect(j，：)；Beta_vect(j，：)]′；

61CV(j，：)＝(AV\BV)′；

62end

63

64％＊＊＊＊＊Solve for U associated parameters  (X，Y，and gamma)

65forj＝1：numProj

66for i＝1：numMarkers

67％The following are″pseudo-constants″in the linearized equations

68U1=u_mm(j，i)；

69X0＝RefLoc(1，i)；

70Y0＝RefLoc(2，i)；

71XD＝CU(j，2)；

72YD＝CU(j，2)；

73ZS＝RefLoc(3，i)+CV(j，1)；

74c1＝cos(CV(j，2))；

75s1＝sin(CV(j，2))；

76c2=cos(CV(j，3))；

77s2＝sin(CV(j，3))；

78a3＝CU(j，3)；

79g3＝a3-sin(a3)；

80h3＝1-cos(a3)；

81s4＝sin(theta(j))；

82c4＝cos(theta(j))；

83

84delX_vect(j，i)＝(U1＊c2＊c4-DSD＊c2＊s4)；

85delY_vect(j，i)=(U1＊c1＊s4-U1＊s1＊s2＊c4+DSD＊s1＊s2＊s4+DSD＊c1＊c4)；

86gamma_vect(j，i)=(U1＊(X0＊c2＊s4-Y0＊s1＊s2＊s4-Y0＊c1＊c4+ZS＊s1＊c4-ZS＊c1＊s2＊s4)...

87+DSD＊(X0＊c2＊c4+Y0＊c1＊s4-Y0＊s1＊s2＊c4-ZS＊s1＊s4-ZS＊c1＊s2＊c4))；

88errorU_vect(j，i)=U1＊DSO-U1＊(X0＊c2＊c4-Y0＊s1＊s2＊c4+Y0＊c1＊s4-ZS＊s1＊s4...

89-ZS＊c1＊s2＊c4)...

90-DSD＊(-X0＊c2＊s4+Y0＊c1＊c4+Y0＊s1＊s2＊s4+ZS＊c1＊s2＊s4-ZS＊s1＊c4)...

91...

92+U1＊(XD＊(-c2＊s4＊a3+c2＊s4＊g3+c2＊c4＊h3)...

93+X0＊(c2＊c4＊h3+c2＊s4＊g3)...

94+YD＊(c1＊c4＊a3-s2＊s1＊c4＊h3+s2＊s1＊s4＊a3+s2＊s1＊s4＊g3-

c1＊c4＊g3...

95+c1＊s 4＊h3)...

96+Y0＊(-s2＊s1＊s4＊g3-s2＊s1＊c4＊h3-c1＊c4＊g3+c1＊s 4＊h3)...

97+ZS＊(-s1＊s4＊h3-s2＊c1＊c4＊h3+s1＊c4＊g3-s2＊c1＊s4＊g3))...

98...

99+DSD＊(XD＊(c2＊c4＊g3-c2＊c4＊a3-c2＊s4＊h3)+X0＊(c2＊c4＊g3-c2＊s4＊h3)...

100+YD＊(-c1＊s4＊a3+c1＊c4＊h3+c1＊s4＊g3+s2＊s1＊c4＊a3+s2＊s1＊s4＊h3...

101-s2＊s1＊c4＊g3)...

102+Y0＊(c1＊c4＊h3+s2＊s1＊s4＊h3-s2＊s1＊c4＊g3+c1＊s4＊g3)...

103+ZS＊(-s1＊s4＊g3-s2＊c1＊c4＊g3+s2＊c1＊s4＊h3-s1＊c4＊h3))...

104；

105end

106

107AU＝[delX_vect(j，：)；delY_vect(j，：)；gamma_vect(j，：)]′；

108BU=errorU_vect(j，：)′；

109CU(j，：)＝(AU\BU)′；

110end

111end

112

113figure；plot([CV(：，2)，CV(：，3)，CU(：，3)]＊180/pi)；

114title(′Motion Angles：Alpha，Beta，and Gamma′)；

115

116figure；plot([CU(：，1)，CU(：，2)，CV(：，1)])；

117title(′Motion Translations：X，Y，and Z′)；

Claims (13)

1.一种用于确定图像中的患者运动的方法，所述方法由包括带有电子处理单元和存储可由所述电子处理单元执行的运动矢量提取模块的存储器模块的计算机的成像系统来执行，所述方法包括：

在所述计算机处，基于在扫描期间由扫描器产生的图像数据获得图像，所述图像包含被假定为具有刚性或半刚性结构的至少三个标志，所述至少三个标志中的各个具有在第一维度和第二维度上的、在扫描器的检测器面板上的实际测量的位置；

利用所述电子处理单元，确定用于所述至少三个标志中各个的基准三维位置；

利用所述电子处理单元，定义用于描述所述至少三个标志中各个的所述基准三维位置和所述实际测量的位置之间的关系、所述扫描器的几何参数、以及患者运动的方程式；以及

利用所述电子处理单元，对所述方程式进行求解，以为所述图像得出用于描述患者运动的六分量运动矢量，所述六分量运动矢量对所述至少三个标志中各个的所述基准三维位置和所述至少三个标志中相应各个的所述实际测量的位置之间的差加以考虑。

2.根据权利要求1所述的方法，其中，所述至少三个标志中各个具有基准三维位置X_i、Y_i和Z_i和在第一维度和第二维度上在所述扫描器的检测器面板上的测量位置U_m和V_m，其中，所述六分量运动矢量包括ΔX、ΔY、ΔZ、α、β和γ，并且其中，对所述方程式进行数字求解以得出六分量运动矢量包括：

（a）假定ΔX、ΔY和γ的值，

（b）使用ΔX、ΔY和γ的假定值以及定义所述第二维度上的位置的所述方程式其中之一，计算用于ΔZ、α和β的值，以及

（c）使用ΔX、ΔY和γ的计算值以及定义所述第一维度上的位置的所述方程式其中之一，计算用于ΔX、ΔY和γ的值。

3.根据权利要求2所述的方法，其中，确定所述运动参数进一步包括：重复步骤（b）到（c）至少一次。

4.根据权利要求2所述的方法，其中，对所述方程式进行数字求解以得出六分量运动矢量包括：使用由线性第一项和估计的残差组成的泰勒展开近似旋转变换中的正弦函数（sine）和余弦函数（cosine）。

5.根据权利要求4所述的方法，其中，对所述方程式进行数字求解以得出六分量运动矢量进一步包括：在每次重复步骤（b）到（c）期间，更新所近似的正弦函数和余弦函数的残差。

6.根据权利要求4所述的方法，其中，使用所述泰勒展开近似所述旋转变换中的正弦函数和余弦函数包括使用作为常量的估计的残差。

7.根据权利要求1所述的方法，其中，在所述计算机处获得的所述图像包括投影图像。

8.一种用于确定图像中的患者运动的系统，所述系统包括：

扫描器；和

计算机，所述计算机包括电子处理单元、存储器模块、和输入/输出接口，其中，所述存储器模块被配置为存储由所述电子处理单元执行的运动矢量提取模块；

其中，所述扫描器被配置为基于在扫描期间由所述扫描器所产生的图像数据获得图像并将所述图像发送到所述计算机，所述图像包含被假定为具有刚性或半刚性结构的至少三个标志，所述至少三个标志中的各个具有在第一维度和第二维度上的、在所述扫描器的检测器面板上的实际测量的位置，

其中，所述电子处理单元被配置为确定用于所述至少三个标志中各个的基准三维位置，

其中，所述电子处理单元被配置为定义用于描述所述至少三个标志中各个的所述基准三维位置和所述实际测量的位置之间的关系、所述扫描器的几何参数、以及患者运动的方程式，

其中，所述电子处理单元被配置为对所述方程式进行求解，以为所述图像得出用于描述患者运动的六分量运动矢量，所述六分量运动矢量对所述至少三个标志中各个的所述基准三维位置和所述至少三个标志中相应各个的所述实际测量的位置之间的差加以考虑。

9.根据权利要求8所述的系统，其中，所述至少三个标志中的各个具有基准三维位置X_i、Y_i和Z_i和在第一维度和第二维度上的在所述扫描器的检测器面板上的测量位置U_m和V_m，其中，所述六分量运动矢量包括ΔX、ΔY、ΔZ、α、β和γ，以及所述电子处理单元被配置为通过如下步骤对所述方程式进行数字求解以得出六分量运动矢量：

（a）假定ΔX、ΔY和γ的值，

（b）使用ΔX、ΔY和γ的假定值以及定义所述第二维度上的位置的所述方程式其中之一，计算用于ΔZ、α和β的值，以及

（c）使用ΔX、ΔY和γ的计算值以及定义所述第一维度上的位置的所述方程式其中之一，计算用于ΔX、ΔY和γ的值。

10.根据权利要求9所述的系统，其中，所述电子处理单元进一步被配置为通过重复步骤（b）到（c）至少一次来确定所述运动参数。

11.根据权利要求9所述的系统，其中，所述电子处理单元被配置为通过使用由线性一次项和估计的残差组成的泰勒展开近似旋转变换中的正弦函数和余弦函数，来对所述方程式进行数字求解以得出六分量运动矢量。

12.根据权利要求11所述的系统，其中，所述电子处理单元进一步配置为通过在每次重复步骤（b）到（c）期间更新所近似的正弦函数和余弦函数的残差，来对所述方程式进行数字求解以得出六分量运动矢量。

13.根据权利要求11所述的系统，其中，所述电子处理单元被配置为通过使用作为常量的估计的残差来使用所述泰勒展开近似所述旋转变换中的正弦函数和余弦函数。

Images