CN101449590A - 用于去除光学断面图像线条的系统和方法 - Google Patents

Abstract

披露了一种用于产生图像的设备、系统和方法。处理器可以基于多幅输入图像产生第一输出图像，以及如果有任何伪影，从该第一输出图像除去伪影以产生第二输出图像。例如，在一实施例中，该处理器可计算该伪影对图像强度值的贡献并从该图像强度值减去所计算的贡献。在另一实施例中，该处理器可删除代表通过将图像变换应用到第一输出图像所获得的变换数据的变换图像的预定部分，由此修正该变换数据，且可以基于经修正的变换数据来产生未变换图像。

Description

用于去除光学断面图像线条的系统和方法

相关申请的交叉引用

本申请主张2006年5月22日提出的美国专利申请11/419,566的优先权，其全部内容引用结合于此。

技术领域

本发明涉及用于从光学图像的断层去除线条的系统。

背景技术

获得三维对象的二维图像经常是期望的，例如，用于器官研究。对对象成像经常通过显微镜来进行。通过成像三维对象的特定二维平面-切片，则图像的清晰度得到增强。

传统系统通过多种不同方法产生三维对象内的二维平面的图像，包括去卷积、共焦激光扫描和光学断层(optical sectioning)。对于光学断层，传统系统将格栅图案投影在三维图像内的特定平面，并仅由该格栅图案落在的那些像素来构造图像。该平面是相对于对象来选择的一个平面。待成像对象的平面取决于对象相对于所选择的平面的位移。格栅图案是指改变光强度的图案，其可以绘制成针对像素测量的正弦波，使得峰和最低强度每隔给定数目的像素周期性出现。图1为说明用于执行光学断层的传统系统(例如显微镜)的部件的示意图。灯100发射光，光照射在水平线的格栅102上且随后作为格栅图案被分束器104反射到待成像对象上。由该对象反射的光，包括该格栅图案，随后被相机106捕获为图像。该图像由处理器108处理以产生输出图像。具体而言，处理器108提供仅由格栅图案落在的那些像素构成的输出图像。

尽管将格栅图案投影在对象上使得可以除去不属于该对象期望平面的那些像素，但是这也在得到的图像中增加了不期望的格栅图案。相应地，格栅102移动到多个位置，在这些位置的每一个得到图像，且这些图像被组合以形成无格栅线的单幅图像。提供压电驱动致动器110以移动格栅102。压电驱动致动器110响应于输入电压。该电压例如可以由处理器108产生。压电驱动致动器110移动格栅102的程度取决于施加到压电驱动致动器110的具体电压。格栅图案的特定强度投影在上面的对象的特定部分，取决于格栅102的位置。压电驱动致动器110被移动以在三个位置之间移动格栅。这些位置设置为使得所得到的相应格栅图案的强度可以绘制成相应正弦波，其中该正弦波内特定点在该三个格栅图案之间被相移相同相位角，即，0度、120度和240度的相位角，每个分隔120度。对于格栅102的三个位置的每一个，相机106捕获相应图像。图2示出相互叠加的3幅图像及其相应格栅线强度曲线。

对于每个像素，处理器108使用公式 $I_P=\mathrm{\α}\sqrt{{(I_1-I_2)}^2+{(I_2-I_3)}^2+{(I_3-I_1)}^2}$ 来组合从三幅图像的每一幅得到的值，其中I_p代表组合像素值，I₁、I₂和I₃分别代表该三幅图像的相应一幅图像的像素值，且α等于

由于格栅图案相位定为等量的120

，即，相位角为0

、120

和240

，该三幅图像内特定像素处的格栅图案的正弦波相互抵消，即，它们的值平均为零。此外，广视野图像，即，格栅图案未聚焦处的图像的部分被I₁-I₂、I₂-I₃和I₃-I₁抵消。相应地，由该三幅图像的组合确定的I_p的值不包含格栅线内相应点的值。输出图像因此不包含格栅线。

为了保证施加到压电驱动致动器110的电压将致使压电驱动致动器110把格栅102移动正确量(此时格栅图案相移120度)，某些或所有传统系统需要校准。对于校准，具有大致上均匀表面的对象，例如平滑镜子，作为待成像对象被插入，且如上所述捕获三幅图像。如果相位不正确，伪影出现在组合图像内，其中该伪影为格栅图案频率的谐波。相应地，施加到压电驱动致动器110的电压，且因此相位，被反复改变。对于每一改变，记录三幅图像，且通过快速傅立叶变换(FFT)来测量组合图像内的伪影的信号功率。这些改变被重复，直到该信号功率确定为小于特定阈值，表明伪影的基本除去，这对应于近似正确的相移。一旦得到近似正确的相移，校准完成。

该过程需要组合每组三幅图像的像素值，用于进行伪影分析。该过程通常花45秒，但是可长达5分钟。此外，相位角不是直接确定。相反，得到的是近似对应于图像位于期望相位角的情形的相位角，即，低于伪影信号阈值的减小。该过程无法精确得到期望相位角。此外，伪影信号低于阈值的情形无法使用FFT来精确地确定，特别是考虑FFT低的精确度，这至少部分可归因于以离散值测量信号功率。因此，未完全从图像除去格栅线和/或伪影。

此外，就图像强度值而言，由相机106返回的像素值经常不精确。相应地，伪影的强度的测量经常不准确。压电驱动致动器110因此被不准确地校准。

再者，与校准过程的精确度无关地，通过组合三幅图像得到的输出图像经常包含伪影。伪影经常为与格栅图案类似的图像强度的正弦变化。尽管伪影的正弦变化不一定与格栅图案频率相同，但是其通常是格栅图案的产物且位于格栅图案的正弦波的某个谐波。伪影存在诸多可能起因。示例性起因为例如压电驱动致动器110的部件的未对准，这导致图像的像素值之间的强度变化(而非格栅图案本身导致的强度变化)。这种强度变化导致三幅图案在组合时格栅图案并不抵消。其它因素也对伪影产生贡献。

此外，尽管组合三幅图像可以除去格栅线，但是该过程不产生最优图像。

相应地，本领域需要一种系统和方法，用于有效地校准格栅102的移动并提供无格栅线或伪影的最优图像。

发明内容

本发明的实施例涉及，通过确定连续地投影在待成像对象上的格栅图案的相位角，经由光学断层来产生图像的设备、计算机系统和方法。本发明的实施例涉及基于格栅图案的相位角来产生图像的设备、计算机系统和方法，其中格栅图案的相位角是参照为实际记录像素值的对数值或近似对数像素值的像素值来设置或确定的。本发明的实施例涉及基于多幅图像的值来产生图像的设备、计算机系统和方法，该图像包含超过三幅图像组合，特别是其中该多幅图像的每对连续图像中的图像是以不同相位角来得到的，即，每一幅图像的相位角均与其前一幅图像不同。此处使用的连续图像是指就格栅图案相位角而言是连续的，而非记录时间上的连续。本发明的实施例涉及从通过光学断层产生的输出图像除去伪影的设备、计算机系统和方法。

该计算机系统可包含以任何传统计算机语言写入的计算机程序。可用于实施本发明的计算机系统和方法的示例性计算机语言可以为C和/或MATLAB。

附图说明

图1为说明用于执行光学断层的传统成像系统的部件的框图。

图2说明传统系统内记录的三幅图像的叠加及其相应格栅图案强度。

图3为说明本发明示例性实施例的成像系统的示例性部件的框图。

图4为说明本发明示例性实施例的用于产生光学断面图像(optical section image)的过程的流程和数据流图。

图5说明本发明示例性实施例的用于确定输出图像像素值的，像素值的同相及正交分量与输出像素值的关系。

图6为说明本发明示例性实施例的用于产生光学断面图像的第二过程的流程和数据流图。

图7说明一示例性成像区域，该成像区域具有对同一像素行产生不同影响的格栅图案。

图8说明本发明示例性实施例的可用于确定格栅图案倾斜的频率集。

图9说明分量r、a、b和相位角的关系，其中a和b分别为幅值的余弦和正弦分量。

图10说明本发明示例性实施例的用于产生输出图像的超过三幅图像的相位角。

图11a和11b示出未变换图像和本发明实施例的变换图像的图像强度正弦变化之间的差异。

图12为说明本发明示例性实施例的通过减去代表伪影的正弦波来除去伪影的过程的流程。

图13为说明本发明示例性实施例的通过修正图像压缩数据来除去伪影的过程的流程。

图14为包括通过光学断层产生的伪影的示例性图像。

图15为依据本发明示例性实施例可产生的示例性压缩图像。

图16为本发明示例性实施例的通过从光学断面图像除去伪影可产生示例性图像。

具体实施方式

直接计算相位角

