CN101660945B - 快速图像重构方法 - Google Patents

Abstract

本发明提出了一种快速图像重构方法，事先通过傅立叶域光学相干层析系统特性对波长信息进行提取，得到一组在波数空间中均匀分布的波长向量，从而得到这个波长向量在CCD中的虚拟位置系数，由离散傅立叶变换补零插值传递函数产生加权系数矩阵，在系统运行时，由加权系数矩阵和采集到的数据进行插值或者对加权系数矩阵进行加窗截断后和采集到的数据进行插值，得到所要求的插值数据。本发明方法简单、易于实现，能够提高傅立叶域光学相干层析数据处理的精度和速度，提高了傅立叶域光学相干层析系统的实时图像重构能力。

Description

快速图像重构方法

技术领域

本发明涉及一种快速图像重构方法，在所述快速图像重构方法中，采用了一种新型的、可用于傅立叶域光学相干层析等有插值需要的仪器的可变插值间隔的插值方法以及基于一切可利用的窗函数进行加窗截断的可变插值间隔的插值方法，从而实现快速图像重构。

背景技术

在快速图像重构领域，傅立叶域光学相干层析系统是一种新型的非接触、高分辨率的光学探测系统，它通过光学干涉的方法对目标进行纵向扫描，最后通过二维或三维重建，得出目标的结构信息、多普勒信息以及偏振信息，因此它可广泛应用于医学成像、各种工业损伤探测之中。傅立叶域光学相干层析技术中，参考光和信号光在光学分光器3中形成干涉，干涉信号由衍射光栅9进行分光，由透镜10聚焦到线扫描CCD11上，由CCD11将模拟信号转换为数字信号，如图2，光谱仪8由衍射光栅9，透镜10和线扫描CCD11构成；CCD11采集到的光栅出射波长按波长线性分布，但数据重建要求波长信息在k空间中线性分布，需要对数据进行插值。傅立叶域光学相干层析系统中可用于快速图像重构的插值方法有多种，如离散傅立叶变换补零插值，B样条拟合，直接线性插值等，而大多数的傅立叶域光学相干层析系统采用离散傅立叶变换补零插值和直接线性插值相结合的方法，即N点数据通过离散傅立叶变换运算之后，产生N点的频率域数据，在高频点补充M*N点的零值，产生M*N+N的数据，然后通过反傅立叶变换得到M*N+N点的数据，其中，M为补零的倍数，最后采用线性插值的方法得到N点的插值数据。由傅立叶域光学相干层析系统扫描采集的一组数据向量

根据传统的离散傅立叶变换补零插的步骤为：

1)将数据进行离散傅立叶变换，得到一组新数据：

$X_1\left(i\right)=\underset{n=0}{\overset{N-1}{\mathrm{\Σ}}}{x}_{n-1}\exp (-j\frac{2\mathrm{\π}}{N}\mathrm{in})$

2)对新数据进行补零插值，得到M倍补零后的数据：

3)对M倍补零后的数据进行离散傅立叶逆变换，得到拓展了M+1倍后的数据；

4)对拓展后的数据根据线性分布K空间进行线性插值，可以得到插值后的数据。

这种方法简单，成熟，但具有运算量大，远远达不到实时图像重构的要求，插值间隔与插值精度由补零倍数M固定，不能随意改变插值间隔等缺点，并且由于进行离散傅立叶变换补零插值后再进行线性插值，因而插值精度也受到了影响，这些局限性严重制约了傅立叶域光学相干层析系统的快速图像重构应用。

发明内容

本发明要解决的技术问题是：克服现有傅立叶域光学相干层析技术的插值方法中插值精度低、运算速度慢，插值间隔固定且不能任意改变等问题，本发明提供一种快速图像重构方法，其中采用了一种新型的插值方法，该方法精度高，运算速度快，并具有可变的插值精度和插值间隔等特点，能有效得提高傅立叶域光学相干层析系统的运算速度和插值精度。

