用于快照光谱成像的设备和方法

摘要

本发明涉及用于在快照中获得来源对象的多幅光谱图像的设备和方法，其中包括使用被应用于发散-漫射快照图像的二维压缩感知数据立方体重建(2D？CS-SCR)。在一些实施例中，所述快照图像是通过RIP漫射器获得的。在一些实施例中，使用随机化器来对发散-漫射快照图像进行进一步随机化。2D？CS-SCR包括对表示从快照图像导出的光谱立方体数据的不同波段的各个阵列单独应用2D小标架变换。对表示不同波段的各个阵列单独应用2D小标架变换包括对各个阵列应用直接和逆2D小标架变换。在一些实施例中，所述直接和逆小标架变换被包括在分裂Bregman迭代中。

说明书

相关申请的交叉引用

本申请涉及2013年6月18日提交的具有相同标题的美国临时专利申请号61/836, 218并且要求其优先权，其被全文合并在此以作参考。

技术领域

这里所公开的实施例总体上涉及光谱成像(SI)设备、系统和方法，并且特别涉及 利用数字摄影机和基于压缩感知(CS)的空间-光谱立方体重建(SCR)算法的快照(“单个镜 头”)超光谱成像(HSI)。

背景技术

“光谱成像”通常被理解成指的是利用给定波长范围内的有限数目的波长(例如多 达10个波长)的成像，所述给定波长范围例如是可见(“VIS”)范围或近红外(“NIR”)范围。 “超光谱成像”通常被理解成指的是利用给定波长范围内的数目更多的波长的成像，例如10 个到数百个或数千个之间的波长。“快照”SI和HIS成像器在单个快照中实施对于空间和光 谱数据的同时(瞬时)采集。所采集的数据形成来源对象(其在这里也被简称作“对象”或“场 景”)的“空间-光谱立方体”(其在这里也被简称作“光谱立方体”或“数据立方体”)。“空间- 光谱立方体”、“光谱立方体”和“数据立方体”在后文中可以互换使用。数据立方体包括两个 空间维度和一个光谱维度中的光强度数据，并且被表达成三维(3D)矩阵。

标题为“Snapshotspectralimagingbasedondigitalcameras(基于数字摄 影机的快照光谱成像)”的共同著作和转让的美国专利申请号13/752,56(其被公开为美国 专利公开号20130194481)被全文合并在此以作参考，其教导了包括插入在来源图像与像素 化(例如在数字摄影机中)图像传感器之间的光径中的成像透镜、发散图像传感器和受限等 距属性(RIP)漫射器的设备中的基于压缩感知(CS)的快照光谱成像(CS-SSI)。RIP漫射器可 以是一维的(1D)。其在发散图像传感器处提供经过发散和漫射的图像(“DD图像”)。由于1D RIP漫射器光学属性，DD图像中的每一个像素包括来自DD图像中的相应一列的所有像素的 光谱和空间信息的线性混合。在US20130194481中，利用基于CS的优化处理来实施数据立 方体的完全重建，以便补偿所述问题的欠定性质(underdeterminednature)。实施线性射 影的运算符可以被描述成“感知矩阵”，其行数少于列数并且在数据立方体上进行运算以形 成DD图像。如果感知矩阵满足RIP条件，则所述重建处理确保对于来源对象的完全重建。RIP 漫射器被设计成使得包括所述漫射器的光学成像系统的转移函数(其与所述感知矩阵完全 相同)在每一个单个波长处(或者在围绕某一单个波长选择的波段处)都满足RIP条件。

在US20130194481中提供的解决方案利用块Toeplitz矩阵实施1DCS-SCR，以便 实施顺序地应用于包括在垂直方向上串联的所有波段图像的一个阵列当中的各列的单一 1D变换。已经发现，在希望重建的信号所要求的稀疏度方面，对应于块Toeplitz矩阵的RIP 条件比对应于随机矩阵的RIP条件更难满足。

发明内容

在各个实施例中公开了针对基于CS原理的HIS的方法和设备。在后面的描述中， “快照光谱成像器”、“SSI设备”和“SSI摄影机”可以互换使用。虽然这里所公开的设备和方 法也可以被用于具有更少数目的波长的光谱成像，但是所提到的始终是“HSI”。因此，在后 面的描述中，“光谱成像”可以被使用作为利用从三个波长(即红色(R)、绿色(G)和蓝色(B)) 到数百个以及甚至数千个波长的成像的一般术语。

这里所公开的一种设备是快照超广谱成像器。因此，“用于HIS的设备”、“HIS设备” 和“快照超光谱成像器”可以互换使用。虽然集中在对于来源对象的光谱成像上，但是这里 所公开的方法和设备也可以被应用于在对象上的不同空间位置处具有多种光谱的灰度或 彩色图像。

我们发现，通过对于表示不同波段的各个阵列单独应用二维(2D)小标架 (framelet)变换而不是这里所使用的顺序地成列应用1D变换，可以进一步改进在共同拥有 的US20130194481中所描述的SCR。对于表示光谱立方体数据的不同波段的各个阵列单独 应用2D小标架变换在后文中被称作“2DCS-SCR”，其中包括对各个阵列应用直接和逆2D小 标架变换。直接和逆小标架变换被示例性地包括在分裂Bregman迭代中。Bregman迭代内部 的小标架变换使用分裂到波长(split-to-wavelength)波段而不是分裂到一个空间坐标末 端。这里所公开的2DCS-SCR提供了比US20130194481中所描述的1DCS-SCR更快并且更好 的(在峰值信噪比或“PSNR”方面)的SCR。

我们还发现，对于特定对象，通过在图像采集期间添加利用所添加的硬件(HW)光 学元件(“HW随机化器”)实施或者通过软件以算法方式(“SW随机化器”)实施的随机化运算， 可以进一步改进2DCS-SCR结果。在后面所公开的一些设备和方法实施例中，随机化器被添 加到US20130194481中所描述的SSI设备。随机化器帮助重建非稀疏(“常规”)对象的光谱 图像。HW随机化器可以被实施成图像传感器平面处的纤薄光学元件。随机化器使得到达图 像传感器的DD图像数据(对于HW随机化器)或者利用图像传感器获得的DD图像数据(对于SW 随机化器)变成完全随机，并且得到“随机化图像”。这里所使用的“随机化图像”指的是在随 机化器对于DD图像的动作之后获得的图像数据。在一些实施例中，从统计总体 (statisticalensemble)取得随机化器的单个随机矩阵R。

在一些实施例中，提供了具有最小硬件改变的基于数字摄影机的快照光谱成像设 备和方法。SSI处理包括在具有或不具有所添加的随机化的情况下从来源对象的DD图像实 施2DCS-SCR。

在一个实例中，提供了用于在快照中获得来源对象的多幅光谱图像的设备，所述 设备包括：包括透镜和像素化图像传感器的数字摄影机的成像部分，所述成像部分被配置 成获得DD快照图像Y；以及数字处理器，其被配置成从快照图像Y实施2DCS-SCR，从而提供 多个光谱带中的来源对象的图像。

在一个实例中，提供了一种用于在快照中获得来源对象的多幅光谱图像的方法， 其包括获得DD快照图像Y以及从快照图像Y实施2DCS-SCR的步骤，从而提供多个光谱带中 的来源对象的图像。在一些实施例中，通过利用包括透镜和位于图像传感器平面处的像素 化图像传感器的数字摄影机的成像部分对来源对象进行成像而获得快照图像Y，其中所述 DD图像是通过满足与感知矩阵A有关的RIP条件的RIP漫射器形成的。在这样的实施例中，2D CS-SCR包括把感知矩阵A转置成转置矩阵A^T，将A^T应用于Y以获得A^TY，对A^TY应用2D直接稀疏 化变换D以获得重建数据立方体X的稀疏版本d，利用逆变换Ψ从d获得X，以及对X进行处理 以获得所述多个光谱带中的来源对象的图像。

在一些实施例中，利用基于样条的标架从随机化图像重建光谱图像。在一些实施 例中，对来自几幅随机化单色图像的叠加的重建光谱图像施加基于样条的标架。

附图说明

通过结合附图考虑后面的详细描述，这里所公开的各个方面、实施例和特征将变 得显而易见。在不同的附图中可以用相同的附图标记来标示相同的元件，其中：

图1A示意性地示出了基于具有RIP漫射器并且没有随机化器的数字摄影机的这里 所公开的快照光谱成像(SSI)设备的一个实施例；

图1B示出了图1A的设备的更多细节；

图2A示意性地示出了基于具有RIP漫射器和SW随机化器的数字摄影机的这里所公 开的另一个SSI设备的另一个实施例；

图2B示出了图2A的设备的更多细节；

图2C示意性地示出了基于具有RIP漫射器和HW随机化器的数字摄影机的这里所公 开的另一个SSI设备的另一个实施例；

图2D示出了图2C的设备的更多细节；

图3A示出了对应于RIP漫射器设计的一个实施例的锯齿相位的曲线图；

图3B示出了对应于RIP漫射器设计的另一个实施例的像素的排列；

图3C示出了基于图2B的RIP漫射器的设计；

图3D示出了图3C的RIP漫射器设计的灰度表示；

图4示出了具有16x16像素的示例性随机化器设计；

