CN107564068B - 一种针对孔径编码超分辨光学传递函数的标定方法 - Google Patents

Abstract

本发明公开了一种针对孔径编码超分辨光学传递函数的标定方法，首先利用计算机编程生成编码矩阵，然后将得到的编码矩阵作为孔径编码图案，在LCOS空间光调制器上显示，在已知物体亚像素移动步长的情况下，利用孔径迭代算法标定每个编码矩阵对应的OTF，最后在物体亚像素移动步长未知的情况下，标定出物体的亚像素位移。本发明可以精确标定出物体的亚像素位移，故可以保证重构出的光学系统实际OTF的准确性，从而避免由于光学系统误差导致的超分辨率重构质量下降。

Description

一种针对孔径编码超分辨光学传递函数的标定方法

技术领域

本发明属于光学成像技术，特别是一种针对孔径编码超分辨光学传递函数的标定方法。

背景技术

虽然现有的光电成像技术已取得巨大进步，并获得一系列的优异成果，然而，无论是高价专业数码相机还是低廉手机摄像头，都未突破传统的“小孔成像”方式，仍为基于透镜式“所见即所得”的成像模式，虽然这种模式原理简单，操作易行，但是仍然存在难以逾越的局限，即图像像素化——探测器尺寸已达物理极限。现有的光电成像系统，在信息获取、具体功能、性能指标等方面，皆受限于探测器工艺水平及制作成本。通常而言，解决图像像素化问题的最直接方法即为提高探测器阵列密度以提高分辨率。但实际上，由于受到探测器制作工艺的限制，目前探测器尺寸水平已达物理极限，因此若想再通过此法来减少图像像素化问题，是非常困难的。同时，即使在不考虑成本的前提下实现了像元尺寸的进一步减小，然而由于光通量下降带来了灵敏度降低的问题，使得原本微弱的光信号愈加掺杂在严重的噪声中，导致信噪比更低，后期更加难于补偿。

为了突破由CCD探测器几何尺寸造成的分辨率限制，有研究人员提出了多种超分辨成像技术(代少升,张德洲,崔俊杰,等.基于微扫描的红外超分辨率成像系统的设计[J].半导体光电,2017,38(1):103-106.)。超分辨率成像技术能在不改变设备硬件条件的情况下，利用多帧图像间的互补信息提高图像空间分辨率,常见的方法是亚像素扫描超分辨成像法(Peleg S,Keren D,Schweitzer L.Improving image resolution using subpixelmotion[J].Pattern recognition letters,1987,5(3):223-226.)，它可实现多帧同一场景下互有亚像素级位移图像的采集，从而优化最终成像质量。但这种方法需要额外的运动部件或是摆镜，系统十分复杂；如此一来，实验条件要求非常苛刻，实验可操作性差。2005年，Solomon J等人提出了在成像系统的傅里叶平面放置一个掩模，这个掩模对物频谱进行编码成像之后再对像频谱进行解码(Solomon J,Zalevsky Z,Mendlovic D.Geometricsuper resolutionby code division multiplexing[J].Applied optics,2005,44(1):32-40.)。加掩模的方法可以克服由CCD两相邻素中心间距离引起的频谱混叠问题。但该方法忽略了CCD像素大小，将CCD像素看成理想的点，并没有解决由CCD每个像素的大小和形状引起的低通效应问题。并且由于其通过反卷积的方法解出物体的频谱，所以该方法必须精确知道系统的OTF((光学传递函数))，且对拍得的图质量要求比较苛刻(刘晶丹,许廷发,荀显超,等.光学掩模实现几何超分辨成像的仿真[J].光学精密工程,2014,22(8):2026-2031.)。2013年，刘海英等人提出了利用数字微镜阵列(DMD)实现超分辨重构，但是DMD器件存在“衍射光栅效应”，影响超分辨效果(刘海英,李云松,吴成柯.一种数字微镜阵列分区控制和超分辨重建的压缩感知成像法[J].光子学报,2013,43(5):510002-0510002.)。所以，想要更加高效高质地实现超分辨率成像又不引入了机械移动装置，需要作进一步研究。