图3说明本发明实施例的成像系统的部件。图3中参考图1在上文所述的部件使用相同的参考符号。参考图3，在本发明实施例中，为了获得对象的图像，格栅102可以通过压电驱动致动器110移动到三个不同位置。将理解，可以使用压电驱动致动器以外的致动器。每个位置可处于不同相位角。对于该三个位置的每一个，相机106，例如CCD(电荷耦合装置)相机或其他传统相机，可以记录包括格栅线的相应图像。处理器108可以基于三幅记录图像产生输出图像，该处理器可以是任何恰当合适的计算机处理器或其等同物。处理器108可以是任何恰当合适的计算装置，例如计算机、个人数字助理(PDA)、膝上型计算机、笔记本计算机、移动电话、基于硬盘的装置、或者可以接收、发送和存储数据的任何装置。

可使用三个格栅位置和相应图像，以便基于与偏移120°的格栅相位角对应的图像而产生输出图像。备选地，三个格栅位置和相应图像即使未偏移120°，仍可用于为每个像素提供三个方程，每幅图像一个方程。每个方程可包括与像素值的分量相对应的三个未知变量。每个方程可以是I_n＝I_w+I_ccosφ_n+I_ssinφ_n，其中I_n代表三幅图像中的特定图像n的像素值，I_w代表像素值的广视野分量，φ_n代表特定图像n的相位角，I_c代表同相分量，以及I_s代表正交分量。如果三幅图像的相应相位角确定，则由于为仅三个未知数提供了三个方程，因此可以计算得到未知数I_w、I_c和I_s的值。

对于所记录的图像中的每一个，该系统可确定该图像的相位角，其中处理器108基于所记录的图像的组合而产生输出图像。就此而言，处理器108可将这些图像之一例如这些图像的第一幅指定为0°相位角，而不管相应格栅位置，这是因为相位角可对应于图像之间相位偏移，而不考虑格栅线相对于外部对象的移动，即，图像相位相对于彼此被测量。处理器108随后计算其余图像的各自相位角，该相位角表示与被指定为0°相位角的图像的相位之间的相位偏移。为了确定相位角，可以拍摄从大致上均匀表面反射的光的图像。例如，如果待成像的对象不具有大致上均匀表面，为了确定相位角，则需要将具有大致上均匀表面的不同对象插入相机的视线内。

在本发明实施例中，处理器108可校准致动器110以移动格栅102，以使得相位角设置为预定相位角，例如0°、120°和240°的相位角。为了校准致动器110，处理器108可致使相机106反复地记录图像集。对于该图像集中的每幅图像，处理器108可分别确定相应图像相位角并将其与预定相位角比较。基于所确定的实际相位角与预定相位角之间的偏差，处理器108可输出新电压值；依据该新电压值，电压被施加到致动器110以移动格栅102。该循环(即，施加电压到致动器110，捕获图像集，分别确定该图像集中的图像的相位角，将所确定的相位角与预定相位角比较，以及输出新电压值)，可以反复地执行，直至所确定的实际相位角在预定容差范围内与该预定相位角匹配。如果存在匹配，处理器108可终止该校准而不改变电压值。该校准可以快速地执行，因为对于每个循环，由相机106记录的图像的相位角可以直接确定。

在校准之后，处理器108例如响应于使用者指令，通过致使相机106记录三幅图像并根据公式 $I_P=\mathrm{\α}\sqrt{{(I_1-I_2)}^2+{(I_2-I_3)}^2+{(I_3-I_1)}^2}$ 来设置输出图像的每个像素的值，可产生对象的输出图像。

图4为说明本发明实施例的用于获得图像的过程的流程，该过程例如可以由如图4标题所示的处理器108、相机106和致动器110来执行。例如，通过执行作为计算机可执行指令而存储于处理器108内的计算机可读取介质、存储器312、或者本领域普通技术人员已知的成像系统的任何其他部分内的算法、软件或者方程，可以执行图4、6、12和13的过程。该计算机可读取介质可包括软盘、光盘、数字视频盘、计算机盘只读存储器(CD-ROM)等。在400，校准过程开始。在402，处理器108指示相机106记录例如三幅图像的图像集。在404，相机开始图像集的记录。在图像集的图像记录之间，处理器108在406可致使施加电压到压电驱动致动器110。响应于该电压，致动器110在408可移动格栅102。在图像记录之后，相机106在410可将记录的图像传输到处理器108。将理解，相机106可在每幅图像的记录之后传输每幅图像，或者通过单一批次传递来传输这些图像。在414，处理器108可分别确定每幅图像的图像相位角。如果处理器在416确定相位角未偏移120°，则处理器108可继续该校准过程。否则，处理器108在418结束该校准过程。

在校准之后，例如响应于使用者指令，处理器108在420开始针对输出图像的图像产生过程。对于该图像产生过程，最初可执行402-410。如果待成像对象提供足够数据以确定图像相位角，则402-410的再次执行可以省略。就此而言，如果待成像对象本身具有均匀表面，例如镜子，则可使用待成像对象来执行校准。相应地，处理器108可使用在校准过程中使用的图像数据用于该图像产生过程。此外，即使待成像对象具有非均匀表面，也有可能从对象的图像得到的数据对于校准过程是足够的。通过计算每幅图像的角频率(在下文详述)和相位角，可以比较计算结果。如果结果大致上匹配，则假定该对象已经提供足够的数据，即，对具有特定性能的校准片(calibration slide)的成像可以省略。由于待成像对象经常未提供足够数据用于确定相位角，指定对象的分别记录可被执行用于相位角确定。随后，在422，处理器108可将公式 $I_P=\mathrm{\α}\sqrt{{(I_1-I_2)}^2+{(I_2-I_3)}^2+{(I_3-I_1)}^2}$ 应用到每个像素以产生输出图像，处理器108可在424输出该输出图像。该图像可以通过任何传统输出装置例如计算机屏幕、投影仪和/或打印机来输出。

在本发明备选实施例中，校准可省略。根据该实施例，处理器108可致使相机记录具有大致上均匀表面的对象的单个图像集以确定由格栅102移动导致的图像的相位角。处理器108可将所确定的相位角保存在存储器312内。备选地，如果待成像对象具有均匀表面或者包括大量细节使得可以从对象图像获得大量数据，则处理器108可根据待成像对象的图像来确定图像相位角，而不预先成像另一对象，其中该另一对象插入到相机的视线内仅用于确定相位角。

在将所确定的相位角保存在存储器312内之后，例如响应于使用者指令，通过致使相机106记录三幅图像并将输出图像的每个像素的值设置为将所保存的相位角代入方程矩阵并求解像素值的I_c和I_s分量而得到的值，处理器108产生对象的输出图像。如上所述，对于该三幅图像的每一幅，特定像素值为I_n＝I_w+I_ccosφ_n+I_ssinφ_n。相应地，特定像素可定义为：

$\begin{array}{c}{I}_1=I_w+I_c\cos {\mathrm{\φ}}_1+I_s\sin {\mathrm{\φ}}_1\\ I_2=I_w+I_c\cos {\mathrm{\φ}}_2+I_s\sin {\mathrm{\φ}}_2\\ I_3=I_w+I_c\cos {\mathrm{\φ}}_3+I_s\sin {\mathrm{\φ}}_3\end{array}\mathrm{or}\left[\begin{array}{c}{I}_1\\ I_2\\ I_3\end{array}\right]=\left[\begin{array}{ccc}1& \cos {\mathrm{\φ}}_1& \sin {\mathrm{\φ}}_1\\ 1& \cos {\mathrm{\φ}}_2& \sin {\mathrm{\φ}}_2\\ 1& \cos {\mathrm{\φ}}_3& \sin {\mathrm{\φ}}_3\end{array}\right]\left[\begin{array}{c}{I}_w\\ I_c\\ I_s\end{array}\right]$

该方程矩阵可以如下重写以求解变量I_w、I_c和I_s：

$\left[\begin{array}{c}{I}_w\\ I_c\\ I_s\end{array}\right]={\left[\begin{array}{ccc}1& \cos {\mathrm{\φ}}_1& \sin {\mathrm{\φ}}_1\\ 1& \cos {\mathrm{\φ}}_2& \sin {\mathrm{\φ}}_2\\ 1& \cos {\mathrm{\φ}}_3& \sin {\mathrm{\φ}}_3\end{array}\right]}^{-1}\left[\begin{array}{c}{I}_1\\ I_2\\ I_3\end{array}\right]$

一旦计算得到I_c和I_s，处理器108可确定输出图像的像素值I_p，因为I_c和I_s为像素值I_pi的同相和正交聚焦(in focus)分量，如图5所示(I_W为广视野图像)。处理器108可根据公式 $I_p=\sqrt{I_c^2+I_s^2}$ (勾股定理)来确定像素值I_p。尽管像素I₁、I₂和I₃的值部分地基于投影在对象上的格栅线，I_p的值(基于分量I_c和I_s确定)完全基于对象而不是基于投影在对象上的格栅线。此外，由于精确或者大致上精确地确定相位角，根据前述方程确定的像素值I_p的组合所产生的图像不包含伪影。