图5A-D示出了四通道PR滤波器组的脉冲和量值响应的曲线图，其对于APZFrame的 一个示例性实施例实施了Bregman迭代内部的小标架变换；

图6A在流程图中示意性地示出了这里所公开的用于SSI的方法的一个实施例；

图6B示出了图6A的流程图中的步骤508的细节；

图6C示出了对应于分裂Bregman迭代的图6B中的步骤522的细节；

图7A-J示出了在利用装备有1DRIP漫射器的数字摄影机获得的快照图像上进行 的2DCC-SCR仿真的结果；

图8A-J示出了在利用装备有1DRIP漫射器和软件随机化器的数字摄影机获得的 快照图像上进行的2DCC-SCR仿真的结果；

图9A示意性地示出了这里所公开的SSI设备的另一个实施例；

图9B示意性地示出了这里所公开的SSI设备的另一个实施例；

图10示意性地示出了这里所公开的SSI设备的另一个实施例；

图11示意性地示出了这里所公开的SSI设备的另一个实施例；

图12示意性地示出了这里所公开的SSI设备的另一个实施例；

图13示意性地示出了这里所公开的SSI设备的另一个实施例。

具体实施方式

不具有随机化器的SSI设备

图1A和1B示意性地示出了基于具有RIP漫射器的数字摄影机的快照光谱成像 (SSI)设备的一个实施例100。设备100可以被视为光学耦合到位于摄影机与利用摄影机成 像的来源对象(图1B中的106)之间的RIP漫射器104的数字摄影机102。摄影机102可以包括 常规数字摄影机的组件，比如透镜(未示出)、图像传感器108和数字处理器110。所述数字摄 影机可以是任何已知的摄影机，比如SLR摄影机、视频摄影机、蜂窝电话摄影机或者使用在 胃窥镜或内窥镜中的类型的小型化摄影机。这里所使用的“数字摄影机”可以仅包括在像素 化图像传感器上获得来源对象的图像所需要的成像部分的基本元件和功能。因此，例如存 在于智能电话摄影机中的耦合到图像处理器的透镜/传感器组合应当被视为用于这里所阐 述的目的的成像部分。

设备100在原理上类似于US20130194481中所公开的设备(例如其中的设备200)， 其不同之处在于处理器110被配置成实施2DCS-SCR而不是US20130194481中所公开的1D CS-SCR。后面提供2DCS-SCR处理的一个详细的实现实例。可选的是，设备100可以包括被配 置成实施这里所公开的2DCS-SCR的一部分或全部的附加的外部(摄影机外部)数字处理器 105。

图1B提供了设备100的更加详细的视图。除了前面的组件之外，摄影机100包括具 有孔径114的成像透镜112以及可选地还有带通滤波器116。RIP漫射器104可以位于成像透 镜112的入射光瞳处或其附近。正如在US20130194481中详细描述的那样，对象106通过RIP 漫射器104和透镜112的成像提供图像传感器108处的快照DD图像118。与图像118相关联的 图像数据由处理器110处理，其可以被配置成还实现其他功能，比如控制各种摄影机功能和 操作。

光学系统的数学模型

将参照如图1和2中的单透镜SSI设备来描述下面的模型。如果摄影机的整个光学 系统被视为具有从所有透镜的组合导出的有效焦距的单个透镜，则这样的设备的基本原理 同样适用于多透镜摄影机。所述模型具有若干阶段：

将漫射图像表示成原始图像与点扩散函数(PSF)之间的卷积

假设所获得的来源对象的理想原始图像(在不使用漫射器或随机化器的情况下) 具有强度分布I₀(x，y；λ)，其是波长λ处的数据立方体的剖面。RIP漫射器具有复传递函数P：

其中，是波长λ处的漫射器的相位函数。当被安装到光学系统中时，RIP 漫射器把原始图像转换成DD图像，这是因为成像系统不再是理想的。DD图像的形状和特性 可以作为P和原始图像的函数来计算。所述系统的相干点扩散函数可以被计算成P的傅立叶 变换：

$h(y^{\prime };\lambda )=\frac{1}{i\lambda R}\int P({\upsilon }^{\prime };\lambda )\exp \left(i2\pi \frac{{\upsilon }^{\prime }{y}^{\prime }}{\lambda R}\right){\mathrm{d\upsilon }}^{\prime }---\left(2\right)$

并且其在电磁场的复幅度中描述了对应于输入处的δ函数的所述系统的脉冲响 应，其中R是从出射光瞳到成像系统的图像传感器的距离。如果光是不相干的，则可以仅测 量由图像传感器接收到的光的强度。相应地，通过由下式给出的不相干PSFh_I(y′；λ)来描 述系统在强度中的脉冲响应：

h_I(y′；λ)＝|h(y′；λ)|²λ²(3)

空间偏移不变模型成像系统把DD图像强度提供成理想(“非发散”)图像I与不相干 PSFh_I之间的1D卷积

I′(x,y′；λ)＝∫h_I(y′-y；λ)I(x,y；λ)dy(4)

其中，I(x，y；λ_l)是在不使用漫射器或随机化器的情况下由成像系统在波长λ_l处获 得的空间不相干来源对象的理想图像的强度，并且x,y是图像传感器处的笛卡尔坐标。应当 提到的是，对于理想图像的每一个坐标x单独计算1D卷积。

将DD或随机化图像、数据立方体和PSF表示成矩阵

由于DD图像是利用像素化图像传感器取得的，因此其实际上被采样，并且可以被 表示成每一个像素中的强度的矩阵。不相干PSF还可以通过表示卷积等式(4)的Toeplitz矩 阵来表示。图像传感器自然地具有离散像素化结构，其由2D空间间距δ_x×δ_y、像素数目N_x，N_y以及每个像素的比特数目N_b表征。在一个实施例中，成像变焦被选择成使得由RIP漫射器导 致的图像模糊使得发散-漫射图像在图像传感器处占据每一列中的所有N_y个像素以及每一 行中的所有N_x个像素。相应地，在没有RIP漫射器的情况下在相同的变焦下获得的“未漫射- 未发散图像”在图像传感器处仅占据位于每一列的中心部分处的更少数目的N＜N_y个像素 以及每一行中的所有N_x个像素。数据立方体被定义成尺寸为N_x×N×L的3D阵列，其中N_x、N是 空间维度，L是光谱维度，也就是光谱图像中的光谱带或波长的数目。尽管所感知到的像素 的数目N_x×N_y(即实验数据的维度)可能远小于具有维度N_x×N×L的目标3D数据立方体中的 体素的数目，但是我们建议了一种通过求诸于CS方法并且利用图像数据中的隐含的冗余性 的针对3D数据立方体的解决方案。所建议的解决方案提供了数据压缩率N×L/N_y。

遵循CS中的离散标记管理，我们定义以下索引范围：在数据立方体和传感器像素 的y方向上延伸的范围和在数据立方体和传感器像素的x方向上延伸 的范围以及在数据立方体的λ方向上延伸的范围传感器的像素的中 心具有以下笛卡尔坐标：

$x_j=(j-\frac{N_x+1}{2}){\delta }_x,y_{i^{\prime }}=-(i^{\prime }-\frac{N_y+1}{2}){\delta }_y---\left(5\right)$

数据立方体的体素共享传感器的空间间距但是具有不同的索引范围，因此其索引 被如下偏移：

$i_c=\frac{N_y-N}{2}>0$

并且其笛卡尔坐标是x_j和RIP漫射器的平面处的笛卡尔坐标被标记成u′， v′。具有透明孔径的线性维度D_u′×D_v′的RIP漫射器(参见图3B和3C)包括：

N_d＝D_v′/Δυ′(6)

即平行于u′轴延伸的垂直直线条带，其具有宽度Δv′和中心：

$\begin{array}{cc}{\upsilon }_k^{\prime }=(k-\frac{N_d+1}{2}){\mathrm{\Delta \upsilon }}^{\prime }& k=\overline{1,N_d}\end{array}---\left(7\right)$

因此，可以通过仅取决于坐标v′的复分段恒定光瞳函数来描述RIP漫射器：

$\begin{array}{cc}P({\upsilon }^{\prime };{\lambda }_l)=\underset{k=0}{\overset{N_d}{\Sigma }}{P}_{kl}rect\left(\frac{{\upsilon }^{\prime }-{\upsilon }_k^{\prime }}{{\mathrm{\Delta \upsilon }}^{\prime }}\right)& l=\overline{1,L}\end{array}---\left(8\right)$

其中，是RIP漫射器上的第k个条带的宽度Δv′内的相位常数，是光谱带编号的中心波长，并且L是光谱带的总数目。

对应于不相干PSF的等式(2)和(3)提供离散卷积内核以作为对应于每一个波长的 Toeplitz卷积矩阵：

$h_I(y_{i^{\prime }}-y_{i+i_c};{\lambda }_l)={\left(\frac{D_{{\upsilon }^{\prime }}}{\lambda R}\right)}^2{K}_{i^{\prime }-i,l}---\left(10\right)$

其中，卷积内核为如下：