本发明解决其技术问题所采用的技术方案是：用于傅立叶域光学相干层析技术中可变插值间隔的插值方法以及利用一切可利用的窗函数进行加窗截断的可变插值间隔的插值方法。其特征在于以下步骤：

1)衍射光栅9通过透镜10映射在像素点为N的线扫描CCD11上的波长采用光谱仪准确标定，得到一组与CCD11中每个像素对应的均匀分布的波长向量

波长差为Δλ，波长向量在CCD11上的实际位置系数为 $\stackrel{\→}{\mathrm{Index}}1=\{n;n=\mathrm{1,2},......,N\};$

2)由第一个波长λ₁和最后一个波长λ_N，利用公式k＝2π/λ可以得到CCD 11首尾两像素点的波数k′₁＝2π/λ₁和k′_N＝2π/λ_N，利用k′₁和k′_N组成一个长度为N的线性分布的波数向量 ${\stackrel{\→}{k}}^{\′}=\{k_n^{\′}=k_1^{\′}+\frac{k_N^{\′}-k_1^{\′}}{N-1}*(n-1);n=\mathrm{1,2},......,N\},$ 利用公式

反算出相应的波长向量从而利用Δλ可以计算出每个波数k′_n对应的λ′_n在CCD11中的虚拟位置系数 $\stackrel{\→}{\mathrm{Index}}2=\{s_n=\frac{{\mathrm{\λ}}_n^{\′}-{\mathrm{\λ}}_1^{\′}}{\mathrm{\Δ\λ}}+1;n=\mathrm{1,2},......,N\}$ (或者 $\stackrel{\→}{\mathrm{Index}}2=\{s_n=\frac{{\mathrm{\λ}}_n^{\′}-{\mathrm{\λ}}_1^{\′}}{\mathrm{\Δ\λ}};n=\mathrm{1,2},......,N\});$

3)由CCD11采集的数据为实数、以及实信号在离散傅立叶变换中厄米对称的特点，在离散傅立叶变换的过程中，在高频点增加一些数据点，即：

由此，得到一个改进的离散傅立叶变换补零插值传递函数 $\mathrm{TF}(n,s_n)=1+\underset{i=1}{\overset{N/2}{\mathrm{\Σ}}}\cos \left(\frac{2\mathrm{\π}}{N}i(n-s_n)\right),$ 将 $\stackrel{\→}{\mathrm{Index}}1=\{n;n=\mathrm{1,2},......,N\}$ 和 $\stackrel{\→}{\mathrm{Index}}2=\{s_n=\frac{{\mathrm{\λ}}_n^{\′}-{\mathrm{\λ}}_1^{\′}}{\mathrm{\Δ\λ}}+1;n=\mathrm{1,2},......,N\}$ (或者 $\stackrel{\→}{\mathrm{Index}}2=\{s_n=\frac{{\mathrm{\λ}}_n^{\′}-{\mathrm{\λ}}_1^{\′}}{\mathrm{\Δ\λ}};n=\mathrm{1,2},......,N\})$ 中不同的位置n和s_n依次代入

到中产生一个N*N的加权系数矩阵H_N*N(n，s_n)，这样，对插值加权系数的处理完成；

4)傅立叶域光学相干层析系统中的CCD11采集一组纵轴扫描数据向量x＝{x₁，x₂，......，x_N}，进行数据插值，可到插值后的数据x′＝{x_s1，x_s2，......，x_sN}。其计算公式如下

$x^{\′}\left(s_n\right)=\underset{n=1}{\overset{N}{\mathrm{\Σ}}}{x}_n{H}_{N*N}(n,s_n)$

可以用一切可利用的窗函数对加权过程进行截断，根据窗长和 $\stackrel{\→}{\mathrm{Index}}2=\{s_n=\frac{{\mathrm{\λ}}_n^{\′}-{\mathrm{\λ}}_1^{,}}{\mathrm{\Δ\λ}}+1;n=\mathrm{1,2},......,N\}$ (或者 $\stackrel{\→}{\mathrm{Index}}2=\{s_n=\frac{{\mathrm{\λ}}_n^{\′}-{\mathrm{\λ}}_1^{\′}}{\mathrm{\Δ\λ}};n=\mathrm{1,2},......,N\})$ 获得插值起始位置Min和结束位置Max，利用公式：