发明内容

本发明的目的在于提供一种针对孔径编码超分辨光学传递函数的标定方法，既可以通过精确移动物体而拍得的多幅图来重构OTF，又能通过方法准确标定出物体的亚像素位移，具有非常好的稳定性、灵活性和可操作性。

实现本发明目的的技术解决方案为：一种针对孔径编码超分辨光学传递函数的标定方法，步骤如下：

第一步，利用计算机编程生成编码矩阵，其中编码矩阵的生成公式为：

R＝Ascend(N_r,N_c)

R_sort＝Randsort(R)

R_mod＝Mod(R_sort,L²)

其中，R为N_r行N_c列递增矩阵，其矩阵元素为1～N_r×N_c，Ascend(·)为产生连续递增矩阵的函数,N_r表示所使用的编码孔径的个数，N_c为每次编码孔径中同时置1的LCOS像素个数，R_sort为R经过随机排序后的N_r行N_c列的矩阵，Randsort(·)为产生随机排列的函数，R_mod为矩阵R_sort中的每个元素对L²求余后得到的矩阵，Mod(·)为求余操作，L为LCOS一行/一列的像素个数，W表示编码矩阵，W(k,l)为矩阵W中第k行l列的元素；

第二步，将得到的编码矩阵作为孔径编码图案，在LCOS空间光调制器上显示，在已知物体亚像素移动步长的情况下，利用孔径迭代方法标定每个编码矩阵对应的OTF；

第三步，在物体亚像素移动步长未知的情况下，标定出物体的亚像素位移。

本发明与现有技术相比，其显著优点：(1)在物体亚像素移动步长未知的情况下，可以精确得标定出物体的亚像素位移，故可以保证重构出的光学系统实际OTF的准确性。(2)可以准确重构出的光学系统实际的OTF，从而避免由于光学系统误差导致的超分辨率重构质量下降。(3)利用计算机编程生成编码矩阵，然后将得到的编码矩阵作为孔径编码图案，在LCOS空间光调制器上显示，在已知物体亚像素移动步长的情况下，利用孔径迭代算法标定每个编码矩阵对应的OTF，解决了基于孔径编码的超分辨成像方法中需要知道不同孔径编码对应的OTF的难题。

下面结合附图对本发明作进一步详细描述。

附图说明

图1为本发明针对孔径编码超分辨传递函数的标定方法的流程图。

图2为孔径迭代步骤的流程图。

图3为标定亚像素位移步骤的流程图。

图4为不同F#下的理论的光学系统的传递函数。其中，图4的(a)中F#为5.6，图4的(b)中F#为14，图4的(c)中F#为22。

图5为孔径迭代算法迭代出的不同F#下实际的光学系统的传递函数。其中，图5的(a)中F#为5.6，图5的(b)中F#为14，图5的(c)中F#为22。

图6的(a)拍得的低分辨率，图6的(b)为用理论的传递函数进行孔径编码超分辨的结果图，图6的(c)为用孔径迭代算法迭出的传递函数进行孔径编码超分辨的结果图。

具体实施方式

结合图1，本发明针对孔径编码超分辨光学传递函数的标定方法，步骤如下：

第一步，利用计算机编程生成编码矩阵。利用计算机Matlab软件(或其他类似编程软件)编程生成编码矩阵。编码矩阵的生成公式为：

R＝Ascend(N_r,N_c)

R_sort＝Randsort(R)

R_mod＝Mod(R_sort,L²)