$\begin{array}{c}{K}_{{\mathrm{\Delta i}}^{\prime },l}={\mathrm{sinc}}^2\left(\frac{{\mathrm{\Delta i}}^{\prime }-i_c}{N_{\lambda }}\right)\frac{1}{N_d^2}\left|\underset{k=0}{\overset{N_d}{\Sigma }}{P}_{kl}\exp \right[-i\frac{2\pi }{N_{\lambda }}(k-\frac{N_d+1}{2})({\mathrm{\Delta i}}^{\prime }-i_c)]{|}^2\\ {\mathrm{\Delta i}}^{\prime }=i^{\prime }-i,i=\overline{1,N},i^{\prime }=\overline{1,N_y},l=\overline{1,L}\end{array}---\left(11\right)$

$N_{\lambda }=\frac{\lambda R}{{\mathrm{\Delta \upsilon }}^{\prime }{\delta }_y}---\left(12\right)$

并且P_kl由等式(9)定义。应当提到的是，对应于固定的的阵列K_Δi′，l具有 范围因此包括N_y+N-1个元素。数据立方体的体素可以被表 达成：

$X_{i,l}^{\left(j\right)}=I(x_j,y_{i+i_c};{\lambda }_l),i=\overline{1,N},j=\overline{1,N_x},l=\overline{1,L}---\left(13\right)$

数据立方体的每一个光谱带中的理想图像强度的离散版本由一个N_x×N ×L阵列表示：

$X=(X_{i,l}^{\left(j\right)},i=\overline{1,N},j=\overline{1,N_x})---\left(14\right)$

换句话说，X是表示空间-光谱数据立方体的矩阵。

假设光学系统仅允许1D发散，从而使得图像的两个空间维度x,y不会被混合，并且 DD图像的每一列可以被单独考虑。每一列包括图像数据(对应于图像矩阵)和相应的转移函 数(PSF矩阵)。此外，由于发散仅仅是1D的，因此PSF的各列完全相同，从而允许丢弃对应于 PSF的列索引j。因此，在每一个波长处，连续1D卷积的等式(4)可以被重写成对于N_x个图像 列当中的每一个单独应用的离散1D卷积。具有单个波长λ_l的光对于DD图像的各个离散像素 的贡献可以通过离散形式被如下表达：

$I(x_j,y_{i^{\prime }};{\lambda }_l)={\left(\frac{D_{{\upsilon }^{\prime }}}{\lambda R}\right)}^2{\delta }_y\underset{i=1}{\overset{N}{\Sigma }}{K}_{i^{\prime }-i,l}{X}_{i,j,l},i^{\prime }=\overline{1,N_y},j=\overline{1,N_x},l=\overline{1,L}---\left(15\right)$

其中，j是DD图像中以及数据立方体中的列编号，并且K_i′-i，l是Toeplitz“卷积矩 阵”的元素。等式(15)表明，在单个光谱带中，由成像透镜形成在图像传感器上的光强度被 描述成通过等式(10)和(11)定义的数据立方体与Toeplitz矩阵的各个元素的离散卷积。

具有随机化器的SSI设备

图2A示意性地示出了基于数字摄影机202的另一个SSI设备的编号为200的一个实 施例，所述数字摄影机202具有RIP漫射器104以及通过配置作为数字摄影机202的一部分的 数字处理器210而实施的SW随机化器220。在图2B中示出了设备200的更多细节。设备实施例 100和200所共有的元件由相同的附图标记标示，其具有类似的功能并且因此不作详细描 述。设备200还包括SW随机化器220。SW随机化器220是操作在图像传感器平面处的软件实施 的随机化器“功能”。所述随机化器在图像传感器108上提供经过漫射、发散和随机化的图像 118’。数字处理器210’被配置成在随机化图像上实施2D-CSSCR。随后从经过漫射、发散和 随机化的图像118’重建光谱图像。

图2C示意性地示出了基于具有RIP漫射器和HW随机化器的数字摄影机的另一个 SSI设备的编号为200’的一个实施例。在图2D中示出了设备200’的更多细节。设备实施例 100和200’所共有的元件由相同的附图标记标示，其具有类似的功能并且因此不作详细描 述。设备200’包括HW随机化器220’。所述随机化器在图像传感器108上提供随机化图像 204’。数字处理器210’被配置成在随机化图像上实施2D-CSSCR。随后从随机化图像重建光 谱图像。

图2的实施例中的RIP漫射器被显示成位于成像透镜之前(在从来源对象到图像传 感器的光径中)。在其他实施例中，RIP漫射器可以位于成像透镜内部或之后，并且优选地位 于成像透镜的入射光瞳处。

随机化器的贡献

随机化器对初始测量矩阵的Toeplitz结构进行随机化，并且对于合理地稀疏的真 实图像允许重建。随机化器可以被表示成具有随机元素R_i′，j的2D矩阵，其中其 与波长无关，并且仅影响传感器上的信号的幅度。其对于每一个波长处的所测量的光强度 具有相同的效果。这一效果是每一个DD图像像素处的光强度I′(x_j，y_i′；λ_l)乘以相应的元素 R_i′，j以获得随机化的DD图像。随机化器的插入通过添加逐个像素的乘法改动等式(15)：

$R_{i^{\prime },j}{I}^{\prime }(x_j,y_{i^{\prime }};{\lambda }_l),i=\overline{1,N},i^{\prime }=\overline{1,N_y},l=\overline{1,L},j=\overline{1,N_x}---\left(16\right)$

其中，对于DD图像的每一列仅使用随机化矩阵$R=(R_{i^{\prime },j},i^{\prime }=\overline{1,N_y},j=\overline{1,N_x})$的单个列$R=(R_{i^{\prime },j},i^{\prime }=\overline{1,N_y}).$通过设定<OK>，这里所开发的数学等式同样适用于不 具有随机化器的系统：

$R_{i^{\prime },j}\equiv 1,i=\overline{1,N},i^{\prime }=\overline{1,N_y},l=\overline{1,L},j=\overline{1,N_x}---\left(17\right)$

在一个实施例中，随机化器可以被实施成用于照片或视频摄影机的数字处理器或 者处于摄影机膝上型或台式计算机外部的数字处理器的算法和软件代码。在另一个实施例 中，随机化器可以被实施成硬件，特别是被实施成位于照片或视频摄影机的成像透镜与图 像传感器之间的光学元件，其优选地紧邻图像传感器或者被安放在图像传感器上。

在用于光谱成像的CS中，具有整个波长集合的光对于DD图像的各个离散像素的贡 献被标示成并且可以被表达成在每一个图像传感器像素处的所有波长上的DD图像的 强度的总和，从而获得感知强度：

$Y_{i^{\prime }}^{\left(j\right)}={\left(\frac{\lambda R}{D_{{\upsilon }^{\prime }}}\right)}^2\frac{1}{{\delta }_y}\underset{l=1}{\overset{L}{\Sigma }}{\kappa }_l{R}_{i^{\prime },j}{I}^{\prime }(x_j,y_{i^{\prime }};{\lambda }_l)=\underset{l=1}{\overset{L}{\Sigma }}\underset{i=1}{\overset{N}{\Sigma }}{A}_{i^{\prime },i,l}^{\left(j\right)}{X}_{i,l}^{\left(j\right)}---\left(18\right)$

非负数字κ_l(在我们的计算机仿真中低于κ_l＝1)表征波长λ_l处的图像传感器的相 对光谱灵敏度，并且系数描述RIP漫射器与随机化器R_i′，j的组合效果：

$A_{i^{\prime },i,l}^{\left(j\right)}=R_{i^{\prime },j}{B}_{i^{\prime },i,l},i=\overline{1,N},i^{\prime }=\overline{1,N_y},l=\overline{1,L},j=\overline{1,N_x}---\left(19\right)$

$B_{i^{\prime },i,l}={\kappa }_l{K}_{i^{\prime }-i,l},i=\overline{1,N},i^{\prime }=\overline{1,N_y},l=\overline{1,L}---\left(20\right)$

因此，随机化器打破每一个波长处的感知矩阵的Toeplitz结构，从而产生对应于 信号的甚至更加随机的结构。其随机地改变由每一个图像传感器像素接收到的幅度，从而 改进满足RIP条件的能力。

在具有随机化器的单个版本(“单项随机化”动作)的一个实施例中，图像传感器像 素处的所采集的DD图像像素与随机化器的灰度像素的乘法可能使得其中一些图像传感器 像素实际上未被使用。举例来说，随机化器像素的零或接近零的数值将使得对应的图像传 感器像素表现为零或接近零，即实际上丢失。使用这样的随机化器可能导致光通量的一些 损失，以及SI摄影机的吞吐量和灵敏度的降低。但是后者对于具有RIP漫射器和单色图像传 感器的SI摄影机中的噪声减少是非常重要的。通过使用后面所描述的SW随机化器的多个版 本(其也被称作“多项随机化”)，可以实现这里所公开的SI摄影机的高吞吐量和灵敏度。

压缩感知公式

出于数学方面的考虑因素，串联数据立方体中的光谱和垂直空间维度看起来是方 便的，也就是说在阵列和中用单个索引替换两个索引i,l。相应地，我们求诸于 具有扩大的长度NL的1D矢量X^(j)：

$X^{\left(j\right)}=(X_{i,l}^{\left(j\right)},i=\overline{1,N},l=\overline{1,L})=\left(\begin{array}{c}\left(\begin{array}{c}{X}_{1,j,1}\\ .\\ .\\ .\\ X_{N,j,1}\end{array}\right){\lambda }_1\\ .\\ .\\ .\\ \left(\begin{array}{c}{X}_{1,j,L}\\ .\\ .\\ .\\ X_{N,j,L}\end{array}\right){\lambda }_L\end{array}\right),j=\overline{1,N_x}---\left(21\right)$