$x^{\′}\left(s_n\right)=\underset{n=\mathrm{Min}}{\overset{\mathrm{Max}}{\mathrm{\Σ}}}{x}_n{H}_{N*N}(n,s_n)W(n-\mathrm{Min})$

对原始数据进行插值，其中W(*)为所需加窗的窗函数，加窗窗长为L，从而提高新插值方法的运算速度。

在傅立叶域光学相干层析数据处理时，可以用一切可利用的窗函数对加权系数进行截断，将加权系数H_N*N(n，s_n)中的数据进行加窗截断，减少数据处理的数据处理长度和数据处理量。其计算公式如下：

$x^{\′}\left(s_n\right)=\underset{n=\mathrm{Min}}{\overset{\mathrm{Max}}{\mathrm{\Σ}}}{x}_n{H}_{N*N}(n,s_n)W(n-\mathrm{Min})$

公式中W(*)为一切可利用的加窗函数，插值运算起始位置Min和结束位置Max由加窗窗长L和虚拟位置系数 $\stackrel{\→}{\mathrm{Index}}2=\{s_n=\frac{{\mathrm{\λ}}_n^{\′}-{\mathrm{\λ}}_1^{,}}{\mathrm{\Δ\λ}}+1;n=\mathrm{1,2},......,N\}$ (或者 $\stackrel{\→}{\mathrm{Index}}2=\{s_n=\frac{{\mathrm{\λ}}_n^{\′}-{\mathrm{\λ}}_1^{\′}}{\mathrm{\Δ\λ}};n=\mathrm{1,2},......,N\})$ 得到，从而提高可变插值间隔插值运算的处理速度，而且可以将加权系数存储在计算机中，方便运算时调用，避免重复计算。

另外，基于傅立叶变换过程中的能量守恒定律，可以对傅立叶补零插值进行如下修改：

相应地，由此得到改进后的傅立叶插值补零传递函数

$\mathrm{TF}(n,s_n)=1+\underset{i=1}{\overset{N/2}{\mathrm{\Σ}}}\cos \left(\frac{2\mathrm{\π}}{N}i(n-s_n)\right)+(\sqrt{2}-2)\cos \left(\mathrm{\π}(s_n-n)\right),$ 以及相应的加权系数矩阵H_N*N(n，s_n)。

另外，基于傅立叶变换过程中的和相等，也可以对傅立叶补零插值进行如下修改：

相应地，由此得到改进后的傅立叶插值补零传递函数

$\mathrm{TF}(n,s_n)=1+\underset{i=1}{\overset{N/2}{\mathrm{\Σ}}}\cos \left(\frac{2\mathrm{\π}}{N}i(n-s_n)\right)-\cos \left(\mathrm{\π}(s_n-n)\right),$ 以及

相应的加权系数矩阵H_N*N(n，s_n)。

本发明与现有技术相比有如下优点：

1.可在事前提取波长和波数信息，构造与K空间相对应的非线性分布波长向量以及这个波长向量在CCD11像素点所对应的虚拟位置系数向量，从而由传递函数计算出加权系数矩阵H_N*N(n，s_n)，传统离散傅立叶变换补零值插的精度由补零倍数M固定，只能达到1/(M+1)的位置精度，因为虚拟位置系数s_n的位置不是由传统离散傅立叶变换补零插值算法中的补零倍数M固定的，可以为计算机数据精度内的任意一个实数，从而实现可变的插值精度和插值间隔；

2.可利用一切可利用的加窗函数对加权系数矩阵进行加窗截取，并存储到计算机中，方便计算时调用，从而避免重复计算，由传统离散傅立叶变换补零插值方法中有一次N点和一次M*N+N点快速傅立叶变换，因此需要 $\frac{N}{2}{\log }_2\left(N\right)+\frac{(M+1)*N}{2}{\log }_2((M+1)*N)$ 次复数乘法，本方法仅需要N*L次实数乘法，其中，N为CCD11像素点数，L为加窗函数的窗长，从而提高插值方法的处理运算速度，提高离散傅立叶域光学相干层析系统的实时处理能力，由此能够实现快速图像重构。

