CN101027677A

CN101027677A - 用于校正连续视图之间的对象－探测器相对运动的方法

Info

Publication number: CN101027677A
Application number: CNA2005800211923A
Authority: CN
Inventors: J·R·拉恩; A·C·内尔森
Original assignee: VisionGate Inc
Current assignee: VisionGate Inc
Priority date: 2004-06-24
Filing date: 2005-06-06
Publication date: 2007-08-29
Anticipated expiration: 2025-06-06
Also published as: JP4841549B2; ES2559411T3; CA2570730A1; CN101027677B; US7542597B2; EP1769423A4; HK1104633A1; US20070071357A1; AU2005267569B2; EP1769423B1; AU2005267569A1; US20050010108A1; WO2006011945B1; US7260253B2; EP1769423A2; WO2006011945A2; JP2008504060A; WO2006011945A3; CA2570730C

Abstract

用于三维光学断层摄影成像的运动校正。照射感兴趣对象(1)以产生图像(111)。确定图像的侧向偏移校正值(114)。确定图像的轴向偏移校正值(115)。对图像应用侧向偏移校正值和轴向偏移校正值以产生校正的文件图像(116)。

Description

用于校正连续视图之间的对象-探测器相对运动的方法

相关的申请：

本申请是申请人为内尔森、申请日为2002年4月19日、发明名称为“小对象的可变运动光学断层摄影”、序列号为No.10/126,026的共同未决(co-pending)美国专利申请的部分继续申请，上述申请被并入在此以作参考。

技术领域

本申请涉及成像，更特别地，涉及成像系统中对象-探测器的相对运动的探测和校正，其中，通常来自不同位置的连续视图被采集，每个视图代表三维对象的二维投影或伪投影。

背景技术

光学投影断层摄影显微镜(OPTM，optical projection tomographicmicroscopy)系统适合于微观对象，例如生物细胞及其细胞核的高分辨率成像，所述微观对象被嵌入到流体介质中并且被包含在内外径分别为40微米和150微米的微毛细管中。OPTM利用多个视图，每个视图通过围绕与光轴垂直并且与微毛细管轴平行的轴旋转对象及其包含容器而被采集。具有由M×N像素阵列组成的CCD图像传感器的照相机采集通过对象和成像光学器件之后的光，所述成像光学器件在CCD上产生视场(FOV)的放大图像。由于每个视图来自不同的透视，因此每个视图的内容将与其他视图不同。

由于部件的极小尺寸，可能很难将旋转轴(通常与微毛细管的中心轴重合)定位在探测器的FOV的中心。此外，当旋转微毛细管时使其保持同定也非常困难。另外，在视图之间细胞自身可能沿着管轴移动。结果，已经由于管旋转已经改变的每个视图可能额外受到轴向(平行于微毛细管轴)和侧向(垂直于光轴和管轴)平移(translation)。这些侧向平移是除了那些已经存在的用于不在旋转轴上的对象之外。

为了获得精确的3D重建，不管是通过滤过背投影还是其他手段，必需校正轴向运动和不归因于从一个视图到另一个的透视变化的侧向运动部分。此外必需确定探测器FOV中旋转轴所处的位置。

Nowak的美国专利4,858,128描述了一种方法，其中连续场景首先在第一轴中，然后独立地在另一轴中彼此被相关。两个相关的最大值的位置确定两个轴的需要偏移。所述方法不能提供用于区分由于透视变化导致的“自然”侧向平移和由于微毛细管平移导致的“错误”侧向偏移的方法。Nowak专利指出，“估计信号的这利背景分量和从图像数据减去估计值可能是有用的”。

William H.Press等人的《Numerical Recipes in C：The Art ofScientific Computing》，剑桥大学出版社；第二版(1993年1月1日)，描述了通过计算机程序使用快速傅立叶一变换(FFTs)执行两个数据阵列之间的互相关的技术的方法。简单地说，可以通过对每个阵列应用FFT，使结果阵列中的一个与另一个的复共轭相乘，和对结果应用反FFT而获得两个数据阵列(例如图像数据)的互相关。

为了克服当前本领域状态中的不足，本发明的目的是提供一种方法，用于在多视图成像系统中为每个视图寻找微毛细管的中心轴的位置。本发明的进一步目的是提供一种方法，用于在多视图程序系统中探测连续视图之间的对象一探测器相对运动。本发明的进一步目的是提供一种方法，用于在图像数据采集期间校正图像数据以消除由于对象运动引起的误差。本发明的更进一步目的是提供一种成像系统，其产生代表对象的投影或伪投影视图的多个X-Y数据矩阵，用于对象的轴向切片的后续断层摄影重建。可以通过适当地移位后面的数据以使它与前面的数据对准来消除被探测到的运动，反之亦然。

发明内容

本发明提供了一种用于三维光学断层摄影成像的运动校正的装置和方法。照射感兴趣对象以产生图像。为图像确定侧向偏移校正值。为图像确定轴向偏移校正值。将侧向偏移校正值和轴向偏移校正值应用于图像以产生被校正的文件图像。

附图说明

尽管特别地在后附权利要求中阐述了本发明的新颖特征，但是结合下文所描述的附图的以下具体描述将更好地理解和懂得本发明的组织和内容及其其他目的和特征。

图1是根据本发明的教导构造的用于校正连续视图之间的对象一探测器相对运动的方法的实施例的功能框图。

图2是利用图1中所述的实施例的成像系统的侧向校正部分的功能框图。

图3是利用图1中所述的实施例的成像系统的轴向校正部分的功能框图。

图4A示出了在用在本发明的方法的一个例子中的加阈(thresholding)操作之前的细胞图像。

图4B示出了将用在本发明的方法的一个例子中的加阈操作应用到图4A中所示的图像的结果。

图4C示出了图4A-4B中所示的图像的亮度分布的直方图。

图5示意性地示出了用在本发明的一个实施方式中的光学投影断层摄影显微镜(OPTM)系统。

图6A和图6B示意性地示出了包含安装在压电平移装置上的显微镜物镜的光学断层摄影系统的一个实施方式。

图7示出了对于本发明的实施方式所预期的三维(3D)图像重建中使用的图像采集的方法的流程图示例。

图8示意性地示出了用在本发明的实施方式所预期的三维(3D)图像重建中的运动校正偏移的例子。

具体实施方式

在这里参考与生物细胞有关的特定例子描述本发明，然而，应当理解的是这些例子是为了举例说明本发明的原理，但本发明并不限于此。尽管本发明可以用于其他类型的成像系统，例如X射线计算机断层扫描(CT)成像，为了进行具体描述，以下公开直接介绍光学投影断层摄影显微镜(OPTM)系统环境中的本发明。

在以下论述中，当提供数字例子时使用以下假设：

1.每个图像由宽为640像素高为480像素的阵列组成；

2.每个像素包含单一8位(灰度级0-255)亮度值；

3.参考使用微毛细管的OPTM，管轴平行于较短轴(480像素)；

4.参考使用微毛细管的OPTM，管壁分离是530像素；

5.用于寻找侧向偏移(B1)的bins数目为20；

6.用于寻找轴向偏移(B2)的bins数目为2；

7.阵列被零填补(zero-padded)成1024×1024像素。

应当理解的是这些数值仅仅是为了举例说明的目的；在不背离本发明的实质的情况下可以利用其他数值。

现在参考图1，其显示了根据本发明的教导构造的用于校正连续视图之间的对象一探测器相对运动的方法的实施例的功能框图。在该实施例中，生成每个图像的更改(alterd)拷贝，其中最亮的像素被再次指定为零亮度级，而所有其他像素保持与初始图像中相同的亮度。接着，该加阈(thresholded)图像的二维(2D)FFT逐像素地与参考图像的2D FFT的复共轭相乘。然后，沿着平行于微毛细管轴(此后被称为Y轴)的各线对所得到阵列的亮度进行求和，以在垂直于光轴和微毛细管轴的方向(此后被称为X轴)上建立包含关于亮度图案(pattern)的信息的一维(1D)阵列。然后应用1D FFT，并且确定最大值的位置。所述位置确定沿着X轴要施加的偏移量，以便使微毛细管的中心轴的图像位于图像中其期望的位置。