可以把N_x个列矢量的整个集合合并到一个矩阵中：

$X=[X^{\left(j\right)},j=\overline{1,N_x}]=[X^{\left(1\right)},X^{\left(2\right)},...,X^{\left(N_x\right)}]---\left(22\right)$

该矩阵的尺寸为NL×N_x，并且包含所有光谱立方体的数据。矩阵X可以被替换地分 解到L个光谱维度中：

$X=\left(\begin{array}{c}{X}_1\\ X_2\\ .\\ .\\ .\\ X_L\end{array}\right)---\left(23\right)$

从而使得每一个光谱维度由一个尺寸为N×N_x的子矩阵X_l描述。一个N_y×NL维感知 矩阵

$A=(A_{i^{\prime },i,l},i^{\prime }=\overline{1,N_y,}i=\overline{1,N,}l=\overline{1,L})---\left(24\right)$

可以作为一个分块矩形矩阵A＝[A₁,A₂,…A_L]来对待，其由L个尺寸分别为N_y×N的 子矩阵A_l构成。每一个子矩阵A_l对应于单个波长。矩阵A通过集成RIP漫射器和随机化器(当 其存在时)的效果提供了高度随机化。我们还定义DD图像的矢量：

$Y^{\left(j\right)}=(Y_{i^{\prime }}^{\left(j\right)},i^{\prime }=\overline{1,N_y})=\left(\begin{array}{c}{Y}_1^{\left(j\right)}\\ .\\ .\\ .\\ Y_{N_y}^{\left(j\right)}\end{array}\right),j=\overline{1,N_x}---\left(25\right)$

应当提到的是，矢量X^(j)是我们希望从一列感知矢量Y^(j)重建的数据。可以把N_x个 列矢量的整个集合合并到具有尺寸N_y×N_x的感知强度矩阵中：

$Y=[Y^{\left(j\right)},j=\overline{1,N_x}]=[Y^{\left(1\right)},Y^{\left(2\right)},...,Y^{\left(N_x\right)}]$

矩阵X是我们希望从感知强度矩阵(或快照图像)Y重建的光谱数据。等式(17)现在 可以通过矩阵形式被表达成具有长度NL的矢量与具有维度N_y×NL的矩阵的乘法。所述乘法 得到具有更小长度N_y的矢量：

Y^(j)＝AX^(j)(26)

相应地，对于2D数据处理，合并的矢量和矩阵可以通过矩阵形式被表达成尺寸为 N_y×NL的矩阵与具有维度NL×N_x的矩阵的乘法，从而得到具有更小尺寸N_y×N_x的矩阵：

Y＝AX.(27)

等式(26)提供了对应于每一列j的CS模型，并且等式(27)提供了对应于图像传感 器处的整个二维DD图像的CS模型。CS问题在于重建矩阵X从而使得对于给定的矩阵Y满足等 式(27)。

由于来源对象的可压缩性属性，可以在该对象在其中是稀疏的空间中稀疏地表示 该对象。通过实施包括来源对象的重建矢量乘以感知矩阵与发散图像之间的差异的l₂范数 的函数的最小化，可以从发散图像重建来源对象的稀疏表示。因此，所述对象在其中也是稀 疏的空间中的坐标的l₁范数。所述最小化处理(在所述约束下)可以示例性地通过分裂 Bregman迭代处理[Z.Cai等人，SIAMJ.ofMultiscaleModelingandSimulation，Vol.8 (2)，pp.337-369，2009年，后文中称作“Cai”]来实现。已经知道该处理是用于CS重建的高效 的工具。分裂Bregman迭代是涉及闭合环路的迭代算法，其中受到重建约束的l₁误差充当针 对环路的反馈，并且具有确保稀疏重建的收缩操作。

由于A是N_y×NL矩阵并且L>1，因此未知变量的数目NL×N_x大于等式的数目N_y×N_x。 因此，所述问题似乎是不适定(ill-posed)的，并且在通常情况下无法被解决。但是CS感知 理论解决了一种特定情况，也就是当矩阵X是可压缩的并且从而可以在某种(可能是冗余 的)基础上被表示成K稀疏矩阵d的线性变换，其中d＝DX是仅具有已知位置处的K个非零元 素的矩阵，并且D是稀疏化矩阵。稀疏化矩阵是把矢量或阵列转换成稀疏矩阵的矩阵，也就 是说转换成仅具有少数非零元素的矩阵。所述冗余基础可以示例性地通过求诸于2D小标架 变换来实施，正如在后面的“Spline-basedframesforspectralimagereconstruction (用于光谱图像重建的基于样条的标架)”章节中更加详细地描述的那样。由此，我们在基于 CS的光谱成像中第一次把源自二次准内插多项式样条的2D半紧密小波标架(或小标架)应 用于光谱图像重建。关于基于样条的标架的详细描述，并且特别是关于从内插和准内插多 项式样条开发多种低通滤波器h⁰的详细描述，可以在A.Averbuch和V.Zheludev的 “Interpolatoryframesinsignalspace(信号空间中的内插标架)”(IEEE Trans.Sign.Proc.，54(6)，pp.2126-2139，2006年)中找到。关于小标架的描述可以在 A.Averbuch、P.Neittaanmaki和V.Zheludev的“Splinesandsplinewaveletmethods withapplicationtosignalandimageprocessing.VolumeI:Periodicsplines(样 条和样条小波方法对于信号和图像处理的应用。第I卷：周期性样条)”(Springer，2014年) (后文中称作“APZFrame”)中找到。具体来说，我们利用稀疏化矩阵D应用2D直接线性变换， 从而获得数据立方体X的稀疏版本d。在一个实施例中，稀疏化矩阵D可以是直接2D小标架变 换的矩阵，其被单独地应用于矩阵X的每一个子矩阵X_l，l＝1,…,L：

$d=DX\iff \left(\begin{array}{c}{\mathrm{DX}}_1\\ {\mathrm{DX}}_2\\ \mathrm{..}\\ \mathrm{..}\\ \mathrm{..}\\ {\mathrm{DX}}_L\end{array}\right)=\left(\begin{array}{c}{d}_1\\ d_2\\ .\\ .\\ .\\ d_L\end{array}\right).---\left(28\right)$

我们利用矩阵Ψ应用2D逆线性(作为示例是标架)变换，以便从通过K稀疏矩阵d描 述的其稀疏版本获得数据立方体X：

$X=\Psi d\iff \left(\begin{array}{c}{\mathrm{\Psi d}}_1\\ {\mathrm{\Psi d}}_2\\ .\\ .\\ .\\ {\mathrm{\Psi d}}_L\end{array}\right)=\left(\begin{array}{c}{X}_1\\ X_2\\ \mathrm{..}\\ \mathrm{..}\\ \mathrm{..}\\ X_L\end{array}\right).---\left(29\right)$

其是对于矩阵d的每一个子矩阵d_l(l＝1,…,L)单独应用的逆2D小标架变换的矩 阵。现在可以通过以下形式来表达Y：

Y＝AX＝AΨd＝Θd(30)

其中，Θ＝AΨ。众所周知，CS中的K阶的RIP条件(参见E.J.Candes和T.Tao的 “DecodingbyLinearProgramming(通过线性编程进行解码)”，IEEETrans.Information Theory51(12)：4203-4215(2005年))要求通过零化其所有各列而形成的Θ的任何子矩阵 (除了小于K阶的之外)对于任何K稀疏矢量d必须满足以下不等式：

$(1-{\delta }_K)\left|\right|d|{|}_{l_2}\le |\left|\Theta d\right|{|}_{l_2}\le (1+{\delta }_K)\left|\right|d|{|}_{l_2}---\left(31\right)$

其中，δ_K＞0是某一更小的数字，并且已经知道，即使对于减 少的数目N_y(从而使得K＜N_y＜NL)的射影，对于满足K阶的RIP条件的矩阵Θ＝AΨ也存在并 且可以找到等式(28)的稳定解d(并且因此还有X)。

满足RIP条件的感知矩阵的其中一个最为人所知的实例是由具有零均值和1/NL方 差的高斯随机变量形成的随机矩阵或随机Toeplitz矩阵。在这种情况下，如果满足以下不 等式，则各列近似正交并且以高概率满足RIP条件：

$N_y\ge C\frac{K}{log\left(\frac{NL}{K}\right)}---\left(32\right)$

其中，0<C≤1是常数。通过求解下面的受到约束的最小化问题，我们从给定的矩阵 Y重建光谱立方体X：

其中DX是如在等式(25)中的矩阵X的分块2D小标架变换，矢量a的l₁范数是 并且σ是预定的停止阈值。

在Cai中给出了求解最小化问题等式(33)的一种方法。所述处理通过引入被单独 处理的几个附加的变量来工作。更具体来说，遵循这里所实施的分析，通过迭代处理实施对 应于线性运算符A的最小化：