附图说明

图1为傅立叶域光学相干层析系统数据插值流程图；

图2为傅立叶域光学相干层析系统结构图，其中1为光源，2为光学隔离器，3光学分光器，4为偏振控制器，5为PZT转换器，6为扫描控制器，7为采样目标，8为光谱仪，9为衍射光栅，10为透镜，11为线扫描CCD；

图3为插值效果比较图；

图4为二维图像重构。

具体实施方案

下面结合附图及具体实施方式详细介绍本发明。本实施方案采用的是傅立叶域光学相干层析系统，对其采集到的数据进行插值，其流程如图1所示；具体步骤如下：

(1)本实施方案从图2所示的光谱仪准确标定入射到CCD11的波长，中心波长为849.72nm，光谱分辨率Δλ＝0.0674nm，线扫描CCD11上像素点数N＝2048，CCD11上首尾两个像素点的波长分别为λ₁＝780.7024nm，λ_N＝918.6702nm，各个波长在CCD11中位置系数为 $\stackrel{\→}{\mathrm{Index}}1=\{n;n=\mathrm{1,2},......,N\};$

(2)由

可得CCD11首尾两个像素点所对应的两个波数为

和

计算一个在K空间线性分布的波数向量

${\stackrel{\→}{k}}^{\′}=\{k_n^{\′}=k_1^{\′}+\frac{k_N^{\′}-k_1^{\′}}{N-1}*(n-1);n=\mathrm{1,2},......,N\},$

由这个在K空间线性分布的波数向量通过公式

得到一组非均匀的λ′＝{λ′₁，λ′₂，......，λ′_N}，显然，λ₁＝λ′₁，λ_N＝λ′_N，从而可以利用公式

计算波长λ′＝{λ′₁，λ′₂，......，λ′_N}对应在CCD11中的虚拟位置系数 $\stackrel{\→}{\mathrm{Index}}2=\{s_n=\frac{{\mathrm{\λ}}_n^{\′}-{\mathrm{\λ}}_1^{\′}}{\mathrm{\Δ\λ}}+1;n=\mathrm{1,2},......,N\};$ 另外，关于s_n的计算公式，也可以采用这样并不影响最终的技术效果，当然，相应地，波长λ′＝{λ′₁，λ′₂，......，λ′_N}对应在CCD11中的虚拟位置系数也表示为 $\stackrel{\→}{\mathrm{Index}}2=\{s_n=\frac{{\mathrm{\λ}}_n^{\′}-{\mathrm{\λ}}_1^{\′}}{\mathrm{\Δ\λ}};n=\mathrm{1,2},......,N\};$ 在以下的描述中，为简单、明了起见，将仅以上述第一种计算方式为例，进行描述，但这并不妨碍本发明采用第二种计算方式来实现，事实上，本发明也可以采用其他多种计算方式来实现；

(3)由实际波长在CCD11上位置系数向量

$\stackrel{\→}{\mathrm{Index}}1=\{n;n=\mathrm{1,2},......,N\}$ 和

在CCD11的虚拟位置系数向量

$\stackrel{\→}{\mathrm{Index}}2=\{s_n=\frac{{\mathrm{\λ}}_n^{\′}-{\mathrm{\λ}}_1^{\′}}{\mathrm{\Δ\λ}}+1;n=\mathrm{1,2},......,N\},$

依次抽取不同的n和s_n，由传递函数

得到加权系数矩阵H_N*N(n，s_n)；

(4)由图2所示的傅立叶域光学相干层析系统的CCD11采集到的一组干涉信号数据为x＝{x₁，x₂，......，x_N}，通过加窗长度为L＝11的布莱克曼(Blackman)窗函数 $(W\left(l\right)=0.42+0.5*\cos \left(\frac{2\mathrm{\πl}}{L}\right)+0.08*\cos \left(\frac{4\mathrm{\πl}}{L}\right),$ l∈[0，L-1])截断加权系数，由插值公式得到插值之后的数据，其计算公式如下：

