傅里叶单像素成像
文章目录
张子邦, 陆天傲, 彭军政,等.傅里叶单像素成像技术与应用[J]. 红外与激光工程, 2019, 48(6).、
数学原理
通过获取图像的傅里叶谱来重建图像本身。图像的傅里叶谱通过对目标物体投影一些列不同的傅里叶基地图案，并通过单像素探测器测量所得的光强值来获得。图像的重建过程仅仅是依次二维离散逆傅里叶变换。利用自然图像在傅里叶域的稀疏性以实现高效率的成像。
对一个目标物体进行二维成像，本质上就是一个对该目标物体的光强分布的空间信息获取和再现过程。二维单像素成像，就是通过仅适用一个没有空间分辨能力的探测器——单像素探测器——来获取物体的空间信息，并对所获取的信息施以响应的算法来重建物体图像。由此可见，二维单像素成像本身包含了一对关于空间分辨能力的矛盾，即：成像过程对空间分辨能力的需求而单像素探测器缺乏空间分辨能力的矛盾。
基于傅里叶分析，1807年，J.Fourier提出傅里叶分析理论，指出任意一个来纳许的信号都可以通过傅里叶基得线性组合构成(傅里叶正变换的意义)。
换言之，任意一个连续信号都可以适用无穷多个空间频率不同的正弦或余弦波进行线性叠加(傅里叶逆变换的物理意义)
二维傅里叶正变换和逆变换：
C
(
f
x
,
f
y
)
=
∫
−
∞
+
∞
∫
−
∞
+
∞
I
(
x
,
y
)
e
x
p
[
−
j
⋅
2
π
(
f
x
x
+
f
y
y
)
]
d
f
x
d
f
y
C(f_x,f_y)=\int_{-\infty}^{+\infty}\int_{-\infty}^{+\infty}I(x,y)exp\left[-j\cdot2\pi(f_xx+f_yy)\right]df_xdf_y
C(fx,fy)=∫−∞+∞∫−∞+∞I(x,y)exp[−j⋅2π(fxx+fyy)]dfxdfy
I
(
x
,
y
)
=
∫
−
∞
+
∞
∫
−
∞
+
∞
C
(
f
x
,
f
y
)
e
x
p
[
−
j
⋅
2
π
(
f
x
x
+
f
y
y
)
]
d
x
d
y
I(x,y)=\int_{-\infty}^{+\infty}\int_{-\infty}^{+\infty}C(f_x,f_y)exp\left[-j\cdot2\pi(f_xx+f_yy)\right]dxdy
I(x,y)=∫−∞+∞∫−∞+∞C(fx,fy)exp[−j⋅2π(fxx+fyy)]dxdy
- x，y是空间域中的直角坐标；
- f x , f y f_x,\;\;f_y fx,fy是傅里叶域中的直角坐标，对应x,y方向的空间频率；
- I ( x , y ) I(x,y) I(x,y)表示二维图像
- C ( f x , f y ) C(f_x,f_y) C(fx,fy)表示二维图像的傅里叶谱
- j是虚数单位
扩展到离散信号：
C
(
u
,
v
)
=
∑
x
=
0
M
−
1
∑
y
=
0
N
−
1
I
(
x
,
y
)
e
x
p
[
−
j
2
π
(
u
x
M
+
v
y
N
)
]
C(u,v)=\sum_{x=0}^{M-1}\sum_{y=0}^{N-1}I(x,y)exp\left[-j2\pi\left(\frac{ux}{M}+\frac{vy}{N}\right)\right]
C(u,v)=x=0∑M−1y=0∑N−1I(x,y)exp[−j2π(Mux+Nvy)]
I
(
x
,
y
)
=
∑
x
=
0
M
−
1
∑
y
=
0
N
−
1
C
(
u
,
v
)
e
x
p
[
−
j
2
π
(
u
x
M
+
v
y
N
)
]
I(x,y)=\sum_{x=0}^{M-1}\sum_{y=0}^{N-1}C(u,v)exp\left[-j2\pi\left(\frac{ux}{M}+\frac{vy}{N}\right)\right]
I(x,y)=x=0∑M−1y=0∑N−1C(u,v)exp[−j2π(Mux+Nvy)]
任何一个二维图像均为一系列傅里叶基地图案(不同空间频率和初相为的余弦条纹图案)的线性叠加的结果(即线性组合)。额日对应个傅里叶基地图案的权重，就是傅里叶系数 C ( u , v ) C(u,v) C(u,v)
通过单像素探测器获取二维图像的傅里叶谱
图像通过傅里叶正变换，可以得到图像在傅里叶域的图像——傅里叶谱。
获取图像的傅里叶谱的过程就是获取图像对应每一个傅里叶基地图案权重的过程。为获得物体图像的傅里叶系数，首先适用计算机产生傅里叶基地图案。并将所得的图案投影到目标物体上，在适用单像素探测器对所得的光信号强度进行探测。
傅里叶基地图案本质上就是一系列具有空间频率、初相位不同的强度成余弦分布的图像。适用数学表达式可表示为：
P
(
x
,
y
;
f
x
,
f
y
,
ϕ
)
=
a
+
b
⋅
c
o
s
(
2
π
f
x
x
+
2
π
f
y
y
+
ϕ
)
P(x,y;f_x,f_y,\phi)=a+b\cdot cos(2\pi f_xx+2\pi f_yy+\phi)
P(x,y;fx,fy,ϕ)=a+b⋅cos(2πfxx+2πfyy+ϕ)
- a是平均光强、
- b是对比度
- x，y是目标物体所在平面的直角坐标，
- f x f_x fx和 f y f_y fy分别是对应x、y方向的空间频率
- ϕ \phi ϕ是初相位
设物体在单像素探测器方向测量投影仪投影傅里叶基地图案时所采用的光照方向的反射强度为
R
(
x
,
y
)
R(x,y)
R(x,y)，记物体图像为
I
(
x
,
y
)
I(x,y)
I(x,y)，则物体图像和物体的反射强度分布关系为
I
(
x
.
y
)
∝
R
(
x
,
y
)
I(x.y)\propto R(x,y)
I(x.y)∝R(x,y)。目标物体在傅里叶基底图案
P
(
x
,
y
;
f
x
,
f
y
,
ϕ
)
P(x,y;f_x,f_y,\phi)
P(x,y;fx,fy,ϕ)照明下，所得的反射光的光强为：
E
ϕ
(
f
x
,
f
y
)
=
∬
S
R
(
x
,
y
)
⋅
P
(
x
,
y
;
f
x
,
f
y
,
ϕ
)
d
x
d
y
E_{\phi}(f_x,f_y)=\iint_SR(x,y)\cdot P(x,y;f_x,f_y,\phi)dxdy
Eϕ(fx,fy)=∬SR(x,y)⋅P(x,y;fx,fy,ϕ)dxdy
其中，
S
S
S为傅里叶基地图案的投影区域。单像素探测器的光响应值为：
D
ϕ
(
f
x
,
f
y
)
=
D
n
+
β
⋅
E
ϕ
(
f
X
,
f
y
)
D_\phi(f_x,f_y)=D_n+\beta\cdot E_\phi(f_X,f_y)
Dϕ(fx,fy)=Dn+β⋅Eϕ(fX,fy)
- D n D_n Dn是背景照明在探测器位置引起的光响应值，
- β \beta β是一个与单像素探测器放大倍数、探测器与物体的空间关系有关的因子。
为获取对应空间频率为
(
f
x
,
f
y
)
(f_x,f_y)
(fx,fy)的傅里叶系数
C
(
f
x
,
f
y
)
C(f_x,f_y)
C(fx,fy)，j将空间频率为
(
f
x
,
f
y
)
(f_x,f_y)
(fx,fy)而初位相依次位
0
,
π
/
2
,
π
,
3
π
/
2
0,\pi/2,\pi,3\pi/2
0,π/2,π,3π/2的四个傅里叶基底图案，分别投影到物体上。这四个傅里叶基图案分别记为
P
1
,
P
2
,
P
3
P_1,P_2,P_3
P1,P2,P3和
P
4
P_4
P4:
那么单像素探测器依次采集信号可以得到响应值
D
1
,
D
2
,
D
3
D_1,D_2,D_3
D1,D2,D3和
D
4
D_4
D4：
然后根据四步相移法，得到关于反射强度
R
(
x
,
y
)
R(x,y)
R(x,y)的傅里叶正变换积分式：
- j是叙述单位， F F{} F表示傅里叶正变换。
又由于
I
(
x
,
y
)
∝
R
(
x
,
y
)
I(x,y)\propto R(x,y)
I(x,y)∝R(x,y)所以可以认为：
通过上述过程我们就可以得到原图像的每一个傅里叶因子
C
(
f
x
,
f
y
)
C(f_x,f_y)
C(fx,fy)。通过对目标物体投影不同的空间频率的四步相移傅里叶基底图案，可获得物体图像傅里叶谱中对应不同空间频率的傅里叶系数。
对所获得的物体图像傅里叶谱实行逆傅里叶变化，即可实现物体图像的重建：
[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-cGXEb4eb-1621759377110)(leanote://file/getImage?fileId=60a4ae27ab64415fa2000817)]
其中，
F
−
1
{
}
F^{-1}\{\}
F−1{}表示逆傅立叶变换。
要无失真的重建图像分辨率为
M
×
N
M\times N
M×N像素的图像，则需要获取其傅立叶谱中的全部
M
×
N
M\times N
M×N个傅立叶系数。由于每个傅立叶系数需要投影四部相移的傅立叶基底图案，即需要花费4次测量；又由于物体图像再数学上是一个实值得矩阵，其傅立叶谱具有共轭对称性，因此获取
M
×
N
M\times N
M×N个傅立叶系数共需要
2
×
M
×
N
2\times M\times N
2×M×N次测量。换言之，在无失真重建图像的情况下，使用四部相移法获得图像的傅立叶谱，花费的测量数是重建图像像素数的两倍，其根本原因是四步相移方法本质上是一种差分测量。
图像能量主要集中在傅立叶谱的低频部分，而高频的傅立叶系数的模很小甚至接近于零。因此，可以通过只投影低空间频率的傅立叶基底图案，获取低频的傅立叶系数并直接使高频的系数置零来重建物体图像，从而实现较少测量次数的目的。
程序
Single-pixel imaging by means of Fourier spectrum acquisition
- 程序步骤：
- 四步相移正弦条纹图案生成
- P Φ ( x , y , f x , f y ) = a + b ⋅ cos ( 2 π f x x + 2 π f y y + Φ ) P_{\Phi}(x,y,f_x,f_y)=a+b\cdot \cos(2\pi f_x x+2\pi f_y y+\Phi) PΦ(x,y,fx,fy)=a+b⋅cos(2πfxx+2πfyy+Φ)
- 图案投影
- E Φ ( f x , f y ) = ∬ ϕ R ( x , y ) P Φ ( x , y , f x , f y ) d x d y E_{\Phi}(f_x,f_y)=\iint\limits_\phi R(x,y)P_{\Phi}(x,y,f_x,f_y)dxdy EΦ(fx,fy)=ϕ∬R(x,y)PΦ(x,y,fx,fy)dxdy
- 单像素信号获取
- D Φ ( f x , f y ) = D n + k E Φ ( f x , f y ) D_\Phi(f_x,f_y)=D_n+kE_\Phi(f_x,f_y) DΦ(fx,fy)=Dn+kEΦ(fx,fy)
- 得到表面反射率
- [ D 0 ( f x , f y ) − D π ( f x , f y ) ] + j ⋅ [ D π / 2 ( f x , f y ) − D 3 π / 2 ( f x , f y ) ] = 2 b k ⋅ F { R ( x , y ) } \left[D_0(f_x,f_y)-D_\pi(f_x,f_y)\right]+j\cdot \left[D_{\pi/2}(f_x,f_y)-D_{3\pi/2}(f_x,f_y)\right]=2bk\cdot F\{R(x,y)\} [D0(fx,fy)−Dπ(fx,fy)]+j⋅[Dπ/2(fx,fy)−D3π/2(fx,fy)]=2bk⋅F{R(x,y)}
- 2 b k ⋅ R ( x , y ) = F − 1 { [ D 0 ( f x , f y ) − D π ( f x , f y ) ] + j ⋅ [ D π / 2 ( f x , f y ) − D 3 π / 2 ( f x , f y ) ] } 2bk\cdot R(x,y)=F^{-1}\{\left[D_0(f_x,f_y)-D_\pi(f_x,f_y)\right]+j\cdot\left[D_{\pi/2}(f_x,f_y)-D_{3\pi/2}(f_x,f_y)\right]\} 2bk⋅R(x,y)=F−1{[D0(fx,fy)−Dπ(fx,fy)]+j⋅[Dπ/2(fx,fy)−D3π/2(fx,fy)]}粗体
- 四步相移正弦条纹图案生成
具体步骤
- 读取图片
- 参数设置
- 相移步数
- 相移量
- 采样率
- 采样方式（圆形）
- 频率坐标图
- 频率采样顺序
- 获取傅立叶域对称中心
- 生成投影图案
- 进行数据采集
- 根据公式利用采集的数据计算频率因子(在这个过程中要补全对称的数据)
- 通过逆傅立叶变换恢复图像
** 本文只是学习过程中的总结，如有侵权请联系 **
1636
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