$\left(\begin{array}{c}{X}^{k+1}={({\mathrm{\mu A}}^{T}A+\chi I)}^{-1}[{\mathrm{\mu A}}^T(-c^k)+{\mathrm{\chi D}}^T(d^k-b^k)]\\ d^{k+1}=shrink({\mathrm{DX}}^{k+1}+b^k,\frac{1}{\chi })\\ b^{k+1}=b^k+{\delta }_b({\mathrm{DX}}^{k+1}-d^{k+1})\\ c^{k+1}=c^k+{\delta }_c({\mathrm{AX}}^{k+1}-Y)\end{array}\right)---\left(34\right)$

其中，k是迭代次数，d^k和c^k是被用来执行迭代的中间矢量，A^T标示A的转置矩阵，并 且

shrink(x,γ)＝sgn(x)max(|x|-γ,0)(35)

是被应用于矩阵的每一个分量的函数。所述处理的参数μ、χ(其中χ^-1是收缩阈值) 允许为问题中的各项给出不同的重要性或权重：||AX^k+1-Y||_l、|DX^k+1|_l。通过为d的稀疏度水 平指派不同的权重(更大的χ)或者更好地拟合到AX＝Y(更大的μ)，所述参数的改变允许实 质上收敛到解。如果感知矩阵满足在等式(28)中表达的RIP条件，则所述重建处理确保对于 来源对象的令人满意的重建。一旦实现不等式||AX^k+1-Y||_l≤σ(其中σ是停止阈值)，则终止 迭代。

在迭代完成之后，我们具有对应于重建的数据立方体的可压缩分块矩阵：

X＝Ψd(36)

其包括子矩阵X_l＝Ψd_l(l＝1,…,L)，其中稀疏分块矩阵d＝DX。随后通过具有尺寸 N×N_x的矩阵来表示包含每一个光谱带中的重建图像的重建数据立方体分量X_l＝Ψd_l：

$\left(\begin{array}{ccc}{X}_{1,l}^{\left(1\right)}& \mathrm{...}& X_{1,l}^{\left(N_x\right)}\\ .& & .\\ .& & .\\ .& & .\\ X_{Nl}^{\left(1\right)}& \mathrm{...}& X_{Nl}^{\left(N_x\right)}\end{array}\right)$

其中，是相应的光谱带的索引。可以在每一个光谱带处通过计算PSNR来 实施重建的X与理想的图像之间的比较：

${\mathrm{PSNR}}_l=20log\left(\frac{max{\bar{X}}_{i,l}^{\left(j\right)}}{\sqrt{\frac{1}{N_{x}N}\underset{i=1}{\overset{N}{\Sigma }}\underset{j=1}{\overset{N_x}{\Sigma }}{[X_{i,l}^{\left(j\right)}-{\bar{X}}_{i,l}^{\left(j\right)}]}^2}}\right),l=\overline{1,L}---\left(37\right)$

具有锯齿衍射光栅的排列的RIP漫射器设计

在US20130194481中描述了一种示例性1DRIP漫射器实现方式。但是在数学方 面，方便的做法是通过光瞳放大系数将其缩放到出射光瞳的尺寸。US20130194481中的RIP 漫射器的沟槽是针对具有如图9B中示出的相位水平的特定波长λ_des设计的，并且被实施成 微起伏(micro-relief)水平线。编号为k的每一条线由宽度Δv′和随机沟槽深度h_k表征，正 如在图9C和9D中示例性地示出的那样。正如在US20130194481中更加详细地描述的那样， 作为1D元件，RIP漫射器变成具有随机沟槽深度的垂直直线条带的阵列。

RIP漫射器包括平行于u′轴延伸的N_d个垂直直线条带，其具有宽度Δv′和在等式 (7)中定义的中心v′_k。沟槽深度h_k在第k个条带的宽度Δv′内是恒定的。每一个沟槽深度h在 波长λ处导致通过下面的近轴情况等式给出的相应的相位偏移

其中，n(λ)是波长λ处的折射率。由于相位是与波长有关的，因此每一个沟槽深度 为具有不同波长的光添加了不同的相位。对应于两个不同波长的相位添加通过下式相关：

前面的等式中的近似可以被应用是因为折射率n随着波长缓慢地变化。因此，如果 掩模沟槽是针对特定波长λ_des设计的，则掩模对于具有波长λ的光的影响是：

根据等式(36)，由RIP漫射器提供的相位可以被如下描述：

其中，是RIP漫射器上的编号为k的直线条带处的设计波长λ_des处的相位，λ_l是编号为l的光谱带的中心波长，并且L是光谱带的总数目。

与RIP漫射器相关联的相干点扩散函数h(y′；λ)也是1D的，其仅取决于图像平面处 的坐标y′，并且可以被计算成分段恒定光瞳函数的逆傅立叶变换。求诸于作为sinc函数的 rect函数的傅立叶变换的已知结果：

$\sin c\left(\xi \right)=\frac{sin\left(\pi \xi \right)}{\pi \xi }---\left(41\right)$

并且求诸于偏移属性，从而得到：

$\begin{array}{c}h(y^{\prime };{\lambda }_l)=\frac{1}{iR}\underset{k=1}{\overset{N_d}{\Sigma }}{P}_{kl}\underset{-0.5D_y}{\overset{0.5D_y}{\int }}rect\left(\frac{{\upsilon }^{\prime }-{\upsilon }_k^{\prime }}{{\mathrm{\Delta \upsilon }}^{\prime }}\right)\exp \left(i2\pi \frac{{\upsilon }^{\prime }{y}^{\prime }}{\lambda R}\right){\mathrm{d\upsilon }}^{\prime }\\ =\frac{D_{{\upsilon }^{\prime }}}{iR}\sin c\left(\frac{{\mathrm{\Delta \upsilon }}^{\prime }}{\lambda R}{y}^{\prime }\right)\frac{1}{N_d}\underset{k=0}{\overset{N_d}{\Sigma }}{P}_{kl}\exp \left(i2{\mathrm{\pi \upsilon }}_k^{\prime }\frac{y^{\prime }}{\lambda R}\right)\end{array}---\left(42\right)$

其中，v′_k和Δv′分别是1DRIP漫射器中的第k个直线条带的位置和宽度，P_kl在RIP 漫射器上的第k个条带的宽度Δv′内是恒定的并且在等式(9)中通过相位偏移定义， λ_l是编号为的光谱带的中心波长。因此，(对应于幅度的)PSF是以传 感器中心为中心的各个sinc函数的总和，所述各个sinc函数具有不同的与漫射器相关的系 数，以及由漫射器上的恒定相位线的倒宽度定义的宽度。

在这里所公开的一个实施例中，还如下开发RIP漫射器设计。RIP漫射器可以被安 装在成像系统的入射光瞳处或其附近。但是在数学等式中，便利的做法是通过光瞳放大系 数将其缩放到出射光瞳的尺寸。在一个实施例中，RIP漫射器是1D纤薄相位光学元件，其在 单个操作方向上提供光学光场的相位改变，并且包括垂直于操作方向延伸的线沟槽。RIP漫 射器由具有折射率n(λ)的透明材料制作，并且由平行于u′延伸的N_d个垂直直线条带构成。 深度和相位在第k个条带的宽度Δv′内是恒定的，并且按照2π/N_Q的相位差被彼此等距地量 化到N_Q个离散相位水平。在这里给出的一个实施例中，对应于RIP漫射器的设计开始于具有 锯齿轮廓的闪耀衍射光栅，正如图3A中所示出的那样。闪耀衍射光栅在v′轴的方向上具有 周期Λ，并且相应地具有周期性锯齿相位函数：

πmod_2π(2πυ′/Λ)(43)其中，mod_2π(.)函数表示在减去2π的倍数之后的自变量 的最小正残差。所述相位函数被量化到N_Q个离散相位水平，从而使得条带宽度是Δv′＝Λ/ N_Q。条带的总数被如下选择：

N_d＝D_y/Δυ′＝N_QD_y/Λ(44)

在一个实施例中，产生经过量化的锯齿阵列，其在每一个循环中具有对应于沟槽 深度水平的数目的一定数目的N_Q个点，并且具有总数目为N_d个像素。因此，所述锯齿阵列中 的每一个点k代表对应于一个条带的相位数值：

在该实施例中，沟槽深度水平的数目N_Q和闪耀周期Λ受到实际考虑因素的限制， 也就是针对特征尺寸的制作规则以及沟槽深度水平的数目。但是其他实施例可以使用不同 数目的沟槽深度水平和闪耀周期。

通过索引的空间排列来执行1D闪耀锯齿阵列的随机化，也就是通过各 个像素的像素位置的随机改变，正如图3B中所示出的那样。这样就提供了RIP漫射器的“随 机化锯齿结构”，正如图3C中所示出的那样。图3D示出了RIP漫射器的典型灰度表示。通过使 用来自阵列中的可能索引集合的随机(作为示例是均匀)分布来选择像素位置索引的新的 随机数值，作为示例这是通过求诸于MATLAB中的排列函数“randperm”。“randperm”被如下 定义：randperm(n)返回整数1:n的随机排列(利用1:n范围中的整数的均匀分布)。

HW随机化器设计