$x^{\′}\left(s_n\right)=\underset{n=\mathrm{Min}}{\overset{\mathrm{Max}}{\mathrm{\Σ}}}{x}_n{H}_{N*N}(n,s_n)W(n-\mathrm{Min})$

s_n由 $\stackrel{\→}{\mathrm{Index}}2=\{s_n=\frac{{\mathrm{\λ}}_n^{\′}-{\mathrm{\λ}}_1^{\′}}{\mathrm{\Δ\λ}}+1;n=\mathrm{1,2},......,N\}$ 给出，

$\mathrm{Max}=s_n+\frac{L-1}{2},$

$\mathrm{Min}=s_n-\frac{L-1}{2}.$

(5)对傅立叶域光学相干层析系统中CCD11采集到数据重复步骤(3)进行插值，由离散傅立叶变换对插值后的各组数据x′(s)进行傅立叶变换得到X′(s)，令对比度Contrast＝6，亮度偏差Brightness_Bias＝-82，对X′(s)中的各个点进行对数运算可以得到图像的灰度值Intensity，其计算公式如下：

Intensity＝Contrast*(10*log₁₀(X′(s)+Brightness_Bias))+255

计算出灰度值还要进行截取，少于0的赋值为0，高于255的赋值为255，其目的是使其灰度范围在[0，255]之间，以符合计算机图像的灰度输出范围，由扫描控制器6对采样目标7重复线扫描，并对CCD11采集到的干涉信号数据中进行插值和映射，重构出二维图像或都三维图像，图4为一张重构的二维图像。

为了对比，在实施本方案的同时，运用传统的离散傅立叶变换补零插值方法，在本次实验中，取一组线扫描数据4倍抽取后再进行插值，插值后的数据如图3所示，与原始数据相减后，新方法均值为0.1409，方差为0.2524，传统的离散傅立叶变换补零方法均值为0.1448，方差为0.2564。可见均值和方差的对比均为基于可变间隔的插值方法的为佳。用傅立叶域光学相干层析系统对一物体进行扫描，采集到2048*300的数据进行处理及重构图像，重构后的图像如图4。在CPU为酷睿Q9300，内存为4GB时，运算时间由原来的9秒减少到目前的400毫秒，大大提高了处理速度。

基于傅立叶变换过程中的能量守恒定律，可以对傅立叶补零插值进行如下修改：

相应地，由此得到改进后的傅立叶插值补零传递函数

$\mathrm{TF}(n,s_n)=1+\underset{i=1}{\overset{N/2}{\mathrm{\Σ}}}\cos \left(\frac{2\mathrm{\π}}{N}i(n-s_n)\right)+(\sqrt{2}-2)\cos \left(\mathrm{\π}(s_n-n)\right),$ 以及

相应的加权系数矩阵H_N*N(n，s_n)。

基于傅立叶变换过程中的和相等，也可以对傅立叶补零插值进行如下修改：

相应地，由此得到改进后的傅立叶插值补零传递函数

$\mathrm{TF}(n,s_n)=1+\underset{i=1}{\overset{N/2}{\mathrm{\Σ}}}\cos \left(\frac{2\mathrm{\π}}{N}i(n-s_n)\right)-\cos \left(\mathrm{\π}(s_n-n)\right),$ 以及

相应的加权系数矩阵H_N*N(n，s_n)。

虽然通过参照本方法的说明和具体实施方案，已经对本方法进行了图示和描述，但普通的技术人员应该明白，可以在形式上和细节上对其作各种各样的改变，而不偏离所附权利要求书所限定的本方法的精神和范围。

Claims (8)

1.一种快速图像重构方法，包括以下步骤：

对衍射光栅(9)通过透镜(10)映射在像素点为N的线扫描CCD(11)上的波长进行准确标定，得到一组在波长空间中均匀分布的波长向量波长差为Δλ，波长向量在CCD(11)上的实际位置系数为