参考图像利用微毛细管壁之间分离的恒定性，并且由所采集图像中得到的已知量所分离的两条最亮线组成；其余参考图像具有零亮度。当管壁的折射率与管内以及管与载玻片/盖玻片组件间的材料匹配时，管壁在伪投影图像中看上去仅仅是模糊的。直方图操作的作用是增强管壁和图像的其余部分之间的对比度。使用预定的管壁分离，结合沿着图像的X轴的已知数目的像素，由于管的旋转和随之而来的透视变化，使得有可能区分管自身的运动和对象在管内的运动。通过基于不变(constant)特征互相关两个图像，我们的方法使跟踪细胞内变化特征的运动的可能性最小化。

互相关法用于确定从Y轴的轴向偏移量。为此，原始图像的拷贝再次被加阈，但是使用不同的标准以确定哪个像素被重置为零亮度。2D FFT被应用于该图像，并且逐像素地与从前一视图得到的加阈图像的2D FFT的复共轭相乘。对结果应用2D FFT，并且沿着与当前图像的侧向校正和先前图像的侧向校正之差对应的线，X轴偏移被确定为互相关函数中的最大值。这与以前的方法截然不同之处在于由Y轴偏移限制X轴偏移；其依赖于Y轴偏移而得出。

不同于侧向校正，轴向校正是迭代过程，因而会有累积误差。然而，只要连续图像之间的透视变化不太大，轴向互相关就能有效地起作用；这对应于旋转中的小角度增量。通过保持小角度增量，空间内容不会变化太大，从而允许互相关跟踪每个图像中的类似特征。由于角度增量也确定3D断层摄影重建的侧向分辨率，因此保持小角度增量以允许互相关正常工作的要求并不是麻烦的要求。

简单地说，通过适当地移位后续图像以使它们与先前图像对准，或反之亦然，本发明的该实施方式消除了微毛细管的轴向和侧向运动的影响。互相关方法用于寻找侧轴上的偏移，然后寻找管轴上的偏移，限制条件是轴向运动的峰值相关必须在确定侧向运动之后进行。

第一步骤111是生成模板图像。具有例如灰度级65，535的两条白线在它们的理想位置处被创建。每条线具有480像素的长度，平行于短图像维度延伸。两条线的位置由长图像维度(640像素)和管壁分离确定，根据经验被确定为530像素。第一条线位于第0行，而第二条线位于第530行。在该实施方式中，模板图像的尺寸可以从640×480扩展到1024×1024以在两个维度中提供零填补；然而，该行为不是本发明必不可少的。

在模板图像上执行2D FFT，从而使实部和虚部被存储在所得到阵列的交替索引(alternating indices)中。因而，对于零填补阵列而言，阵列大小为2048×1024。模板图像现在处于准备使用的形式。

在步骤114寻找侧向偏移。在步骤114中，图像被加阈以便涂黑(black out)背景像素，然后与两条亮线的二进制图像互相关。感兴趣图像经受侧向偏移确定114。为了帮助轴向校正，为每个图像存储D_LAT。

现在参考图2，其显示了利用图1中所述的实施例的成像系统的侧向校正部分的功能框图。涉及找到侧向偏移114中的步骤包括构造图像的灰度直方图，其中bins数目(B1)可以被设置为从2到255的任何一个整数。对本例而言，假设B1＝20。寻找具有最大数目像素的bin(除了第一bin之外，对应于最黑像素)，并且在原始图像的拷贝中具有所述bin的灰度值或更高灰度值的原始图像中的所有像素被设置为等于零。该步骤21的效果是进一步考虑去除背景像素以产生加阈图像。

作为例子，假设输入图像分别具有最小和最大灰度值31和190，从而使每个bin具有八个灰度级的宽度[(1+190+31)/20＝8]。现在进一步假设直方图中的峰值出现在bin#16(即灰度级从151到158)。那么加阈图像将类似于原始图像，除了初始灰度级大于150的所有像素现在具有零灰度级之外。

图4A-4C示出了将这些步骤21应用到图像的效果。可以参考图4A，然后参考图4B进行原始图像与分割图像的视觉比较，其中图4A示出了在分割和加阈之前的细胞图像的例子，图4B示出了与图4A的原始图像对应的分割和加阈细胞图像的例子。图4C示出了原始图像和应用加阈之后的图像的灰度级的比较的图像例子的直方图。

对加阈图像22应用2D FFT，并且其傅立叶变换与模板图像的傅立叶变换的复共轭相乘23。所得到阵列沿着640行中的每一行求和24以计算新阵列，该新阵列经傅立叶变换(在1D中)25以寻找加阈图像和参考图像的行的互相关。定位26和评估28 1D阵列的最大值。最大值的位置被指定为D_LAT并且其量值被指定为C_MAX。

必要偏移由D_LAT和其理想位置55[(640-530)/2＝55]之差确定。因而，例如如果D_LAT＝63，则必须向上移位8个像素(63-55＝8)，而如果D_LAT＝41，则需要向下移位14个像素(55-41＝14)。

为数据集中的所有图像重复步骤114。注意到每个图像参考于相同的模板，因此没有累积误差。为了帮助轴向校正，为每个图像存储29D_LAT。

现在参考图3，其示出了利用图1中所述的实施例的成像系统的轴向校正部分的功能框图。在除了第一图像之外的所有图像上执行轴向校正115。输入图像的拷贝在步骤31被加阈，然后与先前图像的加阈拷贝互相关。沿着与当前透视的侧向校正[D_LAT(N)]和前一透视的侧向校正[D_LAT(N-1)]之差对应的线，偏移被确定为互相关函数中的最大值。所以，与侧向校正114不同，轴向校正115是迭代过程，因而会有累积误差。

输入图像的拷贝以与侧向校正相同的方式被加阈31，但是在该情况下直方图中的bins数目是B3。在本例子中，B3＝2。因而，灰度级大于中范围灰度级的所有像素被设置为零，而具有更低灰度级的那些像素保持它们的初始值。例如，分别具有最小和最大值31和190的输入图像将导致与初始图像相同的加阈图像，区别在于初始亮于110的所有像素现在为零。

这样涂黑亮像素之后，加阈图像在2D中进行傅立叶变换32。然后它被滤波33以消除可能在互相关上产生伪峰值的最小特征。对应于十个像素或以下的特征尺寸的仅仅高达102周期/像素的空间频率被相乘，并且在更高空间频率处的像素被设置为零。所得到阵列被存储34为S_N并且与从先前图像的加阈拷贝获得的S_N-1的复共轭相乘35。接着对所得到阵列应用2D FFT以寻找两个连续加阈低通滤波图像的互相关。从侧向校正步骤114找到的两个连续图像之间的侧向偏移之差[D_LAT(N)-D_LAT(N-1)]现在是必需的，因为寻找校正阵列的全局最大值是错误的。相反，必须在对应于[D_LAT(N)-D_LAT(N-1)]的行中寻找局部最大值F_MAX。包含F_MAX的列被指定为G_MAX。如果G_MAX大于填补图像尺寸的一半(在该例子中为1024)，则它的值表示相对于先前图像的负位移，具有的量值等于零填补尺寸减去G_MAX的值。如果G_MAX小于零填补尺寸的一半，则相对于先前图像的所需位移是正的并且等于G_MAX。

作为例子，假设D_LAT(N-1)＝45，而D_LAT(N)＝63。那么将在相关阵列的行18上找到F_MAX(63-45＝18)。如果行18的最大值F_MAX出现在第五列中，则G_MAX＝5并且图像必须向先前图像的左边移位5个像素。如果最大值出现在行1019(G_MAX＝1019)，则图像必须向先前图像的右边移位5个像素(1024-1019＝5)，因为1019大于512。

在找到G_MAX之后，位移值与所有先前轴向偏移之和相加以确定D_AXIAL，从第一采集图像到当前图像的累积差。位移可以是正的或负的；因此，一些图像可能不需要轴向位移。对于每个图像而言，四个值被写入文本文件：

1.上管壁的位置D_LAT；

2.C_MAX，当前图像和参考图像之间的互相关的最大值；

3.G_MAX，在当前图像和先前图像之间的互相关的适当行上的F_MAX的位置(对于第一图像，G_MAX＝0)；

4.F_MAX(对于第一图像，F_MAX＝0)。

通过从一个或两个边缘裁切适当数目的像素和使剩余像素移位被裁切的数目生成校正文件。为了保持原始图像维度(640×480)，在与被裁切边缘相对的边缘的空间处由被设置成原始图像的最大灰度级的像素取代。

