专题二:欧拉视频放大(基于相位)---原理解析

前言

由于欧拉线性放大方法会不可避免的放大噪声,所以放大率和放大效果都会受到影响。而基于相位的放大方法不会放大噪声,只会平移噪声,代价是算法的复杂度增加,所以在追求更好的放大效果时,我们可以考虑基于相位的放大。

实施的步骤

空间滤波、提取相位

我们知道图像傅里叶变换的幅度谱代表信号的强度,相位谱代表信号的位置。如果我们直接对相位进行操作,那么整幅图像都会受到影响。
所以我们需要构造在空间域和频率域都具有局部性的滤波器,提取低通的频率作用于有限的空间。Gabor滤波器是其中的典型,但是这里并不能使用Gabor滤波器。因为Gabor滤波器没有良好的重构特性。具体原因待我下回分解。
这里使用的是可操纵金子塔分解,即方向的可操纵性。可操纵金字塔的具体原理待我下下回分解,哈哈。
我们并没有使用金字塔的方向可操纵性,也就是说各个方向的权值是相等的。使用的只是可操纵金字塔的子带非混叠性和充当正交相位滤波器的性质。下面是我取得某个尺度上的一组可操纵滤波器的基。
在这里插入图片描述
注意这个是频域上的表现,是实数,我们对以上某个基再做傅里叶反变换,如下图(左图是实部,右图是虚部)
在这里插入图片描述
可以看见其实部为偶函数,虚部为奇函数。这是因为频域一半有数据一半没数据,假设是被x轴对半分吧,那么空间某点就会在对角的相同位置生成共轭项。又因为频域关于y轴偶对称,那么这两个位置的共轭项也会相等,即空间点不仅对角位置成共轭,关于x轴也成共轭。那么我们就得到了关于x轴偶对称的实部,关于x轴奇对称的虚部。这个x轴是我假定的,其实就是有数无数的分界, 垂直于基滤波器的方向y,注意关于y轴都是偶对称的。
那么为什么在空间中实部为偶对称的,虚部为奇对称的滤波器可以提取相位呢,这里仅提供个人理解。在程序里作者在频域用以上基滤波器乘以信号再求反傅里叶变换即得到基函数的空间表示,其实部对应于偶对称滤波器,虚部对应于奇滤波器。我们知道任何信号都可以用傅里叶级数展开,即
在这里插入图片描述
A ( ω ) = a ( ω ) + j b ( ω ) A(\omega)=a(\omega)+jb(\omega) A(ω)=a(ω)+jb(ω)
而任何一个信号都可以分解成偶信号与奇信号的叠加,我们以一维信号为例, f ( x ) = f e ( x ) + f o ( x ) f(x)=f_e(x)+f_o(x) f(x)=fe(x)+fo(x) { f e ( x ) = f ( x ) + f ( − x ) 2 f o ( x ) = f ( x ) − f ( − x ) 2 \begin{cases} f_e(x)=\frac{f(x)+f(-x)}{2} \\[2ex] f_o(x)=\frac{f(x)-f(-x)}{2} \end{cases} fe(x)=2f(x)+f(x)fo(x)=2f(x)f(x)
而信号通过偶滤波器会消除 f o ( x ) f_o(x) fo(x)信号,只保留 f e ( x ) f_e(x) fe(x)信号,通过奇滤波器会消除 f e ( x ) f_e(x) fe(x)信号,只保留 f o ( x ) f_o(x) fo(x)信号。而作者的基滤波器乘以信号再求反傅里叶变换相当于信号先求反傅里叶变换再与奇偶滤波器卷积。

图像在空间是实信号非奇非偶,那么其频域就是复信号,且实部偶对称,虚部奇对称,我们把该复信号拆开得到频域实部 a ( ω ) a(\omega) a(ω)和虚部 b ( ω ) b(\omega) b(ω)两个信号,分别进行傅里叶反变换,得到空域实偶信号和实奇信号。

滤波器在频域是实信号不完全的偶信号,那么滤波器的空域就是如上图的实部偶对称,虚部奇对称。

注意刚刚频域的复信号对应的空间域的两个实信号其实并没有分离是叠加在一起的一个实数,只是这么理解而已,到这里经过空域的实部虚部滤波器的过滤才得以分开,整体以复数形式存在。即实信号 f ( x , y ) f(x,y) f(xy)经过与复数滤波器的卷积得到 F e ( x , y ) + j F o ( x , y ) F_e(x, y)+jF_o(x,y) Fe(xy)+jFo(xy)此时空间某个点就是复数形式存在了。空间域上的函数 F o ( x , y ) F_o(x,y) Fo(xy)对应于频域的 a ( ω ) a(\omega) a(ω), F e ( x , y ) F_e(x,y) Fe(xy)对应与频域的 b ( ω ) b(\omega) b(ω)注意是二维的,在频域某一个点的相位是求虚部比实部的反正切,那么类比的空域某个点x的相位也是虚部比实部的反正切。我忘了在哪看过广义相位了,不知道是不是我理解的这样。注意该点x的相位其实是各个频率的相位在该点的共同作用。
该基方向就是y轴方向
这样我们就达到了引入可操纵金字塔的目的,求解每个空间位置的相位。

时域滤

接着就是进行时域带通滤波,与线性欧拉放大大同小异,不多说。注意是对第N帧与第1帧的相位差进行滤波。

放大

放大前会有一个振幅加权相位的过程。即先根据周围复值距离目标点的远近对复值的振幅进行高斯加权,再根据复值的振幅即模的大小对相位加权。这样做的原因是相位在局部区域具有一定的连续性,考虑周围信号的影响可以减小误差。
然后我们施加一个放大因子magPhase,这个值是自己设定的。

合成

在这里插入图片描述
其中总和包含金字塔中的所有尺度和方向,产生重建图像。

数学原理推导

在这里插入图片描述
在这里插入图片描述
在这里插入图片描述
在这里插入图片描述
这是作者提到的参考文献里的福值可操纵金字塔的构成框图及公式,我已验证过,但是关于系数有一些疑问。

[himask, lomaskPrev] = getRadialMaskPair(rVals(1), log_rad, twidth);%在getRadialMaskPair取了以2为底的对数
filters{count} = himask;%himask内边界衰减速度变快,外边界衰减速度变慢  filters里有count=1个元素,该元素是一个矩阵himask
count = count + 1;
for k = 2:max(size(rVals))%k从276层循环对应于2^-12^-6
   [himask, lomask] = getRadialMaskPair(rVals(k), log_rad, twidth);
   radMask = himask.*lomaskPrev;%相乘相当于嵌套,以lomaskPrev的内边界为界限,外面等于lomaskPrev(正弦)的值,内部等于himask(余弦)的值
   for j = 1:orientations
      anglemask = getAngleMask(j, orientations, angle);
      filters{count} = radMask.*anglemask/2;
      count = count + 1;
   end
   
   lomaskPrev = lomask;

在该参考文献里说由于发生了采样,所以需要乘以系数2,得到下一次的低通滤波器的幅值是上一次的两倍。但我还没有看到采样在哪里。程序里似乎也没递乘2.

  • 5
    点赞
  • 47
    收藏
    觉得还不错? 一键收藏
  • 5
    评论
以下是基于欧拉方程的Rayleigh-Taylor不稳定性流动Python代码: ```python import numpy as np import matplotlib.pyplot as plt # 定义常数和参数 g = 9.81 # 重力加速度 rho1 = 1.0 # 上层液体密度 rho2 = 2.0 # 下层液体密度 p1 = 1.0 # 上层液体压力 p2 = 2.0 # 下层液体压力 Lx = 1.0 # 液面长度 Ly = 1.0 # 液面宽度 Nx = 101 # 网格数目 Ny = 101 dx = Lx / (Nx - 1) # 网格间距 dy = Ly / (Ny - 1) dt = 0.01 # 时间步长 tmax = 5.0 # 最大时间 # 初始化网格 x = np.linspace(0, Lx, Nx) y = np.linspace(0, Ly, Ny) X, Y = np.meshgrid(x, y) u = np.zeros((Nx, Ny)) # x方向速度 v = np.zeros((Nx, Ny)) # y方向速度 p = np.zeros((Nx, Ny)) # 压力 rho = np.zeros((Nx, Ny)) # 密度 h = np.zeros((Nx, Ny)) # 液面高度 # 设置初始液面高度 h[:, :] = 0.5 * (rho1 - rho2) / (rho1 + rho2) * Ly * np.cos(np.pi * X / Lx) # 循环计算 for n in range(int(tmax / dt)): # 计算密度和压力 rho[:, :] = rho1 * (h > 0) + rho2 * (h <= 0) p[:, :] = p1 * (h > 0) + p2 * (h <= 0) # 计算速度 u[1:-1, 1:-1] = u[1:-1, 1:-1] - dt * (p[1:-1, 2:] - p[1:-1, :-2]) / (2 * rho1 * dx) \ - g * dt * (rho[1:-1, 1:-1] - rho[1:-1, :-2]) / (2 * rho1 * dy) v[1:-1, 1:-1] = v[1:-1, 1:-1] - dt * (p[2:, 1:-1] - p[:-2, 1:-1]) / (2 * rho1 * dy) \ - g * dt * (rho[1:-1, 1:-1] - rho[:-2, 1:-1]) / (2 * rho1 * dx) # 更新液面高度 h[1:-1, 1:-1] = h[1:-1, 1:-1] - dt * (u[1:-1, 2:] - u[1:-1, :-2]) / dx \ - dt * (v[2:, 1:-1] - v[:-2, 1:-1]) / dy # 边界条件 u[:, 0] = 0 u[:, -1] = 0 u[0, :] = 0 u[-1, :] = 0 v[:, 0] = 0 v[:, -1] = 0 v[0, :] = 0 v[-1, :] = 0 h[0, :] = h[1, :] h[-1, :] = h[-2, :] h[:, 0] = h[:, 1] h[:, -1] = h[:, -2] # 绘图 if n % 10 == 0: plt.clf() plt.contourf(X, Y, h, cmap='coolwarm') plt.colorbar() plt.title('t = {:.2f}'.format(n * dt)) plt.xlabel('x') plt.ylabel('y') plt.pause(0.01) plt.show() ``` 请注意,此代码仅提供参考,并且可能需要进行一些修改才能适应特定的问题。
评论 5
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值