其中，R为N_r行N_c列递增矩阵，其矩阵元素为1～N_r×N_c，Ascend(·)为产生连续递增矩阵的函数,N_r表示所使用的编码孔径的个数，N_c为每次编码孔径中同时置1的LCOS(硅基液晶)像素个数，R_sort为R经过随机排序后的N_r行N_c列的矩阵，Randsort(·)为产生随机排列的函数。R_mod为矩阵R_sort中的每个元素对L²求余后得到的矩阵，Mod(·)为求余操作，L为LCOS一行/一列的像素个数，W表示编码矩阵，W(k,l)为矩阵W中第k行l列的元素。

第二步，将得到的编码矩阵作为孔径编码图案，在LCOS空间光调制器上显示，在已知物体亚像素移动步长的情况下，利用孔径迭代算法标定每个编码矩阵对应的OTF，具体步骤如下：

1.在频域中生成一个高分辨率光学传递函数的初始解P₀，一般选择元素值全为0的矩阵初始化高分辨率光学传递函数。

2.选择相对于所用的超分辨成像光学系统可看成点光源的物体作为目标物o，在水平和竖直方向亚像素平移目标物o得到平移后的目标物o_m,n,同时拍摄一组对应的低分辨率光强图

直到在水平和垂直方向都平移了一个像素。其中o_m,n为第m次水平位移第n次竖直位移后的目标物，

为第m次水平位移第n次竖直位移后拍得的低分辨率光强图，上标c表示拍摄到的图像。m,n＝1,2,...N,N≥Mag,Mag为超分辨的倍数。

3.对平移后的目标物体o_m,n进行傅里叶变换得到目标物体光强的频谱O_m,n，利用传递函数P_i截取频谱O_m,n中某一子孔径里的频谱信息，生成一个目标光强分布

目标光强生成公式为

下标i表示第i次更新，上标e表示目标图像，F^-1{·}为傅里叶逆变换算子，|·|为取模运算。

4.对目标光强图

进行降采样Mag倍得到

并生成对应的系数矩阵C，生成公式为

其中上标ds表示降采样。

5.更新目标光强图

更新公式为

其中，C^us为最邻近插值Mag倍后的系数矩阵，上标us表示上采样，上标u表示更新后的图。

6.更新光学传递函数P_i,更新公式为

其中，max{·}表示取矩阵元素最大值,F{·}为傅里叶变换算子。

7.如果还有拍得低分辨率光强图未被用于更新光学传递函数,则重复迭代步骤3～6来利用其他低分辨率光强图更新光学传递函数；

8.当所有低分辨率光强图都被用于更新后，再重复迭代步骤3～7，直到高分辨率光学传递函数收敛，从而获得高分辨率光学传递函数最优解。

图4为不同F#下的理想的光学系统的传递函数。其中，图4的(a)中F#为5.6，图4的(b)中F#为14，图4的(c)中F#为22。图5为孔径迭代算法迭代出的不同F#下的光学系统的传递函数。其中，图5的(a)中F#为5.6，图5的(b)中F#为14，图5的(c)中F#为22。可以看出(1)传递函数从低频到高频呈衰减趋势，且F数越小，传递函数衰减的速度越快。(2)理想传递函数和实际测得的传递函数略有不同。

第三步，在物体亚像素移动步长未知的情况下，标定出物体的亚像素位移，具体步骤如下：

1.对参考图r_m×n和待标定图g_m×n做离散傅里叶变换，得到参考图的频谱R_m×n和待标定图的频谱G_m×n，生成插值2倍后的参考图和待标定图的互相关矩阵C，生成公式为

R_m×n＝DFT{r_m×n}

G_m×n＝DFT{g_m×n}

其中，下标m、2m为对应矩阵的行数，n、2n为对应矩阵的列数，DFT{·}为离散傅里叶变换算子，IDFT{·}为离散逆傅里叶变换算子，P_m×n为参考图频谱与待标定图频谱的共轭的乘积矩阵，上标*表示取共轭运算，

为乘积矩阵P_m×n在四周填0后得到的矩阵。

2.计算参考图的半高md2和半宽nd2，并找到互相关矩阵C的最大值所在行rloc所在列cloc，生成亚像素移动步长的初始估计，生成公式为

