9. 图像处理的应用 - 欧式视频放大的实现
上一次我介绍了一个计算摄影技术构成的"动作放大器"，它能够高效的将视频中的难以用肉眼察觉的变化分离出来，并在重新渲染过程中进行放大，生成新的视频。这里面的典型代表是欧式视频动作放大。
这一次我们首先回顾这个算法的基本流程，然后简单给大家讲讲每一步的实现方法中的要点。在CMU的原始课程中将欧式视频放大作为了课程作业之一，我也将我的一个简单实现提交到了github上，供各位参考。
一 欧式视频动作放大的基本流程
原作者在文章中用下面这张图形象的描述了算法的实现流程，我已经在上一讲中做了基本的描述：
- 将视频看做一个四维信号I(x,y,c,t)，第3维是颜色，第4维是时间
- 对视频的每一帧进行拉普拉斯金字塔分解（我们前面的章节讲过，参见zhuanlan.zhihu.com/p/53)，这样金字塔的每一层也是一个四维信号，并且代表着原始视频中的不同空间频率。例如金字塔的最上一层代表着最细微的边缘、细节、噪声，最下一层代表着直流信号。
- 接下来对金字塔的每一层信号中的每一个像素点进行带通滤波，这种带通滤波可以在频域上进行，这样我们可以获取到感兴趣的频段的信号
- 滤波后的每一层信号可以乘以特定的放大倍率αk，并与频域滤波前原始信号相加，得到新的金字塔
- 合成新的金字塔中的每一层，得到最终的放大后的视频
下面我展开讲解下每一步中我认为的要点
二、实现步骤
2.1 加载视频
读者可以用任何一种自己熟悉的语言和图像库来加载视频。需要特别注意的是加载过程中的存储空间消耗。有两个特点使得本算法消耗的存储空间很多。
第一是由于需要对时域信号做处理，意味着我们需要加载所有需要处理的视频帧。
第二是是按照论文的建议，需要将图像加载后转换到合适的颜色空间，在CMU原始课程作业要求中，需要将图像数据转换到YIQ颜色空间，数据类型会变为浮点型，存储空间消耗会进一步增加。
下面简单给大家计算一下。
以下面这个视频为例子，每一帧的尺寸是528x592，一共10秒钟共300帧。
face.mp4
将它完全加载到内存中将消耗 528x592x3*300 = 281,318,400 字节内存
如果将数据变为浮点型，则所消耗内存变为281318400*4 = 1125273600 字节，即约1GB字节
由此可见，直接按原始方法实现是很消耗内存的。我看有的朋友在讨论是否能开发手机上的相应应用，我想如果真的要做的话，一定要从工程上去想一些节约内存的方法。在此我给大家一些我实际开发过程中的建议：
1. 确定感兴趣的图像区域，仅对每一帧此区域中的内容做处理
我们假定现在需要开发一个“魔镜”应用，只需要人站在魔镜面前就可以显示出心跳次数数据? 那么我们只去人的脸部或额头一小块区域进行算法的处理。如果这个区域大小只是50x50，那么按上述计算我们仅仅消耗：50*50*3*300*4 = 9,000,000 字节，即9MB
2. 选取合适的视频长度
有时候并不需要很长时间的数据，可以仅仅选取较短的时间的数据用于计算。这样也能节约很多时间
2.2 构建视频金字塔
构建视频金字塔
构建视频金字塔的第一步是构建图像的金字塔，这一点我已经在第5讲，图像采样与金字塔中讲过，这里给大家回忆一下：
图像金字塔构建算法
从图像金字塔中恢复图像
金字塔的层数可以根据实际的输入视频和实际应用需要而调整。
很多时候我们想形象的展示金字塔构建的成果，在OpenCV所带例程中有一段代码做得特别好，这里我稍加整理作为了一个函数提出，并分解和显示了视频中的第7帧，请注意直流帧（就是最小的那幅图像）的颜色显得非常奇怪，因为我们都是在YIQ空间所做的分解。
def testShowPyramid(image, maxLevel):
rows, cols = image.shape[0:2]
imgPyramid = buildLaplacianPyramid(image, maxLevel)
composite_image = np.zeros((rows, cols + cols // 2, 3), dtype=np.double)
composite_image[:rows, :cols, :] = normFloatArray(imgPyramid[0])
i_row = 0
for p in imgPyramid[1:]:
n_rows, n_cols = p.shape[:2]
composite_image[i_row:i_row + n_rows, cols:cols + n_cols] = normFloatArray(p)
i_row += n_rows
plt.figure(figsize=(15,15))
plt.title("Laplacian Pyramid Show for Frame 7")
plt.imshow(composite_image)
plt.show()
当构建图像金字塔的函数写好后，构建视频金字塔就非常容易了。下面我展示了相关的函数，及构造视频金字塔和从视频金字塔中重建帧的结果，可以看到重建后和重建前的视频帧几乎一致。
def buildVideoLapPyr(frames, maxLevel):
"""
Build Laplacian pyramid for input video frames
Parameters
----------
frames: input video frames
maxLevel: upper limit of the Laplician pyramid layers
Returns
-------
Laplacian pyramid for input video frames, which is a list of
videos that each video mapping to a layer in the laplacian pyramid.
So each video has the shape (frameCount, videoFrameHeight, videoFrameWidth, channelCount).
"""
pyr0=buildLaplacianPyramid(frames[0], maxLevel)
realMaxLevel=len(pyr0)
resultList=[]
for i in range(realMaxLevel):
curPyr = np.zeros([len(frames)]+list(pyr0[i].shape), dtype=np.float32)
resultList.append(curPyr)
for fn in range(len(frames)):
pyrOfFrame = buildLaplacianPyramid(frames[fn], maxLevel)
for i in range(realMaxLevel):
resultList[i][fn]=pyrOfFrame[i]
return resultList
def recreateVideoFromLapPyr(pyrVideo):
"""
Recreate video from input video Laplacian Pyramid
Parameters
----------
pyrVideo: input Laplacian Pyramid for video, returned by buildVideoLapPyr
Returns
-------
Video recreated from inut Laplacian Pyramid
"""
maxLevel=len(pyrVideo)
fNumber, H, W, chNum=pyrVideo[0].shape
videoResult=np.zeros(pyrVideo[0].shape, dtype=np.float32)
for fn in range(videoResult.shape[0]):
framePyr=[pyrVideo[i][fn] for i in range(maxLevel)]
videoResult[fn]=recreateImgsFromLapPyr(framePyr)
return videoResult
2.3 时域信号的频域滤波
2.3.1 观察时域信号
让我们先看看视频金字塔中第1层和第4层的时域信号长什么样，我选取了视频中人像额头上一点来观察：
2.3.2 选择滤波器
对于我们当前这个观察人脸颜色随心跳变换而变换的简单应用来说，到底应该选择什么样的滤波器呢？作者的原论文中有这样一段话，我直接贴到这里：
The Choice of filter is generally application dependent. For motion magnification, a filter with a broad passband is prefered; for color amplification of blood flow, a narrow passband produce a more noise-free result.
即对于颜色放大，我们应该用简单的窄带带通滤波器。而对于一些动作放大的应用，我们更倾向于采用带宽更宽的滤波器。下面我展示了论文中提到的几种滤波器的形态：
我从github.com/flyingzhao/P中找到了相关滤波器的实现，这里稍加改编后采用，并构建了视频金字塔的滤波器
def temporal_ideal_filter(tensor,low,high,fps,axis=0):
"""
Apply Ideal bandpass filter on input data on specified axis.
Parameters
----------
tensor : data to be filered
low, hight: the cut frequency of the bandpass filter
fps : sample rate
axis : The axis of the input data array along which to apply the linear filter.
The filter is applied to each subarray along this axis.
Returns
-------
None
"""
fft=fftpack.fft(tensor,axis=axis)
frequencies = fftpack.fftfreq(tensor.shape[0], d=1.0 / fps)
bound_low = (np.abs(frequencies - low)).argmin()
bound_high = (np.abs(frequencies - high)).argmin()
if (bound_low==bound_high) and (bound_high<len(fft)-1):
bound_high+=1
fft[:bound_low] = 0
fft[bound_high:-bound_high] = 0
fft[-bound_low:] = 0
iff=fftpack.ifft(fft, axis=axis)
return np.abs(iff)
def idealFilterForVideoPyr(videoPyr, low, high, fps, roi=None):
"""
Apply Ideal bandpass filter on input video pyramid
Parameters
----------
videoPyr : video pyramid
low, hight: the cut frequency of the bandpass filter
fps : sample rate
roi : if specified, only filter roi of frames. [TODO]Not implemented yet.
Returns
-------
None
"""
resultVideoPyr=[]
for layer in range(len(videoPyr)):
filteredPyr = temporal_ideal_filter(videoPyr[layer], low, high, fps, axis=0)
resultVideoPyr.append(filteredPyr)
return resultVideoPyr让我们观察视频金字塔中人像额头处某个像素的时域信号在滤波前后的状况，这里显示了第1层和第4层的信号：
可以看到滤波后信号被滤除得非常干净，成为了正弦（余弦）信号
2.3.3 放大信号
接下来对每一层信号做不同程度的放大。
放大过程很简单，是由一定比例的滤波后的视频 + 原始视频信号，论文里面还提到可以对亮度通道和色度通道进行不同比例的放大。
下面展示了放大前后放大后的视频信号，它具有明显的规律性，但又不像之前滤波后的图像那么干净的正弦（余弦）信号。
2.3.4 重建视频
重建视频
调用前面写好的视频金字塔重建函数可以重建视频，具体代码可以参看我给出的代码链接，这里是我得到的视频：
基本结果
2.3.5 比较不同的滤波器
如果前面不选用理想窄带滤波器，而是选用Butterworth滤波器，会是什么结果呢？下面是我的处理结果，可以看到视频完全紊乱了，也许需要仔细的调参才能勉强得到一个更好的视频。当然按照原论文说法，颜色放大的应用需要的是窄带滤波器，本来就不应该使用Butterworth这样的滤波器
Buttworth滤波器结果
2.3.6 采用不同的金字塔层数
下面看看分别采用2层金字塔和6层金字塔的结果：
2层金字塔结果
6层金字塔结果
可以看到当金字塔层数不够时，相当于空间滤波不够，此时大量的噪声被放大了。而金字塔层数太多，放大效应则不明显。
三、总结
这一篇简单讲解了实现欧式视频放大中的一些关键步骤，跟这一系列专题文章相关的Notebook可以从github.com/yourwanghao/获取
也可以通过这个页面直接查看结果：Jupyter Notebook Viewer
参考资料：
这一篇文章的绝大部分素材来自于以下资料，尤其是其中的课程作业介绍
- Eulerian Video Magnification for Revealing Subtle Changes in the World
- graphics.cs.cmu.edu/cou 中的 第二次作业 graphics.cs.cmu.edu/cou
- github.com/yourwanghao/