例如，假设对于其中一个图像，最大灰度级为229，D_LAT＝63，D_AXIAL＝29，并且G_MAX＝1022。那么顶上的八行(63-55＝8)和左边的27列(29-1024+1022＝27)从图像上被删除。因而第28列的第九行占据左上角。八行被加入到图像的底部，并且28列被加入到图像的右边；这些像素具有灰度级229。当这些步骤完成时，校正图像的左上方中的632×453像素区域与原始图像的右下方中的632×453像素区域相同。两个图像都具有640×480的尺寸。

另一个实施例仅仅包含轴向校正115和校正图像的写入116。当微毛细管壁不可见并且已知管的侧向运动小到可忽略时该实施方式是有用的。

在又一个实施方式中，从第一视图(N＝0)自动计算管壁分离。在其他方面它与上文参考图1-3所述的实施方式相同。在本发明的另一个实施方式中，基于一个或多个图像中它们分离的计算确定管壁的分离。这通过使用从单亮线通过2D FFT和复共轭得到的图像作为参考来实现。与第一实施方式相同，行被求和，并且取最大值的位置作为管的一个壁的位置，相对于参考图像从其得到的图像中它的位置而言。下一个最高相关值的位置给出相对于第一个壁的另一个管壁的位置。如果需要的话，搜索该第二最大值可以被限制于一定的范围，相对于第一管壁的该范围的中心位置在假定的管宽度附近。该实施方式也包含对于所有采集的视点使用单线参考的可能性。当管壁分离未知时，或者当管的内壁并不形成圆时，当管的内轮廓是方形或椭圆形时，这种布置可能是有用的。

在本发明的另一个实施方式中，加阈步骤的特征可以根据来自相关的反馈而变化。这样的迭代方法可以在用于侧向校正的第一加阈步骤中，在用于轴向校正的第二加阈步骤中，或同时在两者中被使用。可以变化的一个特征是用在直方图中的划分(division)或bins的数目。可以变化的另一特征是包含在每个直方图bin中的灰度级的数目。例如，直方图可以基于亮度级的平方根。

根据本发明的特征，所述方法的输出是输入文件的裁切拷贝，未裁切部分垂直地和/或水平地移位，并且在一个或两个边界处插入附加的空白像素以保持输入图像尺寸。

根据本发明的进一步特征，所利用方法的结果被存储到数字文件，该数字文件可以使用计算机字处理应用程序进行改变和编辑。然后，改变的文本文件可以用于生成两个轴上的偏移，因而回避了上述的许多计算。在该实施方式中，步骤114-116的侧向校正过程被重复以寻找最大值C_MAX。如果C_MAX具有小于临界值C_CRIT的量值，则整个过程被重复，从加阈27开始，但是直方图中的bins数目从B1变为B2。再次定位26和评估28C_MAX。

现在参考图5，其示意性地示出了本发明的实施方式所预期的填塞在毛细管中的细胞的示例。在该实施例中，毛细管3的一部分充满严密填塞在管中的感兴趣对象1，例如细胞。每个细胞可能包括细胞核2。毛细管3具有参考坐标系6定向的中心轴4，所述坐标系具有在x、y和z方向上的坐标。在一些情况下，至少一个分子探针153可以被束缚在细胞内。计算机7被连接以向旋转电机5和平移电机8提供控制信号。应当认识到一个或多个电机、齿轮或流控装置或产生运动的其他装置的等效布置也可以用于实现毛细管或其他基片必要的平移或旋转运动。在一些情况下，一个或多个电机可以由手动定位装置或齿轮或由产生运动的其他装置例如液压或压电装置代替。平移轴是z轴，旋转围绕z轴。定位电机9被连接，从而必要时为了居中而在x、y轴限定的平面内移动细胞，所述平面基本上垂直于中心轴。

应当认识到毛细管的曲面将充当柱面透镜并且该聚焦效应在投影系统中可能是非理想的。本领域技术人员应当理解当点源和管之间以及管和探测器表面之间的空间充满折射率与毛细管匹配的材料10时可以消除由管导致的光子弯曲，并且管可以光耦合(例如用油或凝胶)到空间填充材料。

考虑填塞在毛细管中的细胞的该例子。细胞优选地可以单行地被填塞，从而它们不会重叠。将直径大约为100微米的完整细胞填塞在直径小于100微米的毛细管中的密度可以大致是每厘米管长度100个细胞。对于直径大约20微米的裸原子核而言，填塞可以大致是每厘米管长度500个细胞核，其中管直径与对象尺寸成比例，在该情况下对象尺寸大约为20微米。因而，在若干厘米的毛细管长度内，可以填塞几千个不重叠的裸细胞核。通过沿着其中心轴4平移管，可以实现沿z方向的运动。必要时，沿x、y方向移动管使管内的对象在光学断层摄影系统的重建圆柱中居中．通过围绕其中心轴4旋转管，可以产生多个径向的投影视图。沿z方向以恒定速度不旋转地移动管模拟了流动光学断层摄影的特殊情况。

移动充满细胞的管但是细胞在管内不动的一个优点在于在允许在逐细胞的基础上近似最佳暴露于光学断层摄影的速度下可以停止，然后旋转感兴趣对象。也就是说，可以提高投影图像的信噪比以比通常在恒速和恒向的典型流动系统中产生更佳的图像。不感兴趣的对象可以很快被移动到成像系统之外，以便在由多个细胞组成的样本中获得分析感兴趣细胞的总速度。另外，多重投影所需的在感兴趣对象上停止然后旋转的能力几乎消除了运动伪影。更进一步地，运动系统可以被引导成以亚微米运动并且以有利地一种方式被应用，该方式允许在比探测器的像素尺寸所提供的分辨率更精细的分辨率下采样细胞。更特别地，例如以填充一半像素宽度的增量移动的运动系统可以使用值为2的尼奎斯特采样因数。类似地，运动系统可以补偿探测器的非理想填充因数。

现在参考图6A，其显示了作为例如单细胞的单样本浸没在光学分度(indexing)材料的介质中的近视图。单样本在微毛细管3(例如，一只这种管由美国LLC.，AZ的Polymicro Technologies公司制造)内被显示，其可以被旋转用于获取多重投影，并且示意性地示出了可以被轴向扫描的物镜40。照射源包括通过孔径51、挡板52和通过位于显微镜载玻片54之前的聚光透镜53投射光的光源50。微毛细管3将细胞1保持在载玻片和薄盖玻片55之间。优选为油浸透镜的物镜40被布置成接收通过微毛细管3的光。物镜沿着光轴由诸如压电元件这样的致动器57平移。盖玻片55必须足够薄以使微毛细管的中心和盖玻片的外表面之间的距离小于物镜的工作距离。聚光透镜53在折射率n₁(例如空气)内。载玻片54和盖玻片55具有折射率n₂。围绕微毛细管3的区域58包含折射率为n₃的指数匹配介质15，例如光学凝胶或浸油。微毛细管3自身具有折射率n₄。在微毛细管内围绕细胞1的区域59包含折射率为n₅的介质10。细胞内的区域60可以充满相同的介质10，或者其折射率可以不同于n₆。优选在由载玻片54和盖玻片55形成的两个平行平表面之间n₃＝n₄＝n₅＝n₆(必须使差值最小化)以避免柱面透镜失真。图像被投影到照相机43上。

现在参考图6A和6B，其示意性地示出了用在本发明中的光学显微镜系统的一个实施方式，其包含安装在压电平移装置上的显微镜物镜。压电换能器57用于使物镜60移动大约40微米或以上的轴向距离。在一个有用的实施方式中，微对象定位系统提供合适的致动器57，该致动器沿着管坐标系6的z轴被上下驱动。在该实施方式中，它可以与高数值孔径物镜和计算机控制光源和聚光透镜组件61一起使用，所述高数值孔径物镜安装在连接有摄影机43的标准传输显微镜64上。计算机控制聚光器和光源50可以有利地是包括一个或多个白炽灯、弧光灯、激光器、或发光二极管的光源。计算机控制信号70链接到计算机控制聚光器和光源50用于控制光调制。