其中，rshift为竖直方向亚像素移动的距离，cshift为竖直方向亚像素移动的距离。此时，得到的亚像素移动距离已精确到0.5个像素。

3.设所需配准精度为1/usfac个像素，利用上采样离散傅里叶变换在互相关矩阵C中以(cshift,rshift)为中心，面积为1.5倍像素×1.5倍像素的邻域U中进行插值，插值公式如下

nr＝nc＝ceil(usfac*1.5)

colvec＝-fix(nc/2):ceil(nc/2)

rowvec＝-fix(nr/2):ceil(nr/2)

N＝kernr·P_m×n·kernc

其中，i为虚数单位，nr表示邻域U插值usfac倍后的行数，nc表示邻域U插值usfac倍后的列数，ceil(·)表示向正无穷方向的取整，fix(·)为向零取整，colvec为一nc维行向量，rowvec为一nr维行向量，双目运算符：生成其左操作数到右操作数的整数行向量，kernc,kernr均为上采样离散傅里叶变换核，上标T表示矩阵的转置，ifftshift将矩阵的第一象限和第三象限互换，第二象限和第四象限互换，V为领域U插值usfac倍后的得到的矩阵，符号“·”表示矩阵乘法运算。

4.找到矩阵V最大值所在行rloc2和所在列cloc2，计算出待配准图相对于参考图亚像素平移距离的精确解，计算公式为

center＝fix(ceil(usfac*1.5)/2)+1

rshift＝rshift+(rloc2-center)/usfac

cshift＝cshift+(cloc2-center)/usfac

其中，center为矩阵V中心所在的行和列。

通过上述步骤，本发明可以精确得标定出物体的亚像素位移，故可以保证重构出的光学系统实际OTF的准确性，从而避免由于光学系统误差导致的超分辨率重构质量下降。

图6的(a)拍得的低分辨率，图6的(b)为用理论的传递函数进行孔径编码超分辨的结果图，图6的(c)为用孔径迭代算法迭出的传递函数进行孔径编码超分辨的结果图。可以看出，与直接利用理论值进行孔径编码超分辨相比，用本方法迭代出的传递函数进行孔径编码超分辨有着更好的效果。

Claims (1)

1.一种针对孔径编码超分辨光学传递函数的标定方法，其特征在于步骤如下：

第一步，利用计算机编程生成编码矩阵，其中编码矩阵的生成公式为：

R＝Ascend(N_r,N_c)

R_sort＝Randsort(R)

R_mod＝Mod(R_sort,L²)

其中，R为N_r行N_c列递增矩阵，其矩阵元素为1～N_r×N_c，Ascend(·)为产生连续递增矩阵的函数,N_r表示所使用的编码孔径的个数，N_c为每次编码孔径中同时置1的LCOS像素个数，R_sort为R经过随机排序后的N_r行N_c列的矩阵，Randsort(·)为产生随机排列的函数，R_mod为矩阵R_sort中的每个元素对L²求余后得到的矩阵，Mod(·)为求余操作，L为LCOS一行/一列的像素个数，W表示编码矩阵，W(k,l)为矩阵W中第k行l列的元素；

第二步，将得到的编码矩阵作为孔径编码图案，在LCOS空间光调制器上显示，在已知物体亚像素移动步长的情况下，利用孔径迭代方法标定每个编码矩阵对应的OTF；

第三步，在物体亚像素移动步长未知的情况下，标定出物体的亚像素位移；

其中，在第二步中，利用孔径迭代步骤如下：

(1)在频域中生成一个高分辨率光学传递函数的初始解P₀；

(2)选择相对于所用的超分辨成像光学系统可看成点光源的物体作为目标物o，在水平和竖直方向亚像素平移目标物o得到平移后的目标物o_m,n,同时拍摄一组对应的低分辨率光强图

直到在水平和垂直方向都平移了一个像素；其中o_m,n为第m次水平位移第n次竖直位移后的目标物，