将在波数空间中、均匀分布在与最短波长λ₁相对应的波数和与最长波长λ_N相对应的波数之间的波数向量转换为一组在波长空间中非均匀分布的波长向量

并根据上述波长差Δλ，得到波长向量

的各个波长在CCD(11)中的虚拟位置系数 $\stackrel{\→}{\mathrm{Index}}2=\{s_n;n=\mathrm{1,2},......,N\};$

根据补零插值传递函数TF(n，s_n)，由

和

产生一个N*N的加权系数矩阵H_N*N(n，s_n)；

在CCD(11)采集到一组数据向量x＝{x₁，x₂，……，x_N}时，利用加权系数矩阵H_N*N(n，s_n)，对数据向量x＝{x₁，x₂，……，x_N}进行插值处理，得到插值后的数据 $x^{\′}\left(s_n\right)=\{x_{s_1}^{\′},x_{s_2}^{\′},......,x_{s_N}^{\′}\};$ 以及

对插值后的数据进行离散傅立叶变换，实现图像重构，

其中

对插值后的数据

进行离散傅立叶变换后，得到强度相关数据X′(s_n)，并根据下述公式实现图像重构：

Intensity＝Contrast*(10*log₁₀(X′(s_n)+Brightness_Bias))+255，

其中Intensity表示灰度值，Contrast表示对比度，Brightness_Bias表示亮度偏差，

如果得到的灰度值Intensity小于0，则为Intensity赋值0；以及如果得到的灰度值Intensity大于255，则为Intensity赋值255。

2.根据权利要求1所述的快速图像重构方法，其中所述补零插值传递函数是根据下述公式得到的：

$\mathrm{TF}(n,s_n)=1+\underset{i=1}{\overset{N/2}{\mathrm{\Σ}}}\cos \left(\frac{2\mathrm{\π}}{N}i(n-s_n)\right).$

3.根据权利要求1所述的快速图像重构方法，其中所述补零插值传递函数是根据下述公式得到的：

$\mathrm{TF}(n,s_n)=1+\underset{i=1}{\overset{N/2}{\mathrm{\Σ}}}\cos \left(\frac{2\mathrm{\π}}{N}i(n-s_n)\right)+(\sqrt{2}-2)\cos \left(\mathrm{\π}(s_n-n)\right).$

4.根据权利要求1所述的快速图像重构方法，其中所述补零插值传递函数是根据下述公式得到的：

$\mathrm{TF}(n,s_n)=1+\underset{i=1}{\overset{N/2}{\mathrm{\Σ}}}\cos \left(\frac{2\mathrm{\π}}{N}i(n-s_n)\right)-\cos \left(\mathrm{\π}(s_n-n)\right).$

5.根据权利要求1所述的快速图像重构方法，其中所述虚拟位置系数

是根据下述公式得到的：

$s_n=\frac{{\mathrm{\λ}}_n^{\′}-{\mathrm{\λ}}_1^{\′}}{\mathrm{\Δ\λ}}+1.$

6.根据权利要求1所述的快速图像重构方法，其中所述虚拟位置系数

是根据下述公式得到的：

$s_n=\frac{{\mathrm{\λ}}_n^{\′}-{\mathrm{\λ}}_1^{\′}}{\mathrm{\Δ\λ}}.$

7.根据权利要求1所述的快速图像重构方法，其中

除了利用加权系数矩阵H_N*N(n，s_n)外，还同时利用窗函数对加权系数矩阵H_N*N(n，s_n)进行加窗截断后对数据向量x＝{x₁，x₂，……，x_N}进行插值处理。

8.根据权利要求7所述的快速图像重构方法，其中

对数据向量x＝{x₁，x₂，……，x_N}的插值处理是根据下述公式执行的：

$x^{\′}\left(s_n\right)=\underset{n=\mathrm{Min}}{\overset{\mathrm{Max}}{\mathrm{\Σ}}}{x}_n{H}_{N*N}(n,s_n)W(n-\mathrm{Min})$

其中W(*)是窗函数，Max和Min的值由窗函数的加窗长度和各个虚拟位置共同确定。

Images