来自照相机43的输出被存储在计算机存储器72中。包含样本的微毛细管3可以沿着管坐标系6的x或y轴平移。另外，微毛细管3可以通过能够受计算机控制的旋转电机5围绕其“θ”轴49旋转。当在这里使用时微毛细管被定义为具有一个直径的毛细管，其中用于显微镜成像的视场相当于毛细管直径。在一个实施例中旋转电机5由计算机7提供的控制信号71进行控制。为了高速应用可以加入其他控制以便减小轴向扫描期间的振动。采集的图像可以显示在监视器73上。

现在参考图7，其示出了在本发明的实施方式所预期的三维(3D)图像重建中使用的图像采集过程的流程图示例。如本发明的一个例子所预期的，3D图像重建过程包括以下步骤：在步骤81装载用细胞填塞的管，在步骤82平移管直到第一感兴趣细胞被定位，必要时在步骤83使感兴趣细胞居中，在步骤84在每个不同的旋转角处生成投影集，在步骤85确定数据集何时完成，以及重复从步骤82到85的过程直到所有感兴趣细胞都被扫描完。在步骤86进行运动校正。所述过程可以在由诸如计算机7这样的个人一计算机运行的计算机软件程序中实现。

现在参考图8，其示意性地示出了用在本发明的实施方式所预期的三维(3D)图像重建中的运动校正偏移的例子。应用运动校正以寻找包含在具有毛细管壁62的毛细管3中的感兴趣对象1(例如细胞或细胞核2)的侧向位置。侧向偏移是沿着垂直于管轴Z的较长图像维度(640像素)的误差。轴向偏移是沿着平行于管轴Z的较短图像维度(480像素)的误差。感兴趣对象1具有侧向位置LP和轴向位置AP。从多种视点采集图像，应用运动校正以便允许在多个视图中重建感兴趣对象，其具有包含在相同平面内的相同特征。

软件代码描述的例子

源代码

下而是执行用于运动校正的本发明方法一个实施方式的源代码例子的文本。可执行文件(regcorr.exe)由两个头文件(rrahn.h和regcorr.h)和由目标文件regcorr.ide链接的四个C++文件(fileinfo.cpp，nr-fft.cpp，multiplycomplexvalues.cpp，和regcorr2.cpp)构成。目标文件由为“Win32控制台”模式而设的Borland C++5.01编译，并且使用静态类库框架。可执行文件大小大致为60kB。(程序，“regcorr.exe”和“regorr_no_lateral.exe”的不同之处仅在于它们是否提供计算侧向偏移和将其写入到文本文件中的选项。这两个程序也可以跳过计算而改为从文本文件读取偏移。)

A.1 rrahn.h(38行)

#include<math.h>

#include<malloc.h>

#include<fstream.h>

#include<dos.h>

#include<iostream.h>

#include<stdio.h>

#include<conio.h>

#include<complex.h>

#pragma hdrstop

#pragma package(smart_init)

#define pi_const 3.14159265358979323846

#define unsigned int

template <class T>

inline const T SQR(const T a) {return a*a}

inline float pow(float x，double y) {return pow(double(x)，y)；}

inline float pow(double x， floaty) {return pow(x，double(y))；}

inline float atan2(float x，double y){return atan2(double(x)，y)；}

inline float atan2(double x，flaot y){return atan2(x，double(y))；}

template <class T>

inline void SWAP(T &a，T &b)

{T dum＝a；a＝b；b＝dum；}

#define SwapFFT SWAP

#define wait while(！kbhit())；return；

unsigned *UIivector(unsigned，unsigned)；

A.2 regcorr.h(40行)

#include″rrahn.h″

#pragma hdrstop

#pragma pakage(smart_init)

const unsigned ImageDimX＝640；

const unsigned ImageDimZ＝480；

const unsigned BigDimX＝1024；

const unsigned BigDimZ＝1024；

const unsigned NumHistoBins＝20；

const unsigned WallSpacingMin＝505，WallSpacingMax＝565；

extern unsigned FileYear， FileMonth， FileDay，FileSet，NumPerspectives，

CurrentPerspectve， BinHist[NumHistoBins]；

extern int SliceNumber；

extern char filenameIn[64]；

extern char ProcessedFilenameOut[64]；

extern float RawImage[ImageDimX*ImageDimZ]；

extern double ShifedImage[BigDimX*BigDimZ*2]，

ShifedImageNew[BigDimX*BigDimZ*2]；

extern double TemplateImage[BigDimX*BigDimZ*2]，

TemplateImageNew[BigDimX*BigDimZ*2]；

extern int ShiftX；

extern int ShiftZ；

extern double MaxVal，MinVal；

extern float BinVal[ImageDimX*ImageDimZ]；

extern unsigned MaxBin， MaxHisto；

extern int BigDims[2]；

extern char* DirectoryName；

extern unsigned WallEdge1，WallEdge2，WallSpacing；

extern double MinValDev1，MinValDev2，MaxValCorr，MaxValCorr2；

extern unsigned MinDevIndexOld；

extern double ShiftedMag[2*BigDimX*BigDimZ]；

extern unsigned ppmax；

extern float CriticalValue；

void MultiDimFFt(double[]，int[]，int，int)；

void FindDeviation()；

void GetFileInfo()；

void WriteCorrectedImage(double[])；

void MultiplyComplexValues(double[]，double[]，double[]，unsigned，unsigned，

unsigned，unsigned)；

void BasicFFT(float[]，unsigned long，int)；

void MakeFileName(char*，unsigned，unsigned)；

A.3 fileinfo.cpp(50行)

#include″rrahn.h″

unsigned NumPerspectives；

unsigned FileYear，FileMonth，fileDay，FileSet；

void GetFileInfo()

{

cout＜＜″Year(YY，3 or 4)：″；

cin＞＞FileYear；

while(FileYear！＝3&&FileYear！＝4)

{

cout＜＜″Year(YY，3 or 4)：″；

cin＞＞FileMonth；

while(FileMonth＞12)

{

cout＜＜″Month(MM)：″；

cin＞＞FileMonth；

}

cout＜＜end1＜＜″Day(DD)：″；

cin＞＞fileDay；

while(FileDay＞31)

{

cout＜＜″Day(DD)：″；

cin＞＞FileDay；

}

cout＜＜end1＜＜″Set(1-99)：″；

cin＞＞FileSet；

while(FileSet＞99)

{

cout＜＜″Set(1-99)：″；

cin＞FileSet；

}

cout＜＜end1；

cout＜＜″Number of Perspectives(1-255)：″；

cin＞＞NumPerspectives；

while(NumPerspectives＞255||NumPerspectives＜1)

{

cout＜＜″Number of Perspectives(1-255)：″；

cin＞＞NumPerspectives；

}

return；

}

A.4 nr-fft.cpp(221行)

#include″rrahn.h″

void BasicFFT(float FFTdata[]，unsigned long FFTnn，int FFTsign)

{

unsigned long FFTm，FFTistep；

unsigned long FFTmmax＝2；

double FFTwtemp，FFTwr，FFTwpr，FFTwpi，FFTwi，FFTtheta；

float FFTtempr，FFTtempi；

unsigned long FFTn＝FFTnn＜＜1；

unsigned long FFTj＝1；

for(unsigned long jj＝1；jj＜FFTn；jj+＝2)

if(FFTj＞jj)

{

SwapFFT(FFTdata[FFTj1]，FFTdata[jj])；

SwapFFT(FFTdata[FFTj+1]，FFTdata[jj+1])；

}

FFTm＝FFTn＞＞1；

while(FFTm＞＝2&&FFTj＞FFTm)

{

FFTj-＝FFTm；

FFTm＝FFTm＞＞1；

}

FFTj+＝FFTm；

}

while(FFTn＞FFTmmax) {

FFTistep＝FFTmmax＜＜1；

FFTtheta＝FFTisign*(2.0*pi_const/FFTmmax)；

FFTwtemp＝sin(0.5*FFTtheta)；

FFTwpr＝-2.0*FFTwtemp*FFTwtemp；

FFTwpi＝sin(FFTtheta)；

FFTwr＝1.0；

FFTwi＝0.0；

for(FFTm＝1；FFTm＜FFTmmax；FFTm+＝2)

{

for(unsigned long jj＝FFTm；jj＜＝FFTn；jj+＝FFTistep)

{

FFTj＝jj+FFTmmax；

FFTtempr＝FFTwr*FFTdata[FFTj]-FFTwi*FFTdata[FFTj+1]；

FFTtempi＝FFTwr*FFTdata[FFTj+1]+FFTwi*FFTdata[FFTj]；

FFTdata[FFTj+1]＝FFTdata[jj+1]-FFTtempi；

FFTdata[jj]+＝FFTtempr；

FFTdata[jj+1]+＝FFTempi；

}

FFTwtemp＝FFTwr；

FFTwr＝FFTwtemp*FFTwpr-FFTwi*FFTwpi+FFTwr；

FFTwi＝FFTwi*FFTwpr+FFTwtemp*FFTwpi+FFTwi；

}

FFTmmax＝FFTistep；

}

return；

}