为第m次水平位移第n次竖直位移后拍得的低分辨率光强图，上标c表示拍摄到的图像，m,n＝1,2,...N,N≥Mag,Mag为超分辨的倍数；

(3)对平移后的目标物体o_m,n进行傅里叶变换得到目标物体光强的频谱O_m,n，利用传递函数P_i截取频谱O_m,n中某一子孔径里的频谱信息，生成一个目标光强分布

目标光强生成公式为

下标i表示第i次更新，上标e表示目标图像，F^-1{·}为傅里叶逆变换算子，|·|为取模运算；

(4)对目标光强图

进行降采样Mag倍得到

并生成对应的系数矩阵C，生成公式为

其中上标ds表示降采样；

(5)更新目标光强图

更新公式为

其中，C^us为最邻近插值Mag倍后的系数矩阵，上标us表示上采样，上标u表示更新后的图；

(6)更新光学传递函数P_i,更新公式为

其中，max{·}表示取矩阵元素最大值,F{·}为傅里叶变换算子；

(7)如果还有拍得低分辨率光强图未被用于更新光学传递函数,则重复迭代步骤(3)～(6)来利用其他低分辨率光强图更新光学传递函数；

(8)当所有低分辨率光强图都被用于更新后，再重复迭代步骤(3)～(7)，直到高分辨率光学传递函数收敛，从而获得高分辨率光学传递函数最优解；

所述第三步中，标定出物体的亚像素位移的步骤为：

(1)对参考图r_m×n和待标定图g_m×n做离散傅里叶变换，得到参考图的频谱R_m×n和待标定图的频谱G_m×n，生成插值2倍后的参考图和待标定图的互相关矩阵C，生成公式为

R_m×n＝DFT{r_m×n}

G_m×n＝DFT{g_m×n}

其中，下标m、2m为对应矩阵的行数，n、2n为对应矩阵的列数，DFT{·}为离散傅里叶变换算子，IDFT{·}为离散逆傅里叶变换算子，P_m×n为参考图频谱与待标定图频谱的共轭的乘积矩阵，上标*表示取共轭运算，

为乘积矩阵P_m×n在四周填0后得到的矩阵；

(2)计算参考图的半高md2和半宽nd2，并找到互相关矩阵C的最大值所在行rloc所在列cloc，生成亚像素移动步长的初始估计，生成公式为

其中，rshift为竖直方向亚像素移动的距离，cshift为竖直方向亚像素移动的距离；

(3)设所需配准精度为1/usfac个像素，利用上采样离散傅里叶变换在互相关矩阵C中以(cshift,rshift)为中心，面积为1.5倍像素×1.5倍像素的邻域U中进行插值，插值公式如下

nr＝nc＝ceil(usfac*1.5)

colvec＝-fix(nc/2):ceil(nc/2)

rowvec＝-fix(nr/2):ceil(nr/2)

N＝kernr·P_m×n·kernc

其中，i为虚数单位，nr表示邻域U插值usfac倍后的行数，nc表示邻域U插值usfac倍后的列数，ceil(·)表示向正无穷方向的取整，fix(·)为向零取整，colvec为一nc维行向量，rowvec为一nr维行向量，双目运算符：生成其左操作数到右操作数的整数行向量，kernc,kernr均为上采样离散傅里叶变换核，上标T表示矩阵的转置，ifftshift将矩阵的第一象限和第三象限互换，第二象限和第四象限互换，V为领域U插值usfac倍后的得到的矩阵，符号·表示矩阵乘法运算；

(4)找到矩阵V最大值所在行rloc2和所在列cloc2，计算出待配准图相对于参考图亚像素平移距离的精确解，计算公式为

center＝fix(ceil(usfac*1.5)/2)+1

rshift＝rshift+(rloc2-center)/usfac

cshift＝cshift+(cloc2-center)/usfac

其中，center为矩阵V中心所在的行和列。

Images