HW随机化器具有像素化结构，其像素尺寸和设置与图像传感器精确地匹配。因此， 图像传感器上的每一个像素接收由漫射器产生的随机化图像的其强度部分乘以一个随机 系数。HW随机化器可以由在感兴趣的波长范围内是透明的任何材料制成，例如玻璃或塑料。 每一个HW随机化器像素具有处于完全透明与完全不透明之间的随机透明度数值。HW随机化 器被放置在RIP漫射器与图像传感器之间的光径中，并且优选地邻近(或者替换地是精确地 处在)图像传感器平面处。随机化器打破测量矩阵的块Toeplitz结构，从而产生对应于信号 的随机结构。其随机地改变由图像传感器的每一个像素接收到的幅度，从而改进满足RIP条 件的能力。

在一些实施例中，随机化器设计使用来自函数的伪随机数产生与图像传感器像素 矩阵具有相同尺寸的矩阵。对应于随机化器矩阵的各个元素的数值是随机的，其优选地由 具有标准正态分布的独立高斯随机变量给出，并且其概率密度具有以下形式：

$P_{normal}\left(x\right)=\frac{1}{\sqrt{2\pi }}\exp (-\frac{x^2}{2})---\left(46\right)$

应当提到的是，对于该目的也可以使用其他概率密度。在其他实施例中，对应于随 机化器矩阵的各个元素的数值可以是不相关的随机变量，或者是通过回归等式描述的1D或 2D随机过程或场的像素，这对于本领域技术人员是众所周知的。在其他实施例中，对应于随 机化器矩阵的各个元素的数值可以是确定性的，并且由闭型(closedform)等式定义，例如 在各个单独的小透镜之间具有相等或可变的相位偏移的小透镜阵列。

图4示出了具有16x16像素的示例性随机透明度随机化器设计。这种设计当然可以 被扩展到任何尺寸，以便匹配所有图像传感器像素。所述随机化器与波长无关，这是在于其 仅改变到达传感器处的光的强度而不改变其相位。因此其对于每一个波长具有相同的效 果，并且按照相同的方式改变测量矩阵的每一个块，正如前面所描述的那样。可以通过产生 具有对应于传感器上的像素数目的一定数目的元件的随机矩阵的模拟照片(analogue photograph)相当容易地制作这样的随机化器。其结果可以是胶片上的具有正确的尺寸(为 了覆盖图像传感器)和正确的分辨率的图像。

SW随机化器

对应于SW随机化器矩阵的各个元素的数值可以与前面对于HW随机化器所描述的 相同。此外，具有SW随机化器的设备和方法实施例可以使用单项随机化或多项随机化。后者 使用SW随机化器的多个随机版本，其中在检测器上的相同位置处，一个版本的零像素后面 跟着另一个版本的非零像素。这样就导致随机化动作的统计求平均，并且允许高效地使用 由检测器所采集的全部光通量。在一个实施例中，随后对于SW随机化器的每一个版本单独 运行数据立方体重建算法，从而提供重建的数据立方体的几个版本。随后可以通过简单的 每像素求平均或者通过图像融合算法将这些版本合并。在另一个实施例中，可以利用一种 算法来实施数据立方体重建，所述算法采用仅在SW随机化器版本方面存在不同的一类测量 矩阵，并且同时依赖于相同的RIP漫射器。在该实施例中，基于所有检测器像素，具有随机化 器的多个随机版本的单个迭代处理将直接提供重建的数据立方体。

“多项随机化”处理还可以被如下描述：通过求诸于从统计总体取得的随机化矩阵 R的几个独立版本可以提供更多随机化。随后可以对于每一个 单独运行CS迭代。所述迭代得到图像传感器平面处的通过软件根据 等式(16)从相同的DD图像等式(15)获得的N_R幅不同的随机化图像的集合。随后可以通过软 件对几个(即“N_R”个)重建的数据立方体求平均或者将其融合在一起，以便具有对应于数据 立方体的单个解。虽然这样可能需要更多的处理资源，但是通过使用由摄影机采集的全部 光通量可以得到SI摄影机的更高的有效灵敏度。

用于光谱图像重建的基于样条的标架

在这里我们将详细描述基于CS的光谱成像中的用于光谱图像重建的基于样条的 标架(或小标架)。之前已经报告过基于样条的标架，并且特别是从内插和准内插多项式样 条开发多种低通滤波器h⁰，例如参见A.Averbuch和V.Zheludev的“Constructionofbi- orthogonaldiscretewavelettransformusinginterpolatorysplines(利用内插样 条构造双正交离散小波变换)”(AppliedandComputationalHarmonicAnalysis，12，15- 56，2002年)，A.Averbuch、A.B.Pevnyi和V.A.Zheludev的“Bi-orthogonalButterworth waveletsderivedfromdiscreteinterpolatorysplines(从离散内插样条导出的双正 交Butterworth小波)”(IEEETrans.SignalProcessing，49(11)，2682-2692，2001年)，以 及APZFrame。在Bregman迭代的过程中，对于从随机化输入导出的相继的近似X^k应用所述基 于样条的小标架变换。

来自Π[N](其是N-周期性信号的空间)的信号的系统(L>N)形成空 间Π[N]的标架，前提是存在正常数A和B从而使得对于任何信号x＝{x[k]}∈Π[N]，下面的 不等式都成立：如果所述标架限制A＝B，则说所述 标架是紧密的。如果是一个标架，则存在另一个标架(合成)，从而使得任 何信号x＝{x[k]}∈Π[N]都由表示。如果(其也被称作“分析”标 架)是紧密的，则合成标架可以是

具有下采样因数2的分析四通道滤波器组和合成滤波器组 形成完美重建(PR)滤波器组，前提是任何信号x＝{x[k]}∈Π[N]可以被如 下展开：

$x\left[l\right]=\underset{s=0}{\overset{3}{\Sigma }}\underset{k=0}{\overset{N/2-1}{\Sigma }}{h}^s[l-2k]d^s\left[k\right],l=0,...,N-1,---\left(47\right)$

等式(47)提供了信号x的标架展开，其中信号$\left\{{\tilde{h}}^s\right[·-2k\left]\right\},s=0,...,3,k=0,\mathrm{...},N/2-1$构成分析标架，而信号$\left\{h^s\right[·-2k\left]\right\},s=0,...,3,k=0,...,N/2-1$则形成合成标架。

分别通过x₀＝{x[2k]}∈Π[N/2]以及通过x₁＝{x[2k+1]}来标示信号x∈Π[N]的 偶数和奇数多相位分量。标示ω＝e^2πi/N。$\hat{x}\left[R\right]={\Sigma }_{k=0}^{N-1}{\omega }^{-kn}x\left[k\right]$和 ${\hat{x}}_p\left[n\right]={\Sigma }_{k=0}^{N/2-1}{\omega }^{-2kn}{x}_p\left[k\right](p=0,1)$是信号x的离散傅立叶变换(DFT)及其多相位 分量。和${\tilde{h}}_p^s(p=0,1,s=0,...,3)$是滤波器h^s和的多相位分量。和 是其DFT。作出如下标示：

$\tilde{P}\left[n\right]={\left(\begin{array}{c}{\hat{\tilde{h}}}_0^0\left[n\right]\mathrm{...}{\hat{\tilde{h}}}_0^3\left[n\right]\\ {\hat{\tilde{h}}}_1^0\left[n\right]\mathrm{...}{\hat{\tilde{h}}}_1^3\left[n\right]\end{array}\right)}^T$以及$P\left[n\right]=\left(\begin{array}{c}{\hat{h}}_0^0\left[n\right]...{\hat{h}}_0^3\left[n\right]\\ {\hat{h}}_1^0\left[n\right]...{\hat{h}}_1^3\left[n\right]\end{array}\right).---\left(48\right)$

其中，和P[n]分别是滤波器组和H的分析和合成多相位矩阵。符号(…)^T意 味着矩阵转置。长度为N的信号x的直接小标架变换产生四个系数集合d＝{d^s}(s＝0,1,2, 3)，其中每一个系数集合包含N/2个成员，所述直接小标架变换可以被表示成：

$\left(\begin{array}{c}{d}^0\left[n\right]\\ .\\ .\\ .\\ d^3\left[n\right]\end{array}\right)=\tilde{P}[-n]\left(\begin{array}{c}{\hat{x}}_0\left[n\right]\\ {\hat{x}}_1\left[n\right]\end{array}\right)---\left(49\right)$

从系数d^s(s＝0,1,2,3)还原信号的逆小标架变换是：

$\left(\begin{array}{c}{\hat{x}}_0\left[n\right]\\ {\hat{x}}_1\left[n\right]\end{array}\right)=P\left[n\right]\left(\begin{array}{c}{d}^0\left[n\right]\\ .\\ .\\ .\\ d^3\left[n\right]\end{array}\right)---\left(50\right)$

因此，长度为N的信号x变成由来自集合d^s(s＝0,1,2,3)的2N个系数表示。在这个 意义上，该表示是双重冗余的。关系式是使得滤波器组对形成PR滤波器组的条件，其中I₂是2x2单位矩阵。来自PR滤波器组的滤波器和h⁰是低通的。

为了把小标架变换扩展到更低分辨率尺度并且为了提高表示冗余度，利用分析多 相位矩阵将所述变换应用于低频系数阵列d⁰。利用合成多相位矩阵和P[2n] 来还原系数阵列d⁰。类似地，利用矩阵和把所述变换扩展 到另外的分辨率尺度。因此，2D阵列的2D小标架变换包括对所述阵列的各列应用1D变换，随 后对所述阵列的各行应用1D变换。