void ApplyOneDimFFT(float FFTdata[]，unsigned long FFTn，int FFTsign)

{

unsigned long FFTi，FFTi1，FFTi2，FFTi3，FFTi4，FFTnp3；

float FFTc1＝0.5，FFTc2，FFThlr，FFThli，FFTh2r，FFTh2I；

double FFTwr，FFTwi，FFTwpr，FFTwpi，FFTwtemp，FFTtheta；

FFTtheta＝pi_const/(double)(FFTn＞＞1)；

void BasicFFT(float FFTdata[]，unsigned long FFTnn，int FFTisign)；

if(FFTisign＝＝1)

{

FFTc2＝-0.5；

BasicFFT(FFTdata，FFTn＞＞1，1)； }

else

{

FFTc2＝0.5；

FFTtheta＝-FFTtheta；

}

FFTwtemp＝sin(0.5*FFTtheta)；

FFTwpr＝-2.0*FFTwtemp*FFTwtemp；

FFTwpi＝1.0+FFTwpr；

FFTwi＝FFTwpi；

FFTnp3＝FFTn+3；

For(FFTi＝2；FFT＜＝(FFTn＞＞2)；FFTi++)

{

FFTi1＝FFTi+FFTi-1；

FFTi2＝1+FFTi1；

FFTi3＝FFTnp3-FFTi2；

FFTi4＝1+FFi3；

FFTh1r＝FFTc1*(FFTdata[FFTi1]+FFTdata[FFTi3])；

FFTh1i＝FFTc1*(FFTdata[FFTi2]-FFTdata[FFTi4])；

FFTh2r＝-FFTc2*(FFTdata[FFTi2]+FFTdata[FFTi4])；

FFT2i＝FFTc2*(FFTdata[FFTi1]-FFTdata[FFTi3])；

FFTdata[FFTi1]＝FFTh1r+FFTwr*FFTh2r-FFTwi*FFTh2i；

FFTdata[FFTi2]＝FFTh1i+FFTwr*FFTh2i+FFTwi*FFTh2r；

FFTdata[FFTi3]＝FFTh1r+FFTwr*FFTh2r-FFTwi*FFTh2i；

FFTdata[FFTi4]＝-FFTh1i+FFTwr*FFTh2i+FFTwi*FFTh2r；

FFTwtemp＝FFTwr；

FFTwr＝FFTwtemp*FFTwpr-FFTwi*FFTwpi+FFTwr； //递推

FFTwi＝FFTwi*FFTwpr+FFTwtemp*FFT+FFTwi；

}

if(FFTisign＝＝1)

}

FFTh1r＝FFTdata[1]；

FFTdata[1]＝FFTh1r+FFTdata[2]；

FFTdata[2]＝FFTh1r-FFTdata[2]；

}

else