图6为说明本发明该实施例的用于获得图像的过程的流程，该过程例如可以由图4标题所示的处理器108、相机106和致动器110来执行。图6中的参考图4在上文所述的部件使用相同的参考符号。根据本实施例，不执行校准。相反，仅执行相位角确定过程。400和418因此用600和618替代，且416的确定不执行。对于420，如上所述，如果待成像图像提供足够数据以确定图像相位角，402-410的再次执行可以省略。此外，由于根据本实施例并不执行用于获得偏移120°的相位角的校准，422被522替代，其中在522应用公式 $I_p=\sqrt{I_c^2+I_s^2}$ 以产生输出图像。

将理解，即使根据不执行校准过程的本实施例，处理器108可使用公式 $I_p=\sqrt{I_c^2+I_s^2}$ 计算输出图像像素。将理解，即使根据第二实施例，如果处理器108在414确定图像相位角为0°、120°和240°，处理器108也可使用公式 $I_P=\mathrm{\α}\sqrt{{(I_1-I_2)}^2+{(I_2-I_3)}^2+{(I_3-I_1)}^2}$ 来计算输出图像像素。

相应地，通过确定三幅图像的相位角，可以快速地执行校准。此外，通过确定相位角，可以基于处于不同相位角的图像的集合来产生输出图像，即使不校准致动器110以致使该图像集的图像的格栅线位于预定相位角。

参考图4和6，在本发明实施例中，处理器108可在412确定该图像集的图像的格栅线的角频率，并可以基于图像的像素值与所确定的频率的关联来计算该图像集的图像的相位角，如下所述。具体参考图4，尽管在本发明一个实施例中，可通过比较所确定的频率在用于质量控制的校准过程的每次迭代期间执行412，但是在备选实施例中，在校准过程的除了第一次迭代的每次迭代期间可以省略412，因为一旦频率已知，就无需重新计算。将理解，频率不被固定。例如，频率可依赖于反射图像或者反射到对象上的光的放大倍数，该放大倍数可取决于透镜的位置。为了通过像素值与所确定频率的关联来计算相位角，需要高度精确地确定频率。例如，使用FFT对于频率的确定可能是不够的。在本发明示例性实施例中，处理器108可使用贝叶斯频谱分析(Bayesian Spectral Analysis)高精度地估计频率，本领域技术人员公认的是，贝叶斯频谱分析这种分析提供了比使用FFT获得的离散值结果更为流畅的结果。

为了应用贝叶斯频谱分析，可以收集图像的信号数据。每个信号可以使用与图像强度的正弦变化相关的方程来表述。该方程可以是f(y_i)＝rcos(ωy_i+φ)+c，其中r为幅值，ω为所确定的角频率，y为像素位置，φ为相位角，以及c为图像强度的平均值。关于y，将理解，其可以是沿垂直方向或者沿水平方向的像素坐标，取决于格栅线的取向。例如，格栅102的取向可以使得格栅线水平地投影在图像上，由此导致沿垂直方向的图像强度变化。在这种情形下，像素坐标为沿垂直方向的像素坐标。图像强度的正弦变化也可以用f(y_i)＝acosωy_i+bsinωy_i+c表示，其中a和b为幅值的余弦和正弦分量。尽管前述两个方程是等价的，前一个方程仅包含两个未知数ω和φ，而后一个方程仅包含这两个未知数之一，即ω。相应地，使用后一个公式，例如通过下述方式通过贝叶斯频谱分析可以确定角频率ω。

将后一个公式应用到多个数据样本‘d’，可以得到下述矩阵表述：

$\left[\begin{array}{c}{d}_1\\ d_2\\ d_3\\ \·\\ \·\\ \·\\ d_N\end{array}\right]=\left[\begin{array}{ccc}{\cos \mathrm{\ωy}}_1& \sin \mathrm{\ω}{y}_1& 1\\ {\cos \mathrm{\ωy}}_2& \sin \mathrm{\ω}{y}_2& 1\\ {\cos \mathrm{\ωy}}_3& \sin \mathrm{\ω}{y}_3& 1\\ \·& & \\ \·& & \\ \·& & \\ {\cos \mathrm{\ωy}}_N& \sin \mathrm{\ω}{y}_N& 1\end{array}\right]\left[\begin{array}{c}a\\ b\\ c\end{array}\right]+\left[\begin{array}{c}{e}_1\\ e_2\\ e_3\\ \·\\ \·\\ \·\\ e_N\end{array}\right]$

因此可以得到一矩阵，其中：

$G=\left[\begin{array}{ccc}{\cos \mathrm{\ωy}}_1& \sin \mathrm{\ω}{y}_1& 1\\ {\cos \mathrm{\ωy}}_2& \sin \mathrm{\ω}{y}_2& 1\\ {\cos \mathrm{\ωy}}_3& \sin \mathrm{\ω}{y}_3& 1\\ \·& \·& \\ \·& \·& \\ \·& \·& \\ {\cos \mathrm{\ωy}}_N& \sin \mathrm{\ω}{y}_N& 1\end{array}\right]$

可以积分(integrate out)得到线性系数和噪声标准偏差。然后通过将G矩阵应用到贝叶斯公式 $p\left(\mathrm{\ω}\right|d,I)\∝\frac{{[d^{T}d-d^{T}G{\left(G^{T}G\right)}^{-1}{G}^{T}d]}^{\frac{M-N}{2}}}{\sqrt{\det \left(G^{T}G\right)}},$ 用于针对数据集d确定角频率ω的可能值(probable value)，可得到频率。M为G矩阵中包含的列的数目。这些图像中单个图像的样本对于确定频率是足够的。然而，通过输入超过一幅图像的数据，可以获得更高的精确度。

在本发明一个示例性实施例中，图像的窄条可以用作数据源，而不是整幅图像。例如，如果格栅线水平地投影在图像上，导致沿垂直方向的图像强度变化，图像的垂直条可用作数据源。尽管与使用整幅图像相比，使用窄条提供更少的数据用于输入到方程内，但是这可以增加精确度，因为格栅线可以相对于成像区域成角度地投影，由此使数据输入歪斜，如下文所详述。

在本发明备选实施例中，对于用作输入到贝叶斯频谱分析的数据，可以使用所有的图像像素的值。

在又一实施例中，可以对每一像素行(或者像素列，如果格栅线垂直投影)求和。每行的像素值和

可以用作贝叶斯频谱分析的数据输入，尽管这些和为一个像素宽的窄垂直条的实际像素值，即，每行一个值。与前一个实施例相比，后两个实施例在频率测量方面可以提供更高的精确度，因为某些记录像素值例如由于噪声而不正确。因此，考虑越多的像素值，频率估计将更精确。

然而，关于后两个实施例，如果格栅102置为使得格栅线相对于成像区域成角度投影，则由投影格栅图案导致的图像强度的正弦变化将不等地影响沿同一行的像素，如上所述。图7说明这种现象。图7示出格栅图案702成角度地投影在其上的成像区域700。对于选定行704，格栅图案702不同地影响行704的不同像素的值。例如，尽管像素2示为未受影响，但是像素13显示为受很大影响。相应地，如果选定行704的像素值被求和或者以其他方式组合用作频率确定的数据输入，所得到的频率值可能不正确。

因此，在本发明示例性实施例中，在确定格栅图案的频率之前，图像可被旋转使得格栅图案平行于成像区域的横坐标(或者平行于纵坐标，如果格栅线垂直投影)。格栅线相对于横坐标投影的角度(α)可确定，且图像可旋转-α。

在本发明实施例中，格栅线相对于成像区域的横坐标投影的弧度角可以如下确定：(a)确定格栅图案沿两个任意对角线的频率，所述对角线叠加在图像上且相对于成像区域的横坐标以相反的角度45°和-45°布置，以及(b)将所确定的频率应用到公式

其中F₁为沿对角线之一所确定的频率，F₂为沿另一对角线所确定的频率。当对角线布置在45°和-45°时，通过应用前述公式可获得倾斜角α，因为tan(45°)＝1且因为

因此，对角线角为45°时， $\tan \left(\mathrm{\α}\right)=\frac{F_1-F_2}{F_1+F_2}.$ 为了确定沿对角线的频率，沿选定对角线的一条像素的值可以输入到上述的频率确定公式，例如，它们可以用作贝叶斯频谱分析的输入。

例如，图8示出相对于成像区域700的横坐标在45°角绘制的对角线801和在-45°角绘制的对角线802。当格栅图案成角度布置时，格栅线的实际频率可能不同于沿对角线得到的频率。此外，沿两个对角线得到的频率会不同。例如，沿实际频率线800得到的格栅线频率代表图8中的实际格栅频率，约为图8的测量单位的0.5603。然而，沿对角线801和802测量的格栅线频率分别约为测量单位的1.1420和0.6402。

将这些频率应用到

可以确定相对于图8中横坐标804的格栅倾斜角弧度约为

注意，图8示出对角线801和802以及实际频率线800的格栅线周期。然而，通过2π/周期将周期转换成角频率得到相同结果。例如，对于周期1.1420，角频率为5.50191358，以及对于周期0.6402，角频率为9.81441004。将这些频率的周期代入方程

