滤波
统计影像的均值和方差时,常规做法是先循环一遍计算均值,再循环一遍计算方差,但通过推导方差计算公式,可以简化上述操作,缩短统计耗时。
float mean = 0 , std = 0 ; for ( int i = 0 ; i < arr_size ; i ++ ) { mean += arr [ i ] ; std += arr [ i ] * arr [ i ] ; } mean /= valid_num ; std = sqrtf ( std / ( valid_num - 1 ) - valid_num / ( valid_num - 1 ) * mean * mean ) ;
均值、方差的公式
x ˉ = 1 n ∑ x i S t d = 1 n − 1 ∑ ( x i − x ˉ ) 2 = 1 n − 1 ( ∑ x i 2 − n x ˉ 2 ) \begin{aligned}
\bar{x} =& \frac{1}{n}\sum x_i \\
Std =& \sqrt{\frac{1}{n-1}\sum {(x_i - \bar{x})}^2} \\
=& \sqrt{\frac{1}{n-1}(\sum{{x_i}^2}-n{\bar{x}}^2)}
\end{aligned} x ˉ = St d = = n 1 ∑ x i n − 1 1 ∑ ( x i − x ˉ ) 2 n − 1 1 ( ∑ x i 2 − n x ˉ 2 )
方差的推导过程,
S t d 2 = 1 n − 1 ∑ ( x i − x ˉ ) 2 = 1 n − 1 ∑ ( x i 2 − 2 x ˉ ⋅ x i + x ˉ 2 ) = 1 n − 1 ( ∑ x i 2 − 2 x ˉ ⋅ ∑ x i + n x ˉ 2 ) = 1 n − 1 ( ∑ x i 2 − 2 n x ˉ 2 + n x ˉ 2 ) = 1 n − 1 ( ∑ x i 2 − n x ˉ 2 ) ⇒ S t d = 1 n − 1 ( ∑ x i 2 − n x ˉ 2 ) \begin{aligned}
&{Std}^2&=& \frac{1}{n-1}\sum {(x_i - \bar{x})}^2 \\
& &=& \frac{1}{n-1}\sum ({x_i}^2 - 2\bar{x}\cdot x_i+{\bar{x}}^2) \\
& &=& \frac{1}{n-1}(\sum {x_i}^2 - 2 \bar{x} \cdot \sum x_i + n {\bar{x}}^2) \\
& &=& \frac{1}{n-1}(\sum {x_i}^2 - 2 n {\bar{x}}^2 + n {\bar{x}}^2) \\
& &=& \frac{1}{n-1}(\sum {x_i}^2 - n {\bar{x}}^2) \\ \\
\Rarr & Std &=& \sqrt{\frac{1}{n-1}(\sum{{x_i}^2}-n{\bar{x}}^2)}
\end{aligned} ⇒ St d 2 St d = = = = = = n − 1 1 ∑ ( x i − x ˉ ) 2 n − 1 1 ∑ ( x i 2 − 2 x ˉ ⋅ x i + x ˉ 2 ) n − 1 1 ( ∑ x i 2 − 2 x ˉ ⋅ ∑ x i + n x ˉ 2 ) n − 1 1 ( ∑ x i 2 − 2 n x ˉ 2 + n x ˉ 2 ) n − 1 1 ( ∑ x i 2 − n x ˉ 2 ) n − 1 1 ( ∑ x i 2 − n x ˉ 2 )
Qi Zhang, Xiaoyong Shen, Li Xu, and Jiaya Jia, ‘Rolling Guidance Filter’, ECCV 2014, pp. 812–830, 2014.
J 1 ( p ) = 1 K p ∑ q ∈ N ( p ) exp ( − ∥ p − q ∥ 2 2 σ s 2 ) I ( q ) J_1(p)=\frac{1}{K_p}\displaystyle\sum_{q\in N(p)}\exp \bigg(-\frac{ {\|p-q\|}^2} {2\sigma_s^2}\bigg)I(q) J 1 ( p ) = K p 1 q ∈ N ( p ) ∑ exp ( − 2 σ s 2 ∥ p − q ∥ 2 ) I ( q )
K p = ∑ q ∈ N ( p ) exp ( − ∥ p − q ∥ 2 2 σ s 2 ) K_p=\displaystyle\sum_{q\in N(p)}{\exp \big(-\frac{ {\|p-q\|}^2} {2\sigma_s^2}\big)} K p = q ∈ N ( p ) ∑ exp ( − 2 σ s 2 ∥ p − q ∥ 2 )
其中, p p p 是待求点, q q q 是 p p p 周围的所有点的集合, I I I 是原图像, σ s 2 \sigma_s^2 σ s 2 是标准偏差( standard deviation)自己设定。通过这一步可以消除尺度小于 σ s 2 \sigma_s^2 σ s 2 的结构。
(This filter completely removes structures whose scale is smaller than σ s 2 \sigma_s^2 σ s 2 as claimed in the scale space theory. It is implemented efficiently by separating kernels in perpendicular directions. Approximation by box filter is also feasible.)
J t + 1 ( p ) = 1 K p ∑ q ∈ N ( p ) exp ( − ∥ p − q ∥ 2 2 σ s 2 − ∥ J t ( p ) − J t ( q ) ∥ 2 2 σ r 2 ) I ( q ) J_{t+1}(p)=\frac{1}{K_p}\displaystyle\sum_{q\in N(p)}\exp \bigg(-\frac{ {\|p-q\|}^2}{2\sigma_s^2}-\frac{ {\|J_t(p)-J_t(q)\|}^2}{2\sigma_r^2}\bigg)I(q) J t + 1 ( p ) = K p 1 q ∈ N ( p ) ∑ exp ( − 2 σ s 2 ∥ p − q ∥ 2 − 2 σ r 2 ∥ J t ( p ) − J t ( q ) ∥ 2 ) I ( q )
K p = ∑ q ∈ N ( p ) exp ( − ∥ p − q ∥ 2 2 σ s 2 − ∥ J t ( p ) − J t ( q ) ∥ 2 2 σ r 2 ) K_p=\displaystyle\sum_{q\in N(p)}\exp \bigg(-\frac{ {\|p-q\|}^2}{2\sigma_s^2}-\frac{ {\|J_t(p)-J_t(q)\|}^2}{2\sigma_r^2}\bigg) K p = q ∈ N ( p ) ∑ exp ( − 2 σ s 2 ∥ p − q ∥ 2 − 2 σ r 2 ∥ J t ( p ) − J t ( q ) ∥ 2 )
其中,p p p 是待求点,q q q 是 p p p 周围的所有点的集合, t = 1 , 2 , … , n t=1,2,…,n t = 1 , 2 , … , n , I I I 是原图像, σ s 2 \sigma_s^2 σ s 2 和 σ r 2 \sigma_r^2 σ r 2 分别控制时间权重和距离权重。
σ s 2 \sigma_s^2 σ s 2 and σ r 2 \sigma_r^2 σ r 2 control the spatial and range weights respectively
经典的频谱滤波
R. M. Goldstein and C. L. Werner, ‘Radar interferogram filtering for geophysical applications’, Geophysical Research Letters, vol. 25, no. 21, pp. 4035–4038, 1998, doi: 10.1029/1998GL900033
H ( u , v ) = S { ∣ Z ( u , v ) ∣ } α ⋅ Z ( u , v ) H(u,v) = {S\{ \lvert Z(u,v)\rvert \}}^\alpha \cdot Z(u,v) H ( u , v ) = S {∣ Z ( u , v )∣} α ⋅ Z ( u , v )
优化方向主要有两种,一是优化平滑算子s { ⋅ } s\{\cdot\} s { ⋅ } ,另一则是优化滤波系数⋅ α {\cdot}^{\alpha} ⋅ α 。
baran滤波,使用相干图作为滤波系数
zhao滤波,使用伪相干图作为滤波系数
某个软件,使用bessi0函数生成二维矩阵,在时域点乘后通过FFT转换到频域(时域点乘等于频域卷积),使用频域功率谱计算得到系数,将系数附加到频域后在通过IFFT转换回时域,完成滤波操作。
当然,老老实实写个均值滤波算子在频域做卷积运算,也是可以的,但效果确实不如上述方法。