在一个示例性实施例中，我们设计了具有不同系数的一系列4通道PR滤波器组(参 见APZFrame)。其多相位矩阵具有特定结构，所述结构由频率响应为 的低通滤波器决定：

$\tilde{P}\left[n\right]=\left(\begin{array}{cc}{\hat{h}}_0^0\left[n\right]& {\hat{h}}_1^0\left[n\right]\\ -{\hat{h}}_1^0[-n]& {\hat{h}}_0^1\left[n\right]\\ {\tilde{T}}^2\left[n\right]& 0\\ 0& {\tilde{T}}^3\left[n\right]\end{array}\right),P\left[n\right]=\left(\begin{array}{cccc}{\hat{h}}_0^0\left[n\right]& -{\hat{h}}_1^0[-n]& T^2\left[n\right]& 0\\ {\hat{h}}_1^0\left[n\right]& {\hat{h}}_0^0[-n]& 0& T^3\left[n\right]\end{array}\right).---\left(51\right)$

其中，$T^2\left[n\right]{\tilde{T}}^2[-n]=T^3\left[n\right]{\tilde{T}}^3[-n]=1-\left|{\hat{h}}_0^0\right[n]{|}^2+|{\hat{h}}_1^0\left[n\right]{|}^2.$如果 $T^2\left[n\right]={\tilde{T}}^2[-n]=T^3\left[n\right]={\tilde{T}}^3[-n],$则滤波器组生成紧密标架。如果 则滤波器组生成半紧密标架。不同于紧密标架 滤波器组，生成半紧密标架的滤波器组具有线性相位。

在一个实施例中，我们使用从准内插二次样条导出的滤波器组(参见APZFrame)。 该滤波器组生成半紧密标架。分析和合成滤波器的频率响应为如下：

$\begin{array}{cc}{\hat{h}}^0\left[n\right]{\hat{\tilde{h}}}^0\left[n\right]=\frac{1}{\sqrt{2}}{\cos }^4\frac{\pi n}{N}(3-\cos \frac{2\pi n}{N}),& {\hat{h}}^1\left[n\right]{\hat{\tilde{h}}}^1\left[n\right]=\frac{{\omega }^{-n}}{\sqrt{2}}{\sin }^4\frac{\pi n}{N}(3+\cos \frac{2\pi n}{N}),\\ {\hat{h}}^2\left[n\right]=T\left[n\right],{\hat{h}}^3\left[n\right]=-{\omega }^{-n}\tilde{T}[-n],& {\hat{\tilde{h}}}^2\left[n\right]=\tilde{T}\left[n\right],{\hat{\tilde{h}}}^3\left[n\right]=-{\omega }^{-n}T[-n],\end{array}---\left(52\right)$

其中，序列T[n]、和G[n]是：

$\begin{array}{c}\tilde{T}\left[n\right]=\sqrt{2}\frac{{\omega }^{4n}-3{\omega }^{2n}+3-{\omega }^{-2n}}{8}T\left[n\right]=\sqrt{2}\frac{(1-{\omega }^{2n})G\left[n\right]}{1024},\\ G\left[n\right]=-{\omega }^{4n}-12{\omega }^{2n}+346-12{\omega }^{-2n}-{\omega }^{-4n}.\end{array}---\left(53\right)$

图5对于使用在SI的一个实施例中的特定标架显示出滤波器和 的脉冲和量值响应的曲线图。图5A示出了滤波器h⁰(实线)和h¹(虚线)的脉冲 响应。图5B示出了滤波器h⁰(实线)和h¹(虚线)的量值响应。图5C示出了滤波器h²和(实 线)以及h³和(虚线)的脉冲响应。图5D示出了滤波器h²和(实线)以及h³和(虚线) 的量值响应。

在Bregman迭代的过程中，重复等式(33)和(34)，即直接和逆2D小标架变换。每一 项2D小标架变换通过利用快速傅立叶变换对矩阵的各列应用1D小标架变换并且随后对各 行应用1D变换来实施，即等式(49)和(50)。在等式(51)中定义的多相位矩阵和P[n]被 用于一层变换，而多相位矩阵和P[2^mn](m＝2,3,…)则被用于多层变换。

使用方法

图6A在流程图中示意性地示出了这里所公开的用于SSI的方法的一个实施例。在 图6A中，在步骤502中，利用例如设备100之类的设备采集使用满足与感知矩阵A有关的RIP 条件的漫射器的对象的快照图像Y。在一些实施例中，所述采集还可以通过附加的随机化器 来进行。对A^TY应用的2D稀疏化变换D在步骤504中被使用来获得数据立方体X的稀疏版本d， 其中A^T是与感知矩阵A有关的转置矩阵。随后在步骤506中使用逆变换Ψ从d重建X。随后在 步骤508中使用CS迭代方案(作为示例是分裂Bregman迭代)来获得L个光谱带当中的每一个 光谱带中的对象强度的图像。

图6B示出了图6A的流程图中的步骤506的细节。在步骤520中，感知矩阵A被转置成 转置矩阵A^T，其被用来构造与逆小标架变换有关的矩阵Ψ。在步骤522中，利用Y、A、A^T、D和Ψ 实施分裂Bregman迭代，以便近似d。在步骤524中，通过Ψ对d进行变换以获得X。

图6C示出了对应于分裂Bregman迭代的图6B中的步骤522的细节。在步骤540中进 行阵列b、c和d的初始化。在步骤542中实施d-b的逆2D小标架变换Ψ(d-b)，在步骤544中求 解线性等式(μA^TA+χ)X＝μA^T(Y-c)+χΨ(d-b)，并且在步骤546中实施X的小标架变换DX。在步 骤548中把DX+b的每一个元素与预定的收缩阈值χ^-1进行比较以便获得对应于矢量d的新的 数值，在步骤550中计算对应于b的新的数值b^new＝b+(DX-b)，并且在步骤552中计算对应于c 的新的数值c^new＝c+(AX-Y)。只要误差范数||AX-Y||超出预定的停止阈值σ，就从步骤542重 复迭代。

数据立方体重建的仿真

针对利用装配有1DRIP漫射器以及不具有或具有随机化器的数字摄影机感知的 测试多光谱来源对象运行了数据立方体重建的各种计算机仿真。所述仿真是利用Matlab软 件代码进行的。

在包括数字摄影机以及不具有随机化器的1DRIP漫射器的设备中进行的2DCC-SCR多光谱图像的仿真

在前面关于2DCS-SCR的描述的基础上实施利用MATLAB进行的仿真。光谱数据来 源是D.H.Foster的“Hyperspectralimagesofnaturalscenes2004(自然场景的超光谱 图像2004)”中的第7号片段“房屋”场http://personalpages.manchester.ac.uk/staff/david.foster/Hyperspectral_images_of_natural_scenes_04.html，2004年(后文中称 作“Porto的房屋”)。通过包括数字摄影机的光学成像系统的计算机仿真获得了DD图像，所 述数字摄影机具有成像透镜和像素化图像传感器，其中在图像传感器平面处的光瞳中插入 了1DRIP漫射器。DD图像的每一列是相应的来源对象图像列中的所有光谱和空间数据与感 知矩阵的线性组合。应用Bregman迭代处理，以便重建对应于光谱立方体的N列和L＝33个光 谱带当中的每一个中的M个体素的光谱立方体信息。其结果是一个矢量集合，每一个矢量包 括对应于相应的图像列中的每一个像素的所有光谱信息。随后把所有重建的图像列彼此并 排放置，从而提供代表完全光谱立方体的全部光谱信息。最后，通过取得对应于所需的光谱 带的M行的接连集合，对光谱立方体进行处理以便获得对象的L幅单独的光谱图像。通过把 我们的结果的PSNR与在所报告的研究中获得的PSNR进行比较，对2DCS-SCR结果的质量进 行评估。在仿真中执行不具有或具有随机化器的选项。

表1总结了在仿真中使用的光学系统和所设计的RIP漫射器的参数。所述参数适合 于10Mp摄影机。

表1

表2对于具有列尺寸M_image、L光谱带和20％稀疏度(“稀疏”图像中的非零数值的部 分)的图像提供了遵循等式(32)满足RIP条件等式(31)所需要的传感器的图像上的最小行 数M。

N L LN K＝0.2LN N_y128 5 640 128 183 128 9 1152 230 330 128 24 3072 614 879 128 33 4224 844 1209 256 33 8448 1690 2417

表2

图7A-J示出了对应于从利用装备有1DRIP漫射器(图3)的数字摄影机获得的快照 图像进行的2DCS-SCR的仿真结果，即在“不具有随机化器”的一个实施例中。光谱数据立方 体的来源是具有L＝33个波长带的Porto场景中的多光谱房屋的一个片段(灰度版本处在图 7I中，左侧)。在图7A-F中的左侧示出了从来源光谱立方体提取出的33幅单色图像当中的6 幅。其光谱带是编号为l＝12、17、19、26、28和33的波长。图7G示出了感知矩阵A的图形表示， 其由匹配所述波长的33个块构成。在图7H中示出了在图像传感器上接收到的对应的DD图 像。图7A-F在右侧示出了6个采样波段处的重建图像，并且作为比较在左侧示出了对应的原 始图像。计算机仿真结果提供了以下PSNR比值：图7A，l＝12，PSNR＝26.08dB；图7B，l＝17， PSNR＝27.81dB；图7C，l＝19，所计算的PSNR＝29.03dB；图7D，l＝26，所计算的PSNR＝ 28.87dB；图7E，l＝28，所计算的PSNR＝30.3dB；图7F，l＝33，PSNR＝26.46dB。