类似地得到15.73°的倾斜角。使用任何合适恰当的过程，图像可倾斜15.73°，使得格栅线可以平行于横坐标804。

一旦倾斜得到校正，该系统和方法可以按照上述方式确定格栅线频率。

一旦频率确定，则图像的相位角可以确定。对于图像的像素值，acosωy_i+bsinωy_i+c的a和b分量可以通过使用像素值到所确定频率的线性回归来估计。一旦a和b被估计，图像的相位角可以根据图9所示的关系而计算为arctan(b/a)。可以进行任意单个图像的相位角的确定而不使用与图像集中的其他图像有关的数据。例如，参考图4和6，只要从相机106接收到图像即可执行412和414，即使相机106在每幅图像的记录之后分别立即传输该图像。相应地，在致动器110移动格栅102以准备记录后一图像和/或当相机106记录后一图像时，处理器108可对先前接收的图像执行412和414。

使用超过三幅图像

如上所详述，通过基于图像集的相应像素值的组合确定像素值，可以执行图像产生过程，其中对于每一图像，格栅线以不同相位角被投影。尽管用于产生输出图像的图像集通常包含三幅图像，在本发明实施例中，为了获得更好质量图像，处理器108可以基于超过三幅图像的像素值来产生输出图像。例如，相位角之间的偏移可以减小，如图10所示。图10示出图像之间偏移30°相位角。为了清楚，仅示出单幅图像即参考图像的强度曲线。虚线表示其余图像强度曲线的开始。根据本实施例，矩阵公式 $\left[\begin{array}{c}{I}_1\\ I_2\\ I_3\end{array}\right]=\left[\begin{array}{ccc}1& \cos {\mathrm{\φ}}_1& \sin {\mathrm{\φ}}_1\\ 1& \cos {\mathrm{\φ}}_2& \sin {\mathrm{\φ}}_2\\ 1& \cos {\mathrm{\φ}}_3& \sin {\mathrm{\φ}}_3\end{array}\right]\left[\begin{array}{c}{I}_w\\ I_c\\ I_s\end{array}\right]$ 可以用 $\left[\begin{array}{c}{I}_1\\ \·\\ \·\\ \·\\ I_M\end{array}\right]=\left[\begin{array}{ccc}1& \cos {\mathrm{\φ}}_1& \sin {\mathrm{\φ}}_1\\ \·& & \\ \·& & \\ \·& & \\ 1& \cos {\mathrm{\φ}}_M& \sin {\mathrm{\φ}}_M\end{array}\right]\left[\begin{array}{c}{I}_w\\ I_c\\ I_s\end{array}\right]$ 来替代。P16

如上所述确定相位角，超过三幅图像的图像集提供了比未知数多的方程，因为仅I_w、I_c和I_s是未知数。由于噪声，这些方程可能不完全一致。相应地，例如最小二乘方回归的回归分析可应用于I_w、I_c和I_s，这可以减小信号中存在的噪声。具体而言，可以应用下述最小二乘方回归公式：

$\left[\begin{array}{c}{I}_w\\ I_c\\ I_s\end{array}\right]={\left(G^{T}G\right)}^{-1}{G}^T\left[\begin{array}{c}{I}_1\\ I_2\\ I_3\\ \·\\ \·\\ \·\\ I_M\end{array}\right],$

其中

${\left(G^{T}G\right)}^{-1}=\left[\begin{array}{ccc}\frac{1}{M}& 0& 0\\ 0& \frac{2}{M}& 0\\ 0& 0& \frac{2}{M}\end{array}\right]$

$G=\left[\begin{array}{ccc}1& \cos {\mathrm{\φ}}_1& \sin {\mathrm{\φ}}_1\\ \·& & \\ \·& & \\ \·& & \\ 1& \cos {\mathrm{\φ}}_M& \sin {\mathrm{\φ}}_M\end{array}\right]$

且G^T为G的转置矩阵。即使在仅使用三幅图像时，仍可应用该公式。

如果超过三幅图像中的每对连续图像的相位角偏移相同角度数，则可以应用其他公式。与图像集的图像数目(M)无关，I_w、I_c和I_s可以计算如下：

$I_w=\frac{1}{M}{\mathrm{\Σ}}_{k=0}^{M-1}{I}_{k+1}$

$I_c=\frac{2}{M}{\mathrm{\Σ}}_{k=0}^{M-1}{I}_{k+1}\cos \left(\frac{2\mathrm{\πk}}{M}\right)$

$I_s=\frac{2}{M}{\mathrm{\Σ}}_{k=0}^{M-1}{I}_{k+1}\sin \left(\frac{2\mathrm{\πk}}{M}\right)$

即使当M＝3时，也可以应用该公式。一旦I_c和I_s使用前述两个公式中任一个来计算，I_p可以使用公式 $I_p=\sqrt{I_c^2+I_s^2}$ 来计算。此外，如果使用四幅图像且四幅图像中的每对连续图像的相位角偏移相同角度数，则I_p可以使用公式 $I_P=\sqrt{{(I_1-I_3)}^2+{(I_2-I_4)}^2}$ 来计算。

在本发明实施例中，根据用于更新最小二乘方解的传统过程来修正最小二乘方解，所产生的图像的像素值可被递归更新以说明新得到的图像。例如，递归最小二乘方公式可包括传统上使用的公式，例如Sherman-Morrison公式或者Woodbury公式。相应地，在基于三幅或更多幅图像的像素数据的图像被输出之后，使用者可以指示处理器108产生更增强的图像。相应地，处理器108可获得新记录的图像(包括格栅图案)，且可更新已经计算的I_c和I_s的值而不使用先前使用的图像来重新进行该计算。相应地，对于期望更新的情形，不需要存储先前使用的图像。

预处理—像素数据的转换用于估计参数

由相机106返回的图像的像素值经常提供图像强度的不均匀正弦变化。相应地，如果是基于由相机106记录的像素值，则校准致动器110以提供特定相位角，无论是基于用FFT对伪影的测量还是基于相位角的直接计算，和/或计算相位角用以基于 $I_p=\sqrt{I_c^2+I_s^2}$ 来产生输出图像，均可能有缺陷。在本发明实施例中，该系统可以使用通过像素值的对数或近似对数转换得到的值来替代用于校准或用于确定相位角(和/或频率)的每个记录像素值。得到的值可提供图像强度的更均匀正弦变化。图11示出未变换图像和变换图像的图像强度的正弦变化之间的差异。曲线(a)代表未变换图像的图像强度，曲线(b)代表变换图像的图像强度，曲线(a)的正弦波的幅值显著不均匀，但是曲线(b)的正弦波的幅值明显更为均匀。

在转换之后，执行传统的校准或者根据直接计算的相位角的校准。备选地，可以计算相位角而不进行前文所详述的校准。在校准和/或相位角的计算之后，处理器108可根据前文所详述的过程基于未变换的即原始记录的像素值来产生输出图像。

在本发明一个实施例中，为了转换所记录的像素值，可以将每个像素简单变换为其对数值。根据本实施例，可以实现相反的效果，其中在低图像强度下的噪声被放大，这使图像强度值失真。在备选实施例中，可以对每个像素使用逆双曲线正弦函数sinh^-1(x/2)，其中x为原始记录的图像强度值。对大的像素值，后一函数近似以e为底的函数log(x)(自然对数)，但对于较小的像素值，并不如此。根据本实施例，低图像强度下噪声的放大可以避免。将理解，可以使用跨过图像使强度的正弦变化的幅值变得平滑的任何函数来进行像素值的转换。

预处理—从输出图像除去伪影

不管校准过程或者相位角确定的精确度如何，通过组合三幅或者更多幅图像而得到的输出图像可包括伪影。伪影可以是与格栅图案相似的正弦变化。伪影的正弦变化可以是格栅图案的产物，且可以位于格栅图案正弦波的特定谐波。尤其是，可以假定伪影的正弦变化在格栅图案频率的三个谐波之内。

在本发明的示例性实施例中，该系统和方法可以从输出图像除去由伪影引起的图像强度的正弦变化。

在本发明的一个示例性实施例中，可以确定代表由于伪影引起的图像强度的正弦变化的正弦波。从该图像减去所确定的正弦波，得到无伪影的图像。输出图像的像素值可用Q＝I+B表示，其中Q为像素值，I为由被成像对象贡献的像素值的部分，以及B为由伪影贡献的像素值的部分。对于每个像素值，B的相应值可被确定并从Q减去以得到I，即，无任何伪影贡献的像素值。

如上文所详述，由于投影格栅图案引起的图像强度的正弦变化可以用f(y_i)＝acosωy_i+bsinωy_i+c表示。类似地，B(在y行伪影对像素值的贡献，假定垂直投影的格栅图案和伪影)可以用a₁cos(ωy_i)+b₁sin(ωy_i)+a₂cos(2ωy_i)+b₂sin(2ωy_i)+a₃cos(3ωy_i)+b₃sin(3ωy_i)表示。每个cos/sin组可对应于伪影。可以假设伪影为谐波1至3中的一个或多个。因此，该系统和方法可以假定包含三个cos/sin组的上述方程代表图像内的伪影，如果有任何伪影。