{

FFTh1r＝FFTdata[1]；

FFTdata[1]＝FFTc1*(FFTh1r+FFTdata[2])；

FFTdata[2]＝FFTc1*(FFTh1r-FFTdata[2])；

BasicFFT(FFTdata，FFTn＞＞1，-1)；

For(unsignedjj＝0；jj＜FFTn；jj++)

{

FFTdata[jj+1]＝2.0*FFTdata[jj+1]/FFTn；

}

void MultiDimFFT(double data[]，int nn[]，int ndim，int isign)

{

int i1，i2，i3，i2rev，i3rev，ip1，ip2，ip3，ifp1，ifp2；

int ibit，idim，k1，k2，n，nprev，nrem，ntot；

double tempi，tempr；

double theta，wi，wpi，wpr，wr，wtemp；

ntot＝1；

for(idim＝1；idim＜＝ndim；idim++) ntot*＝nn[idim]；

nprev＝1；

for(idim＝ndim；idim＞＝1；idim--) {

n＝nn[idim]；

nrem＝ntot/(n*nprev)；

ip1＝nprev＜＜1；

ip2＝ip1*n；

ip3＝ip2*nrem；

i2rev＝1；

for(i2＝1；i2＜＝ip2；i2+＝ip1)

{

if(i2＜i2rev)

{

for(i1＝2；i1＜＝(i2+ip1-2)；i1+＝2)

{

for(i3＝i1；i3＜＝ip3；i3+＝ip2)

{

i3rev＝i2rev+i3-i2；

SwapFFT(data[i3]，data[i3rev])；

SwapFFT(data[i3+1]，data[i3rev+1])；

}

ibit＝ip2＞＞1；

while((ibit＞＝ip1)&&(i2rev＞ibit))

{

i2rev-＝ibit；

ibit＞＞＝1；

}

i2rev+＝ibit；

}

ifp1＝ip1； while(ifp1＜ip2)

{

ifp2＝ifp1＜＜1；

theta＝isign*2*pi_const/(ifp2/ip1)；

wtemp＝sin(0.5*theta)；

wpr＝-2.0*SQR(wtemp)；

wpi＝sin(theta)；

wr＝1.0；

wi＝0.0；

for(i3＝1；i3＜＝ifp1；i3+＝ip1)

{

for(i1＝i3；i1＜＝(i3+ip1-2)；i1+＝2)

{

for(i2＝i1；i2＜＝ip3；i2+＝ifp2)

{

k1＝i2；

k2＝k1+ifp1；

tempr＝wr*data[k2]-wi*data[k2+1]；

tempi＝wr*data[k2+1]+wi*data[k2]；

data[k2]＝data[k1]-tempr；

data[k2+1]＝data[k1+1]-tempi；

data[k1]+＝tempr；

data[k1+1]+＝tempi；

}

wtemp＝wr

wr＝wtemp*wpr-wi*wpi+wr；

wi＝wi*wpr+wtemp*wpi+wi；

}

ifp1＝ifp2；

}

nprev*＝n；

}

return；

}

A.5 multiplycomplexvalues.cpp(30行)

#include″rrahn.h″

void MultiplyComplexValues(double MCV1[]，double MCV2[]，double MCV3[]，

unsigned MCVDim1，unsigned MCVDim2，unsigned MaxFrequency1，unsigned

MaxFrequency2)

{

for(unsigned mm＝0；mm＜MCVDim2；mm++)

{

for(unsigned nn＝0；nn＜MCVDim1；nn++)

{

if((mm＞＝MaxFrequency2 && mm＜MCVDim2-MaxFrequency2) | |

(nn＞＝MaxFrequency1 && nn＜MCVDim1-MaxFrequency1))

{

MCV3[2*nn*MCVDim2+2*mm]＝(double)0.0；

MCV3[2*nn*MCVDim2+2*mm+1]＝(double)0.0；

}

else

{

MCV3[2*nn*MCVDim2+2*mm] ＝ (double)

(MCV1[2*nn*MCVDim2+2*mm] *MCV2[2*nn*MCVDim2+2*mm] +

MCV1[2*nn*MCVDim2+2*mm+1]*MCV2[2*nn*MCVDim2+2*mm+1])；

MCV3[2*nn*MCVDim2+2*mm+1] ＝ (double)

(MCV2[2*nn*MCVDim2+2*mm]*MCV1[2*nn*MCVDim2+2*mm+1]-

MCV2[2*nn*MCVDim2+2*mm+1]*MCV1[2*nn*MCVDim2+2*mm])；

}

return；

}

A.6 regcorr2.cpp(374行)

#include<strstrea.h>

#include″regcorr.h″

NumPerspectives；

unsigned CurrentPerspective；

unsigned BinHist[NumHistoBins]；

int SliceNumber；

char fienameIn[64]；

char ProcessedFilenameOut[64]，IndexFilenameOut[64]；

float RawImage[ImageDimX*ImageDimZ]；

double ShiftedImage[BigDimX*BigDimZ*2]；

double TemplateImage[BigDimX*BigDimZ*2]，

TemplateImageNew[BigDimX*BigDimZ*2]；

double oneLine[BigDimX*BigDimZ*2]；

float LateralSum[2*BigDimX]；

double MaxVaL，MinVal；

float BinVal[ImageDimX*ImageDimZ]；

unsigned MaxBin，MaxHisto，MinDevIndexOld；

intBigDims[2]，SmallDims[2]；

unsigned MaxCorrIndex，WallFdge1，WallEdge2；

double MaxValCorr；

double ShiftedMag[2*BigDimX*BigDimZ]；

char*DirectoryName＝″E:\\VisionGate\\Projection Images\\″；

unsigned nnc，WallOld1，WallOld2；

float LateralMax1＝0，LateralMax2＝0；

float MaxValOld，MinValOld；

unsigned SecondSeg＝2；

unsigned FirstSeg＝20；

unsigned RedoLateralCorrection＝2.2e17；

intCumulativeShift＝0；

void ReadUncorrectedImage(unsigned RUI1)

{

MaxValOld＝MaxVal；

MinValOld＝MinVal；

MaxVal＝0；

MinVal＝1e22；

MakeFileName(filenameIn，RUI1，2)；

fstream RawFile(filenameIn，ios::in)；

Cout＜＜filenameIn＜＜end1；

for(unsigned mm＝0；mm＜ImageDimZ；mm++)

{

for(unsigned nn＝0；nn＜ImageDimX；nn++)

{

RawFile＞＞(float)RawImage[nn*ImageDimZ+mm]；

if(RawImage[nn*ImageDimZ+mm]＜MinVal)

MinVal＝RawImage[nn*ImageDimZ+mm]；

if(RawImage[nn*ImageDimZ+mm]＞MinVal)

MaxVal＝RawImage[nn*ImageDimZ+mm]；

}

RawFile.close()

return；

}

void CopyRawImage(double CRI1[])

{

for(unsigned mm＝0；mm＜BigDimZ；mm++)

for(unsigned nn＝0；nn＜BigDimX；nn++)

{

CRI1[2*nn*BigDimZ+2*mm+1]＝0;

CRI1[2*nn*BigDimZ+2*mm]＝0；

if(mm＜ImageDimZ && nn＜ImageDimX)

(double)CRI1[2*nn*BigDimZ+2*mm]＝RawImage[nn*ImageDimZ+mm]；

}

return；

}

void MakeOneLine()

{

for(unsigned jj＝0；jj＜ImageDimZ；jj++)

{

OneLine[2*jj]＝65535；

OneLine[2*jj+2*530*BigDimX]＝65535；

}

MultiDimFFT(OneLine-1，BigDims-1，2，1)；

return；

}

void MakeHistogram(double MHdata[]，unsigned MHBins)

{

if(MHBins＞NumHistoBins)

MHBins＝NumHistoBins；

MaxBin＝0；

MaxHisto＝0；

for(unsigned jj＝0；jj＜NumHistoBins；jj++)

BinHist[jj]＝0；

for(unsigned mm＝0；mm＜ImageDimZ；mm++)

{

for(unsigned nn＝0；nn＜ImageDimX；nn++)

{

BinVal[nn*ImageDimZ+mm] ＝ (int)

MHBins*(RawImage[nn*ImageDimZ+mm]-MinVal)/(MaxVal-MinVal)；

if(BinVal[nn*ImageDimZ+mm]＞＝MHBins)；

BinVal[nn*ImageDimZ+mm]＝MHBins-1u；

if(BinVal[nn*ImageDimZ+mm]＜＝0)；

BinVal[nn*ImageDimZ+mm]＝0u；

BinHist[(int)BinVal[nn*ImageDimZ+mm]]+＝1u；

}

for(unsigned jj＝0；jj＜MHBins；jj++)

{

if(BinHist[jj]＞MaxHisto)

{

MaxHisto＝BinHist[jj]；

MaxBin＝jj；

}

for(unsigned jj＝0；jj＜2*BigDimX*BigDimZ；jj++)

Mhdata[jj]＝0；

For(unsigned mm＝0；mm＜ImageDimZ；mm++)

{

for(unsigned nn＝0；nn＜ImageDimX；nn++)

{

if(BinVal[nn*ImageDimZ+mm]＜MaxBin)

MHdata[2*nn*BigDimZ+2*mm]＝RawImage[nn*ImageDimZ+mm]；

}

return；

}

void LateralCorrection(unsigned LC1)

{

MakeHistogram(ShiftedImage，LC1)；

MultiDimFFT(ShiiftedImage-1，BigDims-1，2，1)；

MultiplyComplexValues(ShiftedImage，OneLine，TemplateImageNew，BigDimX，

BigdimZ，BigDimX，BigDimZ)；

for(unsigned mm＝0；mm＜BigDimX；mm++)

{

LateralSum[2*mm]＝0；

LateralSum[2*mm+1]＝0；

for(unsigned nn＝0；nn＜BigDimZ；nn++)

{

LateralSum[2*mm]+＝TemplateImageNew[2*mm*BigDimZ+2*nn]；

LateralSum[2*mm+1]+＝TemplateImageNew[2*mm*BigDimZ+2*nn+1]；

}

BasicFFT(LateralSum-1，BigDimX，-1)；

LateralMax1＝0；

for(utnsigned jj＝0；jj＜110；jj++)

{

if(LateralSum[2*jj]＞LateralMax1)

{

WallEEdge1＝jj；

LateralMax1＝LateralSum[2*jj]；

}

WallEdge2＝WallEdge1+530；

return；

}

void AxialCorrection(unsigned AC1)

{

MaxValCorr＝0.0；

MaxCorrIndex＝0；

MakeHistogram(ShiftedImage，AC1)；

for(unsigned jj＝0；jj＜2*BigDimZ*BigDimX；jj++)

{

if(ShiftedImage[jj]＞0)

ShiftedImage[jj]＝MaxVal+MinVal-ShiftedImage[jj]；

if(TemplateImage[jj]＞0)

TemplateImage[jj]＝MaxValOld+MinValOld-TemplateImage[jj]；

}

MultiDimFFT(TempiateImage-1，BigDims-1，2，1)；

MultiDimFFT(ShiftedImage-1，BigDims-1，2，1)；

MultiplyComplexValues(ShiftedImage， TemplateImage， TemplateImageNew，

BigDimX，BigDimZ，BigDimX/10，BigDimZ/10)；

MultiDimFFT(TemplateImage-1，BigDims-1，2，-1)；

if(MinDevINdexOld＞(WallEdge1+WallEdge2)/2)

nnc＝BigDimX+(WallEdge1+WallEdge2)/2-MinDevIndexOld；

else

nnc＝(WallEdge1+WallEdge2)/2-MinDevIndexOld；

for

{

if(TemplateImageNew[2*nnc*BigDimZ+2*mm]＞MaxValCorr)

{

MaxValCorr＝TemplateImageNew[2*nnc*BigDimZ+2*mm]；

MaxCorrIndex＝mm；

}

return；

}

void CrossCorrelation()

{

//如果需要侧向校正则包括这两行

// WallEdge1＝55；

// WallEdge2＝WallEdge1+530；

MinDevIndexOld＝(WallEdge1+WallEdge2)/2；

WallOld1＝WallEdge1；

WallOld2＝WallEdge2；

//如果不需要侧向校正则删除这三行

LateralCorrection(FirstSeg)；

if(LateralMax1＜RedoLateralCorrection)

LateralCorrection(SecondSeg)；

AxialCorrection(SecondSeg)；

MakeHistogram(TemplateImage，SecondSeg)；

return；

}

void MakeFileName(char*MFN1，unsigned MFN2，unsigned MFN3)

{

if(MFN3＝＝0)

sprintf(MFN1，

″％spp％2.2d％2.2d％2.2d％2.2d\\pp％2.2d％2.2d％2.2d％2.2dindex.txt″，

DirectoryName， FileYear， FileMonth，FileDaay，FileSet，FileYear， FileMonth，

FileDay，FileSet)；

if(MFN3＝＝1)

sprintf(MFN1，

″％spp％2.2d％2.2d％2.2d％2.2d\\pp％2.2d％2.2d％2.2d％2.2％3.3dAVG.crw″，

DirectoryName， FileYear， FileMonth， FileDay，FileSet，FileYear，FileMonth，

FileDay，FileSet，MFN2)；

if(MFN3＝＝2)

sprintf(MFN1，

″％spp％2.2d％2.2d％2.2d％2.2d\\pp％2.2d％2.2d％2.2d％2.2％3.3dAVG.raw″，

DirectoryName， FileYear， FileMonth，FileDay， FileSet，FileYear，FileMonth，

FileDay，FileSet，MFN2)；

return；

}

void WriteCorrectedImage(float WCIdata[])

{

int XCenter＝(WallEdge1+WallEdge2)/2；

if(XCenter＞＝ImageDimX/2)

{

if(CumulativeShift＜＝0)

{

for(unsigned jj＝0；jj＜ImageDimX-XCenter+ImageDimX/2；jj++)

{

for(int kk＝ImageDimZ-1；kk＞＝-CumulativeShift；kk--)

// for(unsigned kk＝-CumulativeShift；kk＜ImageDimz；kk++)

WCIdata[jj*ImageDimZ+kk] ＝

WCIdata[(jj+XCenter-ImageDimX/2)*ImageDimZ+kk+CumulativeShift]；

for(unsigned kk＝0；kk＜-CumulativeShift；kk++)

WCIdata[jj*ImageDimZ+kk]＝MaxVal；//0；

}

for(unsigned jj＝(ImageDimX-XCenter+ImageDimX/2)*ImageDimZ；jj＜

ImageDimX*ImageDimZ；jj++)

WCIdata[jj]＝MaxVal；//0；

}

if(CumulativeShift＞0)

{

for(unsigned jj＝0；jj＜ImageDimX-XCenter+ImageDimX/2；jj++)

{

for(int kk＝0；kk＜ImageDimZ-CumulativeShift；kk++)

WCIdata[jj*ImageDimZ+kk]＝WCIdata[(jj+XCenter-ImageDimX/2)*I

mageDimZ+kk+CumulativeShift]；

for(unsigned kk＝ImageDimZ-CumulativeShift；kk＜ImageDimZ；kk++)

WCIdata[jj*ImageDimZ+kk]＝MaxVal；//0；

}

for(unsigned jj＝(ImageDimX-XCenter+ImageDimX/2)*ImageDimZ；

jj＜ImageDimX*ImageDimZ；jj++)

WCIdata[jj]＝MaxVal；//0；

}

if(XCenter＜ImageDimX/2)

{

if(CumulativeShift＜＝0)

{

for(int jj＝ImageDimX-1；jj＞＝(-XCenter+ImageDimX/2)；jj--)

{

for(int kk＝ImageDimZ-1；kk＞＝-CumulativeShift；kk--)

// for(unsigned kk＝-CumulativeShift；kk＜ImageDimZ；kk+)

WCIdata[jj*ImageDimZ+kk]＝WCIdata[(jj+XCenter-ImageDimX/

2)*ImageDimZ+kk+CumulativeShift]；

for(unsigned kk＝0；kk＜-CumulativeShift；kk++)

WCIdata[jj*ImageDimZ+kk]＝MaxVal；//0；

}

for(unsigned jj＝0；jj＜(-XCenter+ImageDimX/2)*ImageDimZ；jj++)

WCIdata[jj]＝MaxVal；//0；

}

if(CumulativeShift＞0)

{

for(int jj＝ImageDimX-1；jj＞＝(-XCenter+ImageDimX/2)；jj--)

{

for(unsigned kk＝0；kk＜ImageDimZ-CumulativeShift；kk++)

// for(int kk＝ImageDimZ-CumulativeShift-1；kk＞＝0；kk--)

WCIdata[jj*ImageDimZ+kk]＝WCIdata[(jj+XCenter-ImageDimX/2)*I

mageDimZ+kk+CumulativeShift]；

for(unsigned kk＝ImageDimZ-CumulativeShift；kk＜ImageDimZ；kk++)

WCIdata[jj*ImageDimZ+kk]＝MaxVal；//0；

}

for(unsigned jj＝0；jj＜(-XCenter+ImageDimX/2)*ImageDimZ；jj++)

WCIdata[jj]＝MaxVal；//0；

}

MakeFileName(ProcessedFilenameOut，CurrentPerspective，1)；

fstream ProcessedFile；

ProcessedFile.open(ProcessedFilenameOut，ios::binary|ios::out)；

for(unsigned jj＝0；jj＜ImageDimZ；jj++)

{

for(unsigned kk＝0；kk＜ImageDimX；kk++)

ProcessedFile＜＜WCIdata[kk*ImageDimZ+jj]＜＜″\t″；

}

ProcessedFile.close()

return；

}

void WriteToIndex File()

{

MakeFileName(IndexFilenameOut，0，0)；

fstream IndexFile；

IndexFile.open(IndexFilenameOut，ios::binary|ios::app)；

IndexFile＜＜WallEdge 1＜＜″\t″＜＜LateralMax1＜＜″\t″＜＜MaxCorrIndex＜＜″\t″

＜＜MaxValCorr＜＜end1；

IndexFile.close()；

return；

}

void CalculateOffsets()

{

MakeFileName(IndexFilenameOut，0，0)；

fstream IndexFile；

IndexFile.open(IndexFilenameOut，ios::binary|ios::out)；

IndexFile.close()；

BigDims[0]＝BigDimX；

BigDims[1]＝BigDimZ；

MakeOneline()；

for(CurrentPerspective＝0； CurrentPerspective＜＜NumPerspectives；

CurrentPerspective++)

{

ReadUncorrectedImage(CurrentPerspective)；

CopyRawImage(ShiftedImage)；

CrossCorrelation()；

if(MaxCorrIndex＞BigDimZ/2)

CumulativeShift-＝(BigDimZ-MaxCorrIndex)；

else

CumulativeShift+＝MaxCorrIndex；

WriteCorrectedImage(RawImage)；

WriteToIndexFile()；

}

return；

}

void ReadOffsets()

{

double Junk；

MakeFileName(IndexFilenameOut，0，0)；

fstream IndexFile；

IndexFile.open(IndexFilenameOut，ios::binary|ios::in)；

for(unsigned jj＝0；jj＜2*NumPerspectives；jj++)

{

IndexFile＞＞(unsigned)CorrectedFile[jj]；

IndexFile＞＞Junk；

Cout＜＜CorrectedFile[jj]＜＜″\n″；

}

IndexFile.close()；

for(CurrentPerspective＝0； CurrentPerspective＜NumPerspectives；

CurrentPerspective++)

{

ReadUncorrectedImage(CurrentPerspective)；

WallEdge1＝CorrectedFile[2*CurrentPerspective]；

WallEdge2＝WallEdge1+530；

MaxCorrIndex＝CorrectedFile[2*CurrentPerspective+1]；

if(MaxCorrIndex＞BigDimZ/2)

CumulativeShift-＝(BigDimZ-MaxCorrIndex)；

else

CumulativeShift+＝MaxCorrIndex；

Cout＜＜CumulativeShift＜＜end1；

WriteCorrectedImage(RawImage)；

}

return；

}

void main()

{

int Calcoffs；

GetFileInfo()；

printf(″Calculate offsets？(Y/N)\n″)；

Calcoffs＝getch()；

if(Calcoffs＝＝′Y′||Calcoffs＝＝′Y′)

CalculateOffsets()；

if(Calcoffs＝＝′N′||Calcoffs＝＝′n′)

ReadOffsets()；

cout＜＜″Done！″＜＜end1；

wait

}

在这里已经相当详细地描述了本发明以便符合专利法和为本领域的技术人员提供应用本发明的新颖原理所需的信息，以及当需要时构造和使用这种示例性和专门的部件。然而，应当理解的是本发明可以由特别不同的设备、装置和重建算法实现，并且可以实现关于设备细节和工作过程的多种改进而不背离本发明的实质精神和范围。

Claims

1.一种用于校正连续视图之间的对象一探测器相对运动的方法，其包括以下步骤：

照射感兴趣对象(1)以产生图像(111)；

确定图像的侧向偏移校正值(114)；

确定图像的轴向偏移校正值(115)；和

对图像应用侧向偏移校正值和轴向偏移校正值以产生校正的文件图像(116)。

2.根据权利要求1所述的方法，其中，确定图像的侧向偏移校正值的步骤(114)进一步包括以下步骤：

加阈图像(21)；和

互相关所述图像和模板图像(22，23，24，25，26，27)。

3.根据权利要求2所述的方法，其中，通过一种方法创建(111)模板图像，所述方法包括以下步骤：

在预定位置创建至少两条白线以形成预备模板图像；

扩展预备模板图像以在导致扩展模板图像的两个维度中提供零填补；和

在扩展模板图像上执行二维FFT(22)以创建最终模板图像。

4.根据权利要求3所述的方法，其中，加阈图像的步骤(21)进一步包括以下步骤：

寻找包括多个bins的图像的灰度直方图；

识别具有最大数目像素的bin；

在图像拷贝中将具有bin的灰度值或更高灰度值的图像中的所有像素设置成等于零；

对图像拷贝应用二维FFT(22)以产生傅立叶变换；

使傅立叶变换与模板图像的傅立叶变化的复共轭相乘以产生新的图像阵列(23)；

沿着多个行中的每一行对新图像阵列求和(24)以计算侧向和阵列；

计算侧向和阵列的一维傅立叶变换(25)以寻找原始图像拷贝和模板图像的行的互相关；

在侧向和阵列的最大值的位置设置上管壁的未校正位置(26)；和

将上管壁的未校正位置和预定壁边缘位置之差确定为侧向偏移(28)。

5.根据权利要求4所述的方法，其中，确定图像的轴向偏移校正值的步骤进一步包括以下步骤：

加阈图像(31)；

互相关所述图像和先前图像的加阈形式(32，33，34，35，36，37)；和

将沿着与两个图像的侧向校正之差对应的线的互相关函数中的最大值确定为轴向偏移(38)。

6.根据权利要求5所述的方法，其中，确定轴向偏移的步骤进一步包括以下步骤：

寻找包括多个bins的图像的灰度直方图；

识别具有最大数目像素的bin；

在原始图像拷贝中将具有bin的灰度值或以上灰度值的原始图像中的所有像素设置成等于零(31)；

对加阈图像应用低通滤波(33)；

计算加阈低通滤波图像和先前图像的加阈低通滤波形式的互相关(35，36)；

在与两个图像的侧向偏移之差对应的结果互相关的行中寻找最大相关值(37)；和

将相关值加入到所有先前轴向偏移的和中。

7.根据权利要求6所述的方法进一步包括以下步骤：

将上管壁的位置写入到文件；

将当前图像和模板图像之间的互相关的最大值写入到文件；

将互相关的适当行上的最大相关值的位置写入到文件，其中所述适当行是与两个图像的侧向偏移之差对应的行；

将最大校正值写入到文件；和

通过从一个或两个边缘裁切适当数目的像素和使剩余像素移位被裁切的数目而生成校正图像。

8.根据权利要求1所述的方法，其中，感兴趣对象(1)是细胞或细胞核。

9.根据权利要求4所述的方法，进一步包括以下步骤：当互相关的最大值具有小于预定值的量值(28)时，重复所述步骤，并且直方图中的bin数目不同。

10.用于校正成像系统中的定位和运动误差的装置，其包括：

部件(5，43，50，60)，用于采集对象(1)的至少两个视图；所述至少两个视图包括M行和N列的图像亮度数据；

部件(7)，用于对来自每个采集视图的图像亮度数据的图案和图像亮度数据的相应参考图案进行比较以产生两个系数“A”和“B”，所述系数包含关于它们之间的相似性的信息；

用于比较的所述部件包括：

第一部件，用于产生所述采集图像的修正拷贝，其中亮度级的改变被应用于具有特定范围内的亮度级的像素；

第二部件，用于执行采集视图的所述修正拷贝和所述参考视图之间的互相关；

第三部件，用于定位和评估结果互相关函数的最大值，被定位的所述最大值的行或列确定系数“A”的值；

第四部件，用于产生所述采集图像的修正拷贝，其中亮度级的改变被应用于具有特定范围内的亮度级的像素；

第五部件，用于执行连续采集图像的两个所述修正拷贝之间的互相关；

第六部件，用于定位和评估结果互相关函数的受限最大值，所述限制在于所述受限最大值所处的行或列必须与系数“A”所处的行或列正交和交叉；并且

所述受限最大值的位置确定系数“B”。

11.根据权利要求10所述的装置，进一步包括响应于所述系数“A”和“B”的所述值以用于沿两个正交方向移位采集视图的部件，所述量的方向和量值与所述系数“A”和“B”的所述值有关。

12.根据权利要求10所述的装置，其中，用于互相关以确定系数“A”的部件包括求和行或列以获得一维阵列的步骤，所述求和沿着与位移的方向正交的方向被执行，所述位移的量值由所述系数“A”确定。

13.根据权利要求10所述的装置，其中，参考图案包括M行和N列图像亮度数据的阵列，其中除了一个或多个行或列之外的所有像素具有相同的亮度值，所述一个或多个行或列具有一致的，但是实质上不同于其他像素的亮度。

14.根据权利要求10所述的装置，其中，参考图案包括M行和N列图像亮度数据的阵列，其中除了由一个或多个行或列组成的两个组之外的所有像素具有相同的亮度值，所述两个组具有一致的，但是实质上不同于其他像素的亮度，并且具有与微毛细管(3)的相对壁的表观分离对应的分离，所述相对壁都出现在连续的两个或以上图像的视场中，所述微毛细管(3)包含感兴趣对象。

15.根据权利要求14所述的装置，其中，感兴趣对象(1)是生物细胞或细胞核。

16.根据权利要求10所述的装置，其包括重复加阈、互相关以及定位和评估所述互相关的最大值的步骤的部件，所述重复包括基于多重加阈寻找最大相关的最大值的部件。

17.根据权利要求14所述的装置，其中，基于一个或多个采集图像计算两个组或行或列的分离。

18.一种用于重建感兴趣对象的三维(3D)重建方法，其包括以下步骤：

(a)将感兴趣对象的集合填塞在线性容器中(81)；

(b)用至少一个光学投影射束照射所述感兴趣对象的集合中的至少一个对象(50)；

(c)平移线性容器直到所述至少一个对象被定位(82)在所述至少一个光学投影射束的区域内；

(d)必要时使所述至少一个对象居中(83)；

(e)使所述至少一个对象转过多个径向角度(84)；

(f)在所述多个径向角度的每个径向角度处生成投影图像集(84)；

(g)重复步骤(b)到(f)直到感兴趣对象(1)的集合被扫描完(85)；和

(h)通过确定扫描图像的侧向偏移校正值，确定扫描图像的轴向偏移校正值，并对扫描图像应用侧向偏移校正值和轴向偏移校正值以产生校正的文件图像，从而校正被扫描感兴趣对象(1)的运动(86)。

19.根据权利要求18所述的方法，其中，所述至少一个对象(1)是细胞或细胞核。

20.根据权利要求18所述的方法，其中，将感兴趣对象的集合填塞在线性容器中的步骤(81)包括将多个细胞填塞在微毛细管(3)中的步骤。

21.一种用于校正连续视图之间的对象—探测器相对运动的方法，其包括以下步骤：

成像感兴趣对象(1)以产生图像(111)；

确定图像的轴向偏移校正值(115)，其中确定图像的轴向偏移校正值的步骤进一步包括以下步骤：加阈图像(31)，互相关所述图像和先前图像的加阈形式(32，33，35，36)，并将轴向偏移确定为沿着与两个图像的侧向校正之差对应的线的互相关函数中的最大值(37)；和

确定图像的侧向偏移校正值(114)；

22.根据权利要求1所述的方法，其中，产生多个图像以生成输入图像文件并且多个校正的文件图像是输入图像文件的裁切拷贝，其中未裁切部分垂直地和/或水平地移位，并且在一个或两个边界处插入附加的空白像素以保持输入文件图像尺寸。

23.根据权利要求1所述的方法，其中，校正文件图像包括被存储到数字文件的计算出的偏移值。

24.根据权利要求23所述的方法，其中，使用计算机字处理来处理计算出的偏移值以产生更改文本文件，而更改文本文件用于生成两个轴的偏移，并且其中用于确定图像的侧向偏移校正值的步骤被重复以寻找最大值。

25.根据权利要求22所述的方法，其中，计算出的偏移值包括：

上管壁的位置；

当前图像和模板图像之间的互相关的最大值；

当前图像和先前图像之间的互相关的相应行上的最大相关值的位置；和

最大相关值。