图7I示出了从33幅原始来源图像编辑的灰度图像(左侧)相比于从33幅所还原的 波段图像编辑的图像(右侧)。在左侧的椭圆形中示出了四个采样点。在图7J中示出了各个 采样点处的完全光谱。具体来说，所述四个采样点处的对应于33个光谱带的原始(虚线)与 重建(实线)光谱之间的比较表明了高质量重建。

图8A-8J示出了对应于从利用装备有1DRIP漫射器(图3)和SW随机化器的数字摄 影机获得的快照图像进行的2DCC-SCR的仿真结果，即在“具有随机化器”的一个实施例中。 光谱数据立方体来源与图7A-7J中相同。在图8A-F中的左侧示出了从整个来源光谱立方体 提取出的33幅单色图像当中的6幅。其光谱带是编号为l＝12、17、19、26、28和33的波长。在 图4中已经示出了随机化器的一个片段。在仿真中，矩阵A被修改以便考虑到随机化器，正如 前面所描述的那样。图8G示出了随机化感知矩阵A的图形表示，其由匹配所述波长的33个块 构成。自此之后，所有计算都按照类似于图7A-7J的“不具有随机化器”实施例的方式进行。 RIP漫射器和随机化器的组合使用在摄影机传感器处提供了随机化图像。这样的图像在图 8H中示出。图8A-F在右侧示出了6个采样波段处的重建图像，并且作为比较在左侧示出了对 应的原始图像。不具有随机化器的计算机仿真的结果提供了以下所计算的PSNR比值：图8A， l＝12，PSNR＝25.5；图8B，l＝17，PSNR＝26.71；图8C，l＝19，PSNR＝27.44；图8D，l＝26， PSNR＝27.31；图8E，l＝28，PSNR＝28.74；图8F，l＝33，PSNR＝25.8。

图8I示出了从33幅原始来源图像编辑的灰度图像(左侧)相比于在具有随机化器 的情况下从33幅所还原的波段图像编辑的图像(右侧)。在左侧的椭圆形中示出了四个采样 点。在图8J中示出了各个采样点处的完全光谱。具体来说，所述四个采样点处的对应于33个 光谱带的原始(虚线)与重建(实线)光谱之间的比较表明了在具有随机化器的情况下的高 质量重建。尽管随机化器可能会略微降低计算机仿真中的PSNR，但是其对于感知矩阵提供 了附加的随机性。所述附加的随机性可以对于在某些光学系统中满足RIP条件以及改进SCR 有所贡献。

附加的设备实施例

图9A示意性地示出了这里所公开的SSI设备的一个实施例900。设备900包括：数字 摄影机902，RIP漫射器904，图像传感器906，成像透镜908，系统孔径910，以及随机化器918。 可选的是，其可以包括1D(或2D)发散器912和带通光谱滤波器914。图像传感器906提供DD图 像916。RIP漫射器(并且可选地还有发散器)可以被放置在系统孔径光圈的平面处，从而允 许漫射器实质上在平面光波模式下运作，从而最小化固有像差并且改进图像质量。在一个 实施例中，RIP漫射器可以是透镜块的一部分，并且被放置在摄影机内部。在另一个实施例 中，RIP漫射器可以被放置在数字摄影机外部，并且处于摄影机透镜与对象之间。在另一个 实施例中，RIP漫射器可以被放置在数字摄影机透镜内部，并且处于透镜与图像传感器之 间。随机化器918可以是HW或SW随机化器。组件透镜的设计可以与漫射器(如果存在的话还 有发散器)的光发散属性相匹配。对于DD图像和其他信息的处理由数字处理器920实施，其 可以被集成在数字摄影机中或者可以处于摄影机外部。这样的处理器存在于所有设备实施 例图9B-13中，但是为了简单起见没有被示出。由于小型化的数字摄影机是已知的，并且由 于RIP漫射器可以被制作成纤薄漫射-衍射光学元件的形式，因此例如设备900之类的设备 显然可以被小型化，以便例如用在智能电话或诊断药丸(diagnosticpill)中。

图9B示意性地示出了这里所公开的SSI设备的另一个实施例900’。类似于设备 900，设备900’包括数字摄影机902’、RIP漫射器90以及随机化器918，其不同之处在于摄影 机902’是包括用于成像透镜908的两个孔径而不是单个系统孔径910的双孔径摄影机。因 此，图像传感器906被双图像传感器906’替代，其捕获常规图像916’和DD图像916全部二者。 在其他实施例中，可以有多于两个孔径。一般来说，这里所公开的设备还可以使用具有导向 单个图像传感器的两条光径的设计，例如参见美国专利申请号2007/0285554。一条路径可 以经过RIP漫射器以产生DD图像，另一个则直接到达图像传感器以提供常规图像。

图10-13示出了这里所公开包括两条成像通道的SSI设备的其他实施例：一条具有 RIP漫射器，另一条具有如数字摄影机中的常规成像。所有后面的实施例都提供DD图像和常 规图像全部二者。DD图像可以通过随机化器918的动作被转换成随机化图像。在必要时将两 幅图像相减，以用于在不相干光学成像系统中实施负值感知矩阵。

图10示意性地示出了这里所公开的SSI设备的一个实施例1000。该实施例包括反 射-折射光束分离器和两个单独的图像传感器。除了前面所提到的元件之外，设备1000还包 括成像透镜前块1004、光束分离器1006、系统孔径光圈1010、第一成像透镜后块1012、第二 成像透镜后块1019以及用于提供常规图像1022的常规图像传感器1020。可选的是，设备 1000可以包括1D(或2D)发散器1024。分离器1006可以是半透明反射镜或光束分离立方体， 其指向与光轴成大于45度角，以便将光反射到透镜1019。RIP漫射器(并且可选地还有发散 器)优选地被放置在系统孔径1010的平面处。

图11示意性地示出了这里所公开的SSI设备的另一个实施例1100。与设备1000一 样，设备1100包括安装在仅捕获DD图像916的双图像传感器1004的一部分前方的随机化器。 与设备1000不同，设备1100包括附加的转动反射镜1102，捕获发散-漫射图像和常规图像全 部二者的双图像传感器1104，以及第二系统孔径1106。在这里，分离器1006的指向与光轴成 接近45度角，以便把光反射到转动反射镜1102。所述双图像传感器是通过把发散-漫射图像 传感器与常规图像传感器对准而获得的。

图12示意性地示出了这里所公开的SSI设备的另一个实施例1200。与设备1000和 1100不同，设备1200包括衍射发散器1202，其同时充当光束分离器和1D发散器。此外，与图 11中一样还包括双图像传感器。衍射发散器1202可以被实施成衍射光栅，其把光分离到第0 和第1衍射阶次1204和1206，其中0阶没有发散并且形成常规图像1208，1阶具有发散并且在 双图像传感器上形成DD图像1210。随机化器优选地被安装在仅捕获DD图像的双图像传感器 的一部分的前方。

图13示意性地示出了这里所公开的SSI设备的另一个实施例1300。设备1300包括 如图12中的衍射发散器，但是具有单个后块1302而不是两个单独的成像透镜的后块。前块 1304和系统孔径光圈与后块1302被包装在一起，从而形成成像透镜。在一些实施例中，成像 透镜的前块或者(多个)后块可以缺失，从而使得成像透镜分别仅由(多个)后块或者仅由前 块形成。RIP漫射器以及(如果存在的话)发散器的空间位置优选地更靠近系统孔径(或者成 像透镜的入射或出射光瞳)而不是图像传感器。随机化器优选地被安装在仅捕获DD图像的 双图像传感器的一部分的前方。

出于这里所阐述的所有目的，在本申请中提到的每一份出版物都被全文合并在此 以作参考。应当强调的是，对于本申请中的任何参考文献的引用或标识不应当被理解成承 认这样的参考文献可用作或者被承认作为现有技术。虽然本公开内容描述了有限数目的实 施例，但是应当认识到，这样的实施例可以有许多变型、修改和其他应用。举例来说，虽然说 明书特别涉及通过特定RIP漫射器获得的DD图像，但是这里所描述的2DCS-SCR方法也可以 被应用于为之定义了感知矩阵A的任何其他发散-漫射图像。举例来说，虽然说明书特别涉 及小标架变换和Bregman迭代，但是其他类型的变换或算法也可以被用于这里所描述的 SCR。此外，虽然详细描述了Toeplitz矩阵和卷积，但是也可以使用对应于非近轴和空间变 化的光学系统和/或具有像差的光学系统的更加一般的矩阵和线性变化。举例来说，虽然详 细描述了1D漫射/发散，但是也可以使用2D漫射/发散。因此，本公开内容应当被理解成不限 于小标架变换、分裂或其他Bregman迭代以及1D漫射/发散。一般来说，本公开内容应当被理 解成不受限于这里所描述的具体实施例，而是仅由所附权利要求书的范围限制。