因此，在特定行y_i的像素值可以用Q(x，y_i)＝I(x，y_i)+a₁cos(ωy_i)+b₁sin(ωy_i)+a₂cos(2ωy_i)+b₂sin(2ωy_i)+a₃cos(3ωy_i)+b₃sin(3ωy_i)表示。例如，取沿垂直条的像素，像素的值可以表示如下：

Q₁＝I₁+a₁cos(ω1)+b₁sin(ω1)+a₂cos(2ω1)+b₂sin(2ω1)+a₃cos(3ω1)+b₃sin(3ω1)；

Q₂＝I₂+a₁cos(ω2)+b₁sin(ω2)+a₂cos(2ω2)+b₂sin(2ω2)+a₃cos(3ω2)+b₃sin(3ω2)；

Q₃＝I₃+a₁cos(ω3)+b₁sin(ω3)+a₂cos(2ω3)+b₂sin(2ω3)+a₃cos(3ω3)+b₃sin(3ω3)；

.

.

.

Q_n＝I_n+a₁cos(ωn)+b₁sin(ωn)+a₂cos(2ωn)+b₂sin(2ωn)+a₃cos(3ωn)+b₃sin(3ωn)

这些方程可以重新写成q＝b*G的矩阵形式，其中：

$\left[\begin{array}{c}{q}_1\\ q_2\\ q_3\\ \·\\ \·\\ \·\\ q_n\end{array}\right],b=\left[\begin{array}{c}I\\ a_1\\ b_1\\ a_2\\ b_2\\ a_3\\ b_3\end{array}\right],$

以及

$G=\left[\begin{array}{ccccccc}1& \cos \left({\mathrm{\ωy}}_1\right)& \sin \left({\mathrm{\ωy}}_1\right)& \cos \left({2\mathrm{\ωy}}_1\right)& \sin \left({2\mathrm{\ωy}}_1\right)& \cos \left({3\mathrm{\ωy}}_1\right)& \sin \left(3{\mathrm{\ωy}}_1\right)\\ 1& \cos \left({\mathrm{\ωy}}_2\right)& \sin \left({\mathrm{\ωy}}_2\right)& \cos \left({2\mathrm{\ωy}}_2\right)& \sin \left({2\mathrm{\ωy}}_2\right)& \cos \left({3\mathrm{\ωy}}_2\right)& \sin \left(3{\mathrm{\ωy}}_2\right)\\ 1& \cos \left({\mathrm{\ωy}}_3\right)& \sin \left({\mathrm{\ωy}}_3\right)& \cos \left({2\mathrm{\ωy}}_3\right)& \sin \left({2\mathrm{\ωy}}_3\right)& \cos \left({3\mathrm{\ωy}}_3\right)& \sin \left(3{\mathrm{\ωy}}_3\right)\\ \·& & & & & & \\ \·& & & & & & \\ \·& & & & & & \\ 1& \cos \left({\mathrm{\ωy}}_n\right)& \sin \left({\mathrm{\ωy}}_n\right)& \cos \left({2\mathrm{\ωy}}_n\right)& \sin \left({2\mathrm{\ωy}}_n\right)& \cos \left({3\mathrm{\ωy}}_n\right)& \sin \left(3{\mathrm{\ωy}}_n\right)\end{array}\right]$

尽管I的值取决于确切的像素位置(x，y)，然而对于确定由伪影引起的强度的正弦变化，这不重要。因此，对于上述矩阵，I乘以1，而与行无关。

如果格栅图案的角频率确定，则该G矩阵的每个分量的值可以已知。相应地，该系统和方法可首先如上所述确定角频率。例如，像素值可以输入贝叶斯频谱分析。这包括旋转图像以产生与横坐标平行的格栅图案以及对每一行内的像素值求和，如参考图8在上文所述。在角频率被确定以用于确定相位角的实施例中，例如在校准期间，先前确定的频率可以代入方程，而不再次执行角频率确定。在后处理(postprocess)图像之前不确定角频率的实施例中，该系统和方法可以执行角频率确定用于该后处理过程。根据该实施例，该系统和方法可以取回输入图像的存储像素值以确定格栅图案的角频率，并将所确定的频率输入到G矩阵，但是可以取回输出图像的像素值用于输入到q矩阵。

因此，通过将上述矩阵输入回归分析，可以确定(b矩阵的)a1、b1、a2、b2、a3和b3的值。例如，可以使用该G矩阵应用下述最小二乘方回归公式：

$\left[\begin{array}{c}I\\ a_1\\ b_1\\ a_2\\ b_2\\ a_3\\ b_3\end{array}\right]={\left(G^{T}G\right)}^{-1}{G}^T\left[\begin{array}{c}{q}_1\\ q_2\\ q_3\\ \·\\ \·\\ \·\\ q_n\end{array}\right],$

其中G^T为G的转置矩阵。尽管在上文示出了每行(或者列，如果格栅图案是由垂直线形成)一个像素样本，每行(或列)的额外像素样本可以输入到q矩阵内，且G矩阵可相应地包含额外的行。因此G矩阵的多个行可涉及相同的ωy_i值。

一旦矩阵b的值确定，则下述的值：

a₁cos(ω1)+b₁sin(ω1)+a₂cos(2ω1)+b₂sin(2ω1)+a₃cos(3ω1)+b₃sin(3ω1)；

a₁cos(ω2)+b₁sin(ω2)+a₂cos(2ω2)+b₂sin(2ω2)+a₃cos(3ω2)+b₃sin(3ω2)；

a₁cos(ω3)+b₁sin(ω3)+a₂cos(2ω3)+b₂sin(2ω3)+a₃cos(3ω3)+b₃sin(3ω3)；

.

.

.

a₁cos(ωn)+b₁sin(ωn)+a₂cos(2ωn)+b₂sin(2ωn)+a₃cos(3ωn)+b₃sin(3ωn)

可以通过b*G确定。这些值可以是伪影对像素的贡献。伪影在不同行的贡献不同以及因此在不同行方程值不同，原因可能是由于图像强度的垂直正弦变化，其中格栅图案垂直投影。

一旦这些值确定，可以从输出图像像素值减去上述方程的值。例如，对于行1的每个像素，可以从像素值减去a₁cos(ω1)+b₁sin(ω1)+a₂cos(2ω1)+b₂sin(2ω1)+a₃cos(3ω1)+b₃sin(3ω1)的值。对于正弦波减去阶段，该系统和方法可以赋予每行与输入到回归分析时该行被赋予的相同的行号。例如，如果对于回归分析，成像区域第一行的像素被指定为用“行1”标记的行，则可以从成像区域第一行的像素减去a₁cos(ω1)+b₁sin(ω1)+a₂cos(2ω1)+b₂sin(2ω1)+a₃cos(3ω1)+b₃sin(3ω1)的值。否则，会从每一行减去错误的值。

输出图像中的伪影可能是由格栅图案的少于三个谐波构成，这种情况下，a₁cos(ωi)+b₁sin(ωi)+a₂cos(2ωi)+b₂sin(2ωi)+a₃cos(3ωi)+b₃sin(3ωi)的一些部分可能等于0。例如，如果伪影仅为格栅图案的一次谐波，则a₂cos(2ωi)+b₂sin(2ωi)+a₃cos(3ωi)+b₃sin(3ωi)等于0。

在本发明的一个实施例中，总是假定3个谐波。这种情况下，回归分析会得到接近或等于0的a₂、b₂、a₃和b₃的值。

在本发明的另一实施例中，不采用3个谐波假定，可根据用于确定谐波数目的传统过程，最初确定形成伪影的格栅图案频率的谐波数目。该系统和方法可以根据这种确定来改变矩阵结构。例如，如果确定伪影包括格栅图案频率的两个谐波的分量，则不采用上述矩阵结构，该系统和方法可将n×5G矩阵而非n×7矩阵输入该回归分析。该矩阵可具有下述结构：

$G=\left[\begin{array}{ccccc}1& \cos \left({\mathrm{\ωy}}_1\right)& \sin \left({\mathrm{\ωy}}_1\right)& \cos \left({2\mathrm{\ωy}}_1\right)& \sin \left({2\mathrm{\ωy}}_1\right)\\ 1& \cos \left({\mathrm{\ωy}}_2\right)& \sin \left({\mathrm{\ωy}}_2\right)& \cos \left({2\mathrm{\ωy}}_2\right)& \sin \left({2\mathrm{\ωy}}_2\right)\\ 1& \cos \left({\mathrm{\ωy}}_3\right)& \sin \left({\mathrm{\ωy}}_3\right)& \cos \left({2\mathrm{\ωy}}_3\right)& \sin \left({2\mathrm{\ωy}}_3\right)\\ \·& & & & \\ \·& & & & \\ \·& & & & \\ 1& \cos \left({\mathrm{\ωy}}_n\right)& \sin \left({\mathrm{\ωy}}_n\right)& \cos \left({2\mathrm{\ωy}}_n\right)& \sin \left({2\mathrm{\ωy}}_n\right)\end{array}\right]$

该矩阵使得可以确定系数a1、a2、b1和b2，而不确定系数a3和b3，因为基于谐波数目确定可以假定后一对系数等于0。这可以提高确定前两对系数的值的精确度，因为回归分析不将任何值归因于虚构的第三对系数。

在本发明的示例性实施例中，输出图像的窄条例如垂直条可用作q矩阵的输入。备选地，可以使用图像的更多部分。在一个特定实施例中，例如每行的像素值之和可用作q矩阵的输入。在本发明的一个特定实施例中，输出图像可旋转以使得伪影，如果有任何伪影，平行于成像区域的横坐标，如上文参考图8有关角频率确定和格栅图案所描述。

在本发明示例性实施例中，不使用实际输出图像像素值用于q矩阵，该系统和方法可使用预处理输出图像像素值用于q矩阵。这些输出图像像素值可通过对像素值的对数或近似对数转换来预处理，以得到图像强度的更为均匀的正弦变化，如上文参考图11所述。类似地，对于G矩阵的格栅图案角频率确定，可使用输入图像的预处理像素值，而不使用输入图像的实际像素值。

根据本实施例，如果有任何伪影，则可以从经预处理的输出图像像素值减去伪影对输出图像的贡献，而不是从实际输出图像像素值减去。在该减去之后，该系统和方法可将修正的预处理像素值输入到用于逆对数转换的方程，以得到修正的实际输出图像像素值，即，如果有任何伪影，通过除去伪影修正。

图12为说明本发明示例性实施例的用于从输出图像除去伪影的过程的流程，例如可以由处理器108来执行该过程。在1202，包括输入图像1200a和输出图像1200b的图像1200可被预处理以输出经预处理图像1203，该经预处理图像1203包括图像强度的正弦变化更为均匀的经预处理输入图像1203a和经预处理输出图像1203b。在1204，经预处理图像1203可被旋转，使得输入图像1203a的投影格栅线图案和输出图像1203b的伪影平行于其相应成像区域的横坐标(假定水平格栅线的垂直格栅图案)，如果他们之前不是平行的。在1206，包括旋转输入图像1205a和旋转输出图像1205b的旋转图像1205的一些或所有行的像素值(假定水平格栅线的垂直格栅图案)，可被求和以产生像素列1207，该像素列1207包括输入图像像素列1207a和输出图像像素列1207b。

在1206之后，1208和1210可顺序或同时执行，因为二者的执行相互独立。在1208，基于输入图像像素列1207a来确定输入图像的格栅图案的角频率1209。在1210，基于输出图像像素列1207b来确定形成伪影的谐波数目1211。

在1212，基于输出图像像素列1207b(对于q矩阵)、频率1209(对于G矩阵)和谐波数目1211(对于G矩阵)，确定经预处理输出图像像素值的伪影贡献1213。在1214，可以从经预处理输出图像1203b的像素值减去伪影贡献1213以得到修正的经预处理输出图像1215。在1216，对修正的经预处理输出图像1215执行图像预处理的逆处理以产生修正的输出图像，其中该修正的输出图像实际而言与不包含伪影的输出图像1200b相同，如果有任何伪影。

在本发明实施例中，其中通过计算伪影贡献并从输出图像的像素值减去该伪影贡献以除去输出图像的伪影，伪影除去可以通过1208-1214来执行，而不执行1202-1206、1210和1216中的一个或多个，例如全部。例如，未旋转输出图像1200b的实际像素值，而不是经预处理像素值和/或行的和，可用于q矩阵。类似地，可以基于未旋转、未预处理以及未行求和的原始输入图像1200a来确定输入图像1200a的格栅图案的角频率。类似地，可以使用未预处理输出图像1200b、经预处理输出图像1203b或者未预处理输出图像像素列，而不是输出图像像素列1207b，作为用于确定谐波数目1211的输入。此外，不是确定谐波数目1211，而假定一数目，例如3。最后，如果输出图像1200b未被预处理，则预处理的逆处理可略去。例如，1214可产生修正的输出图像1217。

此外，在本发明的示例性实施例中，不在1208确定格栅图案频率用于伪影除去，而可以在1208获得在伪影除去过程之前确定的格栅图案频率。

在本发明备选示例性实施例中，为了从光学断层输出图像除去伪影，该系统和方法可除去该输出图像的代表图像变换数据的图像(这里称为变换图像)的断面(section)，其中该断面位于该变换图像的预定位置，即，位于该预定位置的形成变换图像的该部分的图像变换数据的部分可被除去。

图13为说明本发明示例性实施例的通过相应变换图像以从输出图像除去伪影的过程的流程，例如可以由处理器108来执行该过程。在1300，执行图像变换以形成图像变换数据。可以使用的示例性图像变换的非穷举清单为小波变换、离散余弦变换(DCT)、离散正弦变换(DST)、快速傅立叶变换(FFT)、阿达玛(Hadamard)变换、哈特利(Hartley)变换和离散小波变换(DWT)。示例性的DWT变换可以是Haar小波和Daubechies小波。

在1302，可以产生代表图像变换数据的变换图像。例如，图14为包含多条水平线的伪影1402的示例性输出图像1400，某些水平线被挑选出来。图15为示例性变换图像1500，该变换图像1500代表通过将图像变换应用到输出图像1400而产生的变换数据。

在1304，该系统和方法可以除去所产生的变换图像的预定断面。对于代表输出图像的变换数据的任何变换图像，该预定断面可以是同一限定区域。该预定断面可以是变换图像的一部分，变换图像的该部分对应于被确定为具有低水平频率和高垂直频率(其中格栅线且因此伪影被水平地投影)的输出图像的元素，这是与输出图像的无伪影元素相比而言的。可以假定无伪影元素将不会表现在变换图像的该断面内。具体而言，如果变换图像是下述情形，即，变换图像的像素的行号越小，该像素所对应的图像元素沿垂直方向的频率越低，并且变换图像的像素的列号越小，该像素所对应的图像元素沿水平方向的频率越低，则该断面可包含在变换图像的左手侧开始的约1至2％连续像素列的底部约80％的像素。然而，可以根据伪影的预期频率来选择其他像素和列百分比。示例性伪影断层1502描绘于图15。伪影部分(artefact section)的除去可将与该伪影部分相对应的数据设置为0。

在1306，该系统和方法可执行逆变换以基于修正的变换数据来产生未变换图像。在1306产生的图像基本上与排除伪影的输出图像相同。例如，图16为由与图14输出图像1400相对应的修正的变换数据所产生的示例性修正的输出图像1600。修正的输出图像1600不包含图14所示的伪影1402。

本领域技术人员从前述说明可以理解，本发明可按照各种形式来实施。因此，尽管已经结合其具体示例来描述本发明的实施例，本发明的实施例的真正范围不应限于这些实施例，因为本领域技术人员在研究附图、说明书和所附权利要求书之后可以显见其他变形。

Claims (56)

1.一种图像产生方法，包括：

基于多幅输入图像产生第一输出图像；以及

基于该第一输出图像产生第二输出图像，该第二输出图像与减去伪影的第一输出图像相同。

2.如权利要求1所述的方法，还包括：

确定该伪影对该第一输出图像的图像强度值的贡献；

其中产生该第二输出图像包括从该图像强度值减去该贡献。

3.如权利要求2所述的方法，其中该图像强度值为像素强度值。

4.如权利要求2所述的方法，其中该贡献沿水平和垂直方向之一正弦变化。

5.如权利要求4所述的方法，其中：

确定该贡献包括确定代表该正弦变化的方程的系数的值；

该减去包括：

通过将所确定的系数代入所述方程来确定该方程的值；以及

从该图像强度值减去所述方程的值；以及

所述方程的值在像素行和像素列之一之间变化。

6.如权利要求5所述的方法，其中：

该伪影由所述输入图像中的至少一幅图像的格栅图案的多个谐波m形成；以及

该方程为a₁cos(ωy)+b₁sin(ωy)+...a_mcos(mωy)+b_msin(mωy)，其中a₁...a_m和b₁...b_m为所述系数，ω为该格栅图案的角频率，以及y为对于每一像素行是不同的像素行号和对于每一像素列是不同的像素列号之一。

7.如权利要求6所述的方法，其中：

所述系数通过将回归分析应用到矩阵Q、b和G来确定；

$Q=\left[\begin{array}{c}{q}_1\\ q_2\\ q_3\\ \·\\ \·\\ \·\\ q_x\end{array}\right],$ 其中q为该第一输出图像的采样像素强度值，x为被采样像素数目；

$b=\left[\begin{array}{c}I\\ a_1\\ b_1\\ \·\\ \·\\ \·\\ a_m\\ b_m\end{array}\right],$ ,其中I代表减去该伪影对采样像素强度的贡献的采样像素强度值；以及

$G=\left[\begin{array}{ccccccc}1& \cos \left({\mathrm{\ωy}}_i\right)& \sin \left({\mathrm{\ωy}}_i\right)& \cos \left({2\mathrm{\ωy}}_i\right)& \sin \left(2{\mathrm{\ωy}}_i\right)& \cos \left({3\mathrm{\ωy}}_i\right)& \sin \left({3\mathrm{\ωy}}_i\right)\\ 1& \cos \left({\mathrm{\ωy}}_i\right)& \sin \left({\mathrm{\ωy}}_i\right)& \cos \left({2\mathrm{\ωy}}_i\right)& \sin \left({2\mathrm{\ωy}}_i\right)& \cos \left({3\mathrm{\ωy}}_i\right)& \sin \left(3{\mathrm{\ωy}}_i\right)\\ 1& \cos \left({\mathrm{\ωy}}_i\right)& \sin \left({\mathrm{\ωy}}_i\right)& \cos \left(2{\mathrm{\ωy}}_i\right)& \sin \left({2\mathrm{\ωy}}_i\right)& \cos \left({3\mathrm{\ωy}}_i\right)& \sin \left(3{\mathrm{\ωy}}_i\right)\\ \·& & & & & & \\ \·& & & & & & \\ \·& & & & & & \\ 1& \cos \left({\mathrm{\ωy}}_i\right)& \sin \left({\mathrm{\ωy}}_i\right)& \cos \left({2\mathrm{\ωy}}_i\right)& \sin \left(2{\mathrm{\ωy}}_i\right)& \cos \left({3\mathrm{\ωy}}_i\right)& \sin \left(3{\mathrm{\ωy}}_i\right)\end{array}\right]$

其中G矩阵的每一行与q矩阵的相应值相关联，且其中，对于G矩阵的每一行，y_i为G矩阵行的相应q矩阵值的行号和列号之一。

8.如权利要求7所述的方法，其中通过应用公式b＝(G^TG)^-1G^T使用最小二乘方回归来确定所述系数，其中G^T为G的转置矩阵。

9.如权利要求7所述的方法，还包括：

变换该第一输出图像的图像像素值以跨过该第一输出图像平滑图像强度正弦变化的幅值，该正弦变化代表该伪影，其中：

变换的值被采样以用于回归分析；

从变换的值减去所述贡献；以及

产生该第二输出图像还包括，在从所述图像强度值减去该贡献

之后，执行图像像素值变换的逆变换。

10.如权利要求9所述的方法，其中该第一输出图像的图像像素值变换成其相应对数值。

11.如权利要求9所述的方法，其中：

通过应用逆双曲线正弦函数到(q′/2)，变换该第一输出图像的图像像素值；以及

q′代表未变换像素值。

12.如权利要求7所述的方法，还包括：

确定该第一输出图像内的伪影相对于该第一输出图像的成像区域的倾斜；

旋转该第一输出图像以消除该倾斜；

对于该第一输出图像的多个像素行中的每一行，对该行的像素值求和；以及

将该行像素值之和输入所述回归分析。

13.如权利要求6所述的方法，其中角频率通过贝叶斯频谱分析来确定。

14.如权利要求13所述的方法，其中：

角频率通过应用公式 $p\left(\mathrm{\ω}\right|d,I)\∝\frac{{[d^{T}d-d^{T}G\′{({G\′}^{T}G\′)}^{-1}{G\′}^{T}d]}^{\frac{M-N}{2}}}{\sqrt{\det ({G\′}^{T}G\′)}}$ 来确定；

G′为矩阵 $\left[\begin{array}{ccc}\cos {\mathrm{\ωx}}_1& \sin {\mathrm{\ωx}}_1& 1\\ \cos {\mathrm{\ωx}}_2& \sin {\mathrm{\ωx}}_2& 1\\ \cos {\mathrm{\ωx}}_3& \sin {\mathrm{\ωx}}_3& 1\\ \·& \·& \\ \·& \·& \\ \·& \·& \\ \cos {\mathrm{\ωx}}_N& \sin {\mathrm{\ωx}}_N& 1\end{array}\right];$

x为该至少一幅输入图像以及与该至少一幅图像相关联的从中获得该矩阵的相应数据的图像之一中的像素位置的标识；

M为G′矩阵中包含的列数目；以及

N为G′矩阵中包含的行数目。

15.如权利要求13所述的方法，还包括：

变换所述输入图像中的该至少一幅的图像像素值，以跨过所述输入图像中的该至少一幅平滑图像强度的正弦变化的幅值，该正弦变化代表格栅图案，其中经变换的值输入到贝叶斯频谱分析。

16.如权利要求15所述的方法，其中所述输入图像中的该至少一幅的图像像素值变换为其相应对数值。

17.如权利要求15所述的方法，其中：

通过应用逆双曲线正弦函数到(q′/2)，变换所述输入图像中的该至少一幅的图像像素值；以及

q′代表未变换像素值。

18.如权利要求13所述的方法，还包括：

确定格栅图案相对于所述输入图像中的该至少一幅的成像区域的倾斜；

旋转所述输入图像中的该至少一幅以消除该倾斜；

对于所述输入图像中的该至少一幅的多个像素行中的每一行，对该行的像素值求和；以及

将该行像素值之和输入贝叶斯频谱分析。

19.如权利要求6所述的方法，其中：

该格栅图案包含在该多幅输入图像的每一个内；

产生第一输出图像包括：

记录该多幅输入图像；

计算角频率ω；

对于所述输入图像中的至少一幅，基于所计算的角频率ω来计算该输入图像的格栅图案的相位角；以及

对于该第一输出图像的每个像素，依据所计算的至少一个相位角并基于该多幅输入图像的相应像素值，计算一值，该方法还包括：

保存所计算的角频率ω；以及

取回所保存的角频率ω用于确定所述方程的值。

20.如权利要求6所述的方法，还包括：

记录该多幅输入图像；

校准在每幅输入图像记录之间移动格栅的致动器，格栅图案由该格栅形成；以及

在校准该致动器的过程中计算角频率ω。

21.如权利要求6所述的方法，还包括：

确定谐波数目m。

22.如权利要求6所述的方法，其中该谐波数目m预设为3。

23.如权利要求1所述的方法，其中产生该第二输出图像包括：

将图像变换应用到该第一输出图像以获得变换数据；

删除代表变换数据的变换图像的预定部分，该变换数据通过删除该预定部分而被修正；以及

基于经修正的变换数据来产生未变换图像。

24.如权利要求23所述的方法，其中产生无压缩图像包括应用逆图像变换。

25.如权利要求23所述的方法，其中该图像变换为小波变换、离散余弦变换(DCT)、离散正弦变换(DST)、离散小波变换(DWT)、快速傅立叶变换(FFT)、阿达玛变换、哈特利变换、Haar小波和Daubechies小波之一。

26.如权利要求23所述的方法，其中该预定部分包括该变换图像的最左侧约1至2％连续像素列的底部约80％。

27.如权利要求1所述的方法，其中该第一输出图像通过光学断层产生。

28.一种其上存储适于由处理器执行的指令的计算机可读取介质，该指令在被执行时致使该处理器执行一种图像产生方法，该图像产生方法包括：

基于多幅输入图像产生第一输出图像；以及

基于该第一输出图像产生第二输出图像，该第二输出图像与减去伪影的该第一输出图像相同。

29.一种成像设备，包括：

相机，用于记录多幅输入图像；以及

处理器，配置成：

基于该多幅输入图像产生第一输出图像；以及

从该第一输出图像除去伪影以产生第二输出图像。

30.如权利要求29所述的成像设备，其中该处理器配置成：确定该伪影对该第一输出图像的图像强度值的贡献；以及为了产生该第二输出图像，从所述图像强度值减去该贡献。

31.如权利要求30所述的成像设备，其中所述图像强度值为像素强度值。

32.如权利要求30所述的成像设备，其中该贡献沿水平和垂直方向之一正弦变化。

33.如权利要求32所述的成像设备，其中该处理器配置成：

为了确定该贡献，确定代表所述正弦变化的方程的系数的值，所述方程的值在像素行和像素列之一之间变化；以及

为了从图像强度值减去该贡献：

通过将所确定的系数代入该方程来确定该方程的值；以及

从图像强度值减去该方程的值。

34.如权利要求33所述的成像设备，其中：

该伪影由所述输入图像中的至少一幅的格栅图案的多个谐波m形成；以及

该方程为a₁cos(ωy)+b₁sin(ωy)+...a_mcos(mωy)+b_msin(mωy)，其中a₁...a_m和b₁...b_m为所述系数，ω为该格栅图案的角频率，以及y为对于每一像素行是不同的像素行号和对于每一像素列是不同的像素列号之一。

35.如权利要求34所述的成像设备，其中：

所述系数通过将回归分析应用到矩阵Q、b和G来确定；

$Q=\left[\begin{array}{c}{q}_1\\ q_2\\ q_3\\ \·\\ \·\\ \·\\ q_x\end{array}\right],$ ,其中q为该第一输出图像的采样像素强度值，x为被采样像素数目；

$b=\left[\begin{array}{c}I\\ a_1\\ b_1\\ \·\\ \·\\ \·\\ a_m\\ b_m\end{array}\right],$ 其中I代表减去该伪影对采样像素强度的贡献的采样像素强度值；以及

$G=\left[\begin{array}{ccccccc}1& \cos \left({\mathrm{\ωy}}_i\right)& \sin \left({\mathrm{\ωy}}_i\right)& \cos \left({2\mathrm{\ωy}}_i\right)& \sin \left(2{\mathrm{\ωy}}_i\right)& \cos \left({3\mathrm{\ωy}}_i\right)& \sin \left({3\mathrm{\ωy}}_i\right)\\ 1& \cos \left({\mathrm{\ωy}}_i\right)& \sin \left({\mathrm{\ωy}}_i\right)& \cos \left({2\mathrm{\ωy}}_i\right)& \sin \left({2\mathrm{\ωy}}_i\right)& \cos \left({3\mathrm{\ωy}}_i\right)& \sin \left(3{\mathrm{\ωy}}_i\right)\\ 1& \cos \left({\mathrm{\ωy}}_i\right)& \sin \left({\mathrm{\ωy}}_i\right)& \cos \left(2{\mathrm{\ωy}}_i\right)& \sin \left({2\mathrm{\ωy}}_i\right)& \cos \left({3\mathrm{\ωy}}_i\right)& \sin \left(3{\mathrm{\ωy}}_i\right)\\ \·& & & & & & \\ \·& & & & & & \\ \·& & & & & & \\ 1& \cos \left({\mathrm{\ωy}}_i\right)& \sin \left({\mathrm{\ωy}}_i\right)& \cos \left({2\mathrm{\ωy}}_i\right)& \sin \left(2{\mathrm{\ωy}}_i\right)& \cos \left({3\mathrm{\ωy}}_i\right)& \sin \left(3{\mathrm{\ωy}}_i\right)\end{array}\right]$

其中G矩阵的每一行与q矩阵的相应值相关联，且其中，对于G矩阵的每一行，y_i为G矩阵行的相应q矩阵值的行号和列号之一。

36.如权利要求35所述的成像设备，其中通过应用公式b＝(G^TG)^-1G^T使用最小二乘方回归来确定所述系数，其中G^T为G的转置矩阵。

37.如权利要求35所述的成像设备，其中：

该处理器配置成变换该第一输出图像的图像像素值以跨过该第一输出图像平滑图像强度的正弦变化的幅值，所述正弦变化代表伪影；

变换的值被采样以用于所述回归分析；

从变换的值减去所述贡献；以及

产生该第二输出图像包括，在从图像强度值减去所述贡献之后，执行图像像素值变换的逆变换。

38.如权利要求37所述的成像设备，其中该第一输出图像的图像像素值变换成其相应对数值。

39.如权利要求37所述的成像设备，其中：

通过应用逆双曲线正弦函数到(q′/2)，变换该第一输出图像的图像像素值；以及

q′代表未变换像素值。

40.如权利要求35所述的成像设备，其中该处理器配置成：

确定该第一输出图像内的伪影相对于该第一输出图像的成像区域的倾斜；

旋转该第一输出图像以消除该倾斜；

对于该第一输出图像的多个像素行中的每一行，对该行的像素值求和；以及

将该行像素值之和输入所述回归分析。

41.如权利要求34所述的成像设备，其中角频率通过贝叶斯频谱分析来确定。

42.如权利要求41所述的成像设备，其中：

角频率通过应用公式 $p\left(\mathrm{\ω}\right|d,I)\∝\frac{{[d^{T}d-d^{T}G\′{({G\′}^{T}G\′)}^{-1}{G\′}^{T}d]}^{\frac{M-N}{2}}}{\sqrt{\det ({G\′}^{T}G\′)}}$ 来确定；

G′为矩阵 $\left[\begin{array}{ccc}\cos {\mathrm{\ωx}}_1& \sin {\mathrm{\ωx}}_1& 1\\ \cos {\mathrm{\ωx}}_2& \sin {\mathrm{\ωx}}_2& 1\\ \cos {\mathrm{\ωx}}_3& \sin {\mathrm{\ωx}}_3& 1\\ \·& \·& \\ \·& \·& \\ \·& \·& \\ \cos {\mathrm{\ωx}}_N& \sin {\mathrm{\ωx}}_N& 1\end{array}\right];$

x为该至少一幅输入图像以及与该至少一幅输入图像相关联的从中获得该矩阵的相应数据的图像之一中的像素位置的标识；

M为G′矩阵中包含的列数目；以及

N为G′矩阵中包含的行数目。

43.如权利要求41所述的成像设备，其中：

该处理器配置成变换所述输入图像中的该至少一幅的图像像素值以跨过所述输入图像中的该至少一幅平滑图像强度的正弦变化的幅值；

所述正弦变化代表格栅图案；以及

将经变换的值输入到贝叶斯频谱分析。

44.如权利要求43所述的成像设备，其中所述输入图像中的该至少一幅的图像像素值变换为其相应对数值。

45.如权利要求43所述的成像设备，其中：

通过应用逆双曲线正弦函数到(q′/2)，变换所述输入图像中的该至少一幅的图像像素值；以及

q′代表未变换像素值。

46.如权利要求41所述的成像设备，其中该处理器配置成：

确定格栅图案相对于所述输入图像中的该至少一幅的成像区域的倾斜；

旋转所述输入图像中的该至少一幅以消除该倾斜；

对于所述输入图像中的该至少一幅的多个像素行的每一行，对该行的像素值求和；以及

将该行像素值之和输入贝叶斯频谱分析。

47.如权利要求34所述的成像设备，其中：

格栅图案包含在所述多幅输入图像中的每一幅内；

产生该第一输出图像包括：

记录所述多幅输入图像；

计算角频率ω；

对于所述输入图像中的至少一幅，基于所计算的角频率ω来计算该输入图像的格栅图案的相位角；以及

对于该第一输出图像的每个像素，依据所计算的至少一个相位角并基于所述多幅输入图像的相应像素值来计算一值；以及该处理器配置成：

保存所计算的角频率ω；以及

取回所保存的角频率ω以用于确定所述方程的值。

48.如权利要求34所述的成像设备，还包括：

灯；

格栅，所述格栅图案通过将该灯的光照射在该格栅上而形成；以及

致动器，用于在每幅输入图像记录之间移动格栅，其中该处理器配置成：

校准该致动器；以及

在校准该致动器的过程中计算角频率ω。

49.如权利要求34所述的成像设备，其中该处理器配置成：

确定谐波数目m。

50.如权利要求34所述的成像设备，其中该谐波数目m预设为3。

51.如权利要求29所述的成像设备，其中，为了产生该第二输出图像，该处理器配置成：

将图像变换应用到该第一输出图像以获得变换数据；

删除代表该变换数据的变换图像的预定部分，该变换数据通过删除该预定部分而被修正；以及

基于经修正的变换数据来产生未变换图像。

52.如权利要求51所述的成像设备，其中产生无压缩图像包括应用逆图像变换。

53.如权利要求51所述的成像设备，其中该图像变换为小波变换、离散余弦变换(DCT)、离散正弦变换(DST)、离散小波变换(DWT)、快速傅立叶变换(FFT)、阿达玛变换、哈特利变换、Haar小波和Daubechies小波之一。

54.如权利要求51所述的成像设备，其中该预定部分包括该变换图像的最左侧约1至2％连续像素列的底部约80％。

55.如权利要求29所述的成像设备，其中该第一输出图像通过光学断层产生。

56.一种用于光学section成像的系统，包括：

相机，用于记录成像表面的多幅输入图像；

格栅；

灯，用于将光照射在该格栅上以将格栅图案投影在该成像表面上，使得每幅输入图像包含位于相应相位角的相应格栅图案；

致动器，用于在每幅输入图像记录之间移动该格栅，使得该多幅输入图像中的至少两幅的格栅图案处于不同相位角；以及

处理器，配置成：

对于该多幅输入图像中的每一幅，计算该图像的格栅图案相位角；

依据该多幅输入图像中的每一幅的相应像素值以及所计算的相位角，通过为该第一输出图像的每个像素计算一值，产生第一输出图像；以及

通过除去包含在该第一输出图像内的伪影来产生第二输出图像，其中该伪影通过下述之一被除去：

(a)：

确定该伪影对该第一输出图像的图像强度值的贡献；以及

从所述图像强度值减去该贡献；和

(b)：

将图像变换应用到该第一输出图像以获得变换数据；

删除代表该变换数据的变换图像的预定部分，该变换数据通过删除该预定部分而被修正；以及

基于经修正的变换数据来产生未变换图像。

Images