當(dāng)前位置：首頁 > 编程资源 > 编程问答 >内容正文

编程问答

Eulerian Video Magnification

發(fā)布時(shí)間：2024/3/24 编程问答 39 豆豆

生活随笔收集整理的這篇文章主要介紹了 Eulerian Video Magnification 小編覺得挺不錯(cuò)的,現(xiàn)在分享給大家,幫大家做個(gè)參考.

原文：http://www.hahack.com/codes/eulerian-video-magnification/

引言

人類的視覺感知存在有限的感知域。對(duì)于超出感知域的變化，我們無法感知。然而，這類信號(hào)卻可能蘊(yùn)藏著驚人的秘密。

比如，血液循環(huán)使得人體的皮膚發(fā)生細(xì)微的周期性變化，這個(gè)裸眼無法感知的變化卻和人的心率非常吻合。2011 年，MIT 的一個(gè)亞裔學(xué)生 Mingzhe Poh 就利用這個(gè)微弱的信號(hào)設(shè)計(jì)了一個(gè)“魔鏡”[1]?—— 不僅能照出你的模樣，還能測(cè)出你的心率。

Mingzhe Poh 的這面神奇的鏡子的原理是利用了血液在人體內(nèi)流動(dòng)時(shí)光線的變化?[2]?：心臟跳動(dòng)時(shí)血液會(huì)通過血管，通過血管的血液量越大，被血液吸收的光線也越多，人皮膚表面反射的光線就越少。因此，通過對(duì)圖像的時(shí)頻分析就可以估算出心率。

再比如，樂器之所以會(huì)發(fā)出聲音，是因?yàn)樗l(fā)聲的部分在彈奏過程中發(fā)生了有規(guī)律的形變，而這個(gè)形變的振幅對(duì)應(yīng)著樂器發(fā)聲的響度，快慢對(duì)應(yīng)著樂器的音高。微弱信號(hào)所蘊(yùn)藏的信息量是如此重大，無怪乎禪語有云：

一花一世界，一葉一菩提。

既然如此，能否將影像中的這些肉眼觀察不到的變化“放大”到裸眼足以觀察的幅度呢？這就是本文將要重點(diǎn)討論的問題。

在接下來的篇幅中，我將首先追溯最早的一個(gè)放大變化的實(shí)驗(yàn)——卡文迪許實(shí)驗(yàn)，然后引出適用于現(xiàn)代計(jì)算機(jī)的兩種視角下的影像放大方法。這兩種視角分別稱為拉格朗日視角（Lagrangian Perspective）和歐拉視角（Eulerian Perspective）。最后我將重點(diǎn)探討歐拉視角的算法實(shí)現(xiàn)細(xì)節(jié)。我所實(shí)現(xiàn)的一個(gè)歐拉影像放大算法程序在 Github 上開源。

附圖 1?Henry Cavendish

最早的放大：卡文迪許扭秤實(shí)驗(yàn)

我所能追溯到的最早的將變化“放大”的實(shí)驗(yàn)是 1797 年卡文迪許（H.Cavendish）的經(jīng)典實(shí)驗(yàn)——扭秤實(shí)驗(yàn)。卡文迪許實(shí)驗(yàn)是第一個(gè)在實(shí)驗(yàn)室里完成的測(cè)量?jī)蓚€(gè)物體之間萬有引力的實(shí)驗(yàn)，并且第一個(gè)準(zhǔn)確地求出了萬有引力常數(shù)和地球質(zhì)量。

實(shí)驗(yàn)的裝置由約翰·米切爾設(shè)計(jì)，由兩個(gè)重達(dá)350磅的鉛球和扭秤系統(tǒng)組成。

卡文迪許用兩個(gè)質(zhì)量一樣的鉛球分別放在扭秤的兩端。扭秤中間用一根韌性很好的鋼絲系在支架上，鋼絲上有個(gè)小鏡子。用準(zhǔn)直的細(xì)光束照射鏡子，細(xì)光束反射到一個(gè)很遠(yuǎn)的地方，標(biāo)記下此時(shí)細(xì)光束所在的點(diǎn)。用兩個(gè)質(zhì)量一樣的鉛球同時(shí)分別吸引扭秤上的兩個(gè)鉛球。由于萬有引力作用。扭秤微微偏轉(zhuǎn)。但細(xì)光束所反射的遠(yuǎn)點(diǎn)卻移動(dòng)了較大的距離。他用此計(jì)算出了萬有引力公式中的常數(shù)?GG?。

卡文迪許實(shí)驗(yàn)取得成功的原因，是將不易觀察的微小變化量，轉(zhuǎn)化（放大）為容易觀察的顯著變化量，再根據(jù)顯著變化量與微小量的關(guān)系算出微小的變化量。

卡文迪許的實(shí)驗(yàn)給了我們一個(gè)啟示：放大變化，就是要解決以下兩個(gè)問題：

何為“變”?—— 如何找出不易觀察的微小變化量；

放大“變”?—— 如何放大這個(gè)變化量，使之肉眼可見。

不過，卡文迪許的這個(gè)實(shí)驗(yàn)需要借助一個(gè)龐大的扭秤裝置，并不能直接用來放大影像中的變化。對(duì)于生活在二十一世紀(jì)的我們，最理想的方式當(dāng)然是要借助計(jì)算機(jī)這個(gè)神器了。接下來將介紹兩種現(xiàn)代的技術(shù)方案，能夠讓計(jì)算機(jī)為我們放大影像中細(xì)微的變化，從而使我們具有這樣一對(duì)火眼金睛，去發(fā)現(xiàn)大自然隱藏的秘密。

拉格朗日視角

附圖 2?六祖慧能

時(shí)有風(fēng)吹幡動(dòng)。一僧云：風(fēng)動(dòng)。一僧云：幡動(dòng)。議論不已。能進(jìn)曰：不是風(fēng)動(dòng)，不是幡動(dòng)，仁者心動(dòng)。一眾駭然。

六祖慧能在初見五祖的時(shí)候，恰逢有風(fēng)吹來，吹得幡動(dòng)。于是一個(gè)和尚說是風(fēng)在動(dòng)，另一個(gè)和尚說是幡在動(dòng)，而慧能卻一語道破：不是風(fēng)動(dòng)，也不是幡動(dòng)，而是你的心在動(dòng)。

看待一樣?xùn)|西，視角不同，得出的結(jié)論也就不同。正如看待生活中的“變”，視角不同，也會(huì)得到不同的結(jié)果。

所謂拉格朗日視角，就是從跟蹤圖像中感興趣的像素（粒子）的運(yùn)動(dòng)軌跡的角度著手分析。打個(gè)比方，假如我們要研究河水的流速，我們坐上一條船，順流而下，然后記錄這條船的運(yùn)動(dòng)軌跡。

何為“變” —— 感興趣的像素點(diǎn)隨著時(shí)間的運(yùn)動(dòng)軌跡，這類像素點(diǎn)往往需要借助人工或其他先驗(yàn)知識(shí)來輔助確定；

放大“變” —— 將這些像素點(diǎn)的運(yùn)動(dòng)幅度加大。

2005 年，Liu 等人最早提出了一種針對(duì)影像的動(dòng)作放大技術(shù)[3]，該方法首先對(duì)目標(biāo)的特征點(diǎn)進(jìn)行聚類，然后跟蹤這些點(diǎn)隨時(shí)間的運(yùn)動(dòng)軌跡，最后將這些點(diǎn)的運(yùn)動(dòng)幅度加大。

然而，拉格朗日視角的方法存在以下幾點(diǎn)不足：

需要對(duì)粒子的運(yùn)動(dòng)軌跡進(jìn)行精確的跟蹤和估計(jì)，需要耗費(fèi)較多的計(jì)算資源；
對(duì)粒子的跟蹤是獨(dú)立進(jìn)行的，缺乏對(duì)整體圖像的考慮，容易出現(xiàn)圖像沒有閉合，從而影響放大后的效果；
對(duì)目標(biāo)物體動(dòng)作的放大就是修改粒子的運(yùn)動(dòng)軌跡，由于粒子的位置發(fā)生了變化，還需要對(duì)粒子原先的位置進(jìn)行背景填充，同樣會(huì)增加算法的復(fù)雜度。

歐拉視角

不同于拉格朗日視角，歐拉視角并不顯式地跟蹤和估計(jì)粒子的運(yùn)動(dòng)，而是將視角固定在一個(gè)地方，例如整幅圖像。之后，假定整幅圖像都在變，只是這些變化信號(hào)的頻率、振幅等特性不同，而我們所感興趣的變化信號(hào)就身處其中。這樣，對(duì)“變”的放大就變成了對(duì)感興趣頻段的析出和增強(qiáng)。打個(gè)比方，同樣是研究河水的流速，我們也可以坐在岸邊，觀察河水經(jīng)過一個(gè)固定的地方時(shí)的變化，這個(gè)變化可能包含很多和水流本身無關(guān)的成分，比如葉子掉下水面激起的漣漪，但我們只關(guān)注最能體現(xiàn)水流速的部分。

何為“變” —— 整個(gè)場(chǎng)景都在變，而我們所感興趣的變化信號(hào)藏在其中；

放大“變” —— 通過信號(hào)處理手段，將感興趣的信號(hào)分離，并進(jìn)行增強(qiáng)。

2012 年， Wu 等人從這個(gè)視角著手，提出了一種稱為歐拉影像放大技術(shù)（Eulerian Video Magnification）的方法[4]，其流程如下：

空間濾波。將視頻序列進(jìn)行金字塔多分辨率分解；

時(shí)域?yàn)V波。對(duì)每個(gè)尺度的圖像進(jìn)行時(shí)域帶通濾波，得到感興趣的若干頻帶；

放大濾波結(jié)果。對(duì)每個(gè)頻帶的信號(hào)用泰勒級(jí)數(shù)來差分逼近，線性放大逼近的結(jié)果；

合成圖像。合成經(jīng)過放大后的圖像。

在下一節(jié)我們將重點(diǎn)探討 Wu 等人所提出的這種線性的歐拉影像放大技術(shù)。之所以加上“線性”這個(gè)修飾詞，是因?yàn)?Wadhwa 等人在 2013 年對(duì)這項(xiàng)技術(shù)進(jìn)行了改進(jìn)，提出了基于相位的影像動(dòng)作處理技術(shù)[5]。基于相位的歐拉影像放大技術(shù)在放大動(dòng)作的同時(shí)不會(huì)放大噪聲，而是平移了噪聲，因而可以達(dá)到更好的放大效果。不過對(duì)它的討論超出了本文的篇幅，感興趣的讀者可以自己找來他們的 paper 閱讀。

算法細(xì)節(jié)

空間濾波

如前面所說，歐拉影像放大技術(shù)（以下簡(jiǎn)稱 EVM ）的第一步是對(duì)視頻序列進(jìn)行空間濾波，以得到不同的空間頻率的基帶。這么做是因?yàn)?#xff1a;

有助于減少噪聲。圖像在不同空間頻率下呈現(xiàn)出不同的SNR（信噪比）。一般來說，空間頻率越低，信噪比反而越高。因此，為了防止失真，這些基帶應(yīng)該使用不同的放大倍數(shù)。最頂層的圖像，即空間頻率最低、信噪比最高的圖像，可使用最大的放大倍數(shù)，下一層的放大倍數(shù)依次減小；

便于對(duì)圖像信號(hào)的逼近。空間頻率較高的圖像（如原視頻圖像）可能難以用泰勒級(jí)數(shù)展開來逼近。因?yàn)樵谶@種情況下，逼近的結(jié)果就會(huì)出現(xiàn)混淆，直接放大就會(huì)出現(xiàn)明顯失真1?1對(duì)于這種情況，論文通過引入一個(gè)空間波長(zhǎng)下限值來減少失真。如果當(dāng)前基帶的空間波長(zhǎng)小于這個(gè)下限值，就減少放大倍數(shù)。。

由于空間濾波的目的只是簡(jiǎn)單的將多個(gè)相鄰的像素“拼”成一塊，所以可以使用低通濾波器來進(jìn)行。為了加快運(yùn)算速度，還可以順便進(jìn)行下采樣操作。熟悉圖像處理操作的朋友應(yīng)該很快可以反應(yīng)出來：這兩個(gè)東西的組合就是金字塔。實(shí)際上，線性的 EVM 就是使用拉普拉斯金字塔或高斯金字塔來進(jìn)行多分辨率分解。

時(shí)域?yàn)V波

得到了不同空間頻率的基帶后，接下來對(duì)每個(gè)基帶都進(jìn)行時(shí)域上的帶通濾波，目的是提取我們感興趣的那部分變化信號(hào)。

例如，如果我們要放大的心率信號(hào)，那么可以選擇 0.4 ～ 4 Hz （24～240 bpm ）進(jìn)行帶通濾波，這個(gè)頻段就是人的心率的范圍。

不過，帶通濾波器有很多種，常見的就有理想帶通濾波器、巴特沃斯（Butterworth）帶通濾波器、高斯帶通濾波器，等等。應(yīng)該選擇哪個(gè)呢？這得根據(jù)放大的目的來選擇。如果需要對(duì)放大結(jié)果進(jìn)行后續(xù)的時(shí)頻分析（例如提取心率、分析樂器的頻率），則應(yīng)該選擇窄通帶的濾波器，如理想帶通濾波器，因?yàn)檫@類濾波器可以直接截取出感興趣的頻段，而避免放大其他頻段；如果不需要對(duì)放大結(jié)果進(jìn)行時(shí)頻分析，可以選擇寬通帶的濾波器，如 Butterworth 帶通濾波器，二階 IIR 濾波器等，因?yàn)檫@類濾波器可以更好的減輕振鈴現(xiàn)象。

放大和合成

經(jīng)過前面兩步，我們已經(jīng)找出了“變”的部分，即解決了何為“變”這個(gè)問題。接下來我們來探討如何放大“變”這個(gè)問題。

一個(gè)重要的依據(jù)是：上一步帶通濾波的結(jié)果，就是對(duì)感興趣的變化的逼近。接下來將證明這個(gè)觀點(diǎn)。

這里以放大一維的信號(hào)為例，二維的圖像信號(hào)與此類似。

假定現(xiàn)在有個(gè)信號(hào)?I(x,t)I(x,t)?，在任意時(shí)刻?tt?，有：

eq: 1 ?

I(x,t)I(x,0)=f(x+δ(t))=f(x),t>0,t=0I(x,t)=f(x+δ(t)),t>0I(x,0)=f(x),t=0

其中，δ(t)δ(t)?是變化信號(hào)，在本例中就是一個(gè)位移函數(shù)。

我們希望得到這個(gè)變化放大?αα?倍后的結(jié)果，即：

eq: 2 ?

I^(x,t)=f(x+(1+α)δ(t))I^(x,t)=f(x+(1+α)δ(t))

為了將變化的部分分離出來，我們用一階泰勒級(jí)數(shù)展開來逼近公式 1 表示的信號(hào)2?2想起了高數(shù)老師在第一節(jié)課所說的話：高等代數(shù)，也就是微積分，研究的就是一個(gè)字：“變”。：

eq: 3 ?

I(x,t)≈f(x)+δ(t)?f(x)?xI(x,t)≈f(x)+δ(t)?f(x)?x

公式 2 中標(biāo)藍(lán)的部分恰好就是變化的部分，而這個(gè)部分又和上一步的帶通濾波有著重要的聯(lián)系。下面分兩種情況討論。

理想情況

我們先考慮一種理想情況：假如所有的變化信號(hào)?δ(t)δ(t)?的頻率范圍恰好是我們進(jìn)行帶通濾波時(shí)所選的頻帶范圍，那么帶通濾波結(jié)果?B(x,t)B(x,t)?應(yīng)該恰好等于公式 2 中標(biāo)藍(lán)的部分，即：

eq: 4 ?

B(x,t)=δ(t)?f(x)?xB(x,t)=δ(t)?f(x)?x

對(duì)公式 2 所逼近的信號(hào)進(jìn)行放大，就是將這個(gè)變化的部分乘以一個(gè)放大倍數(shù)?αα?，再加回原來的信號(hào)中。即：

eq: 5 ?

I~(x,t)=I(x,t)+αB(x,t)I~(x,t)=I(x,t)+αB(x,t)

聯(lián)立公式 2～4 ，可以得到

eq: 6 ?

I~(x,t)≈f(x)+(1+α)δ(t)?f(x)?xI~(x,t)≈f(x)+(1+α)δ(t)?f(x)?x

在這種理想情況下，這個(gè)?I~(x,t)I~(x,t)?約等于我們希望得到的?I(x,t)I(x,t)?，即：

eq: 7 ?

I~(x,t)≈f(x+(1+α)δ(t))I~(x,t)≈f(x+(1+α)δ(t))

下圖演示了使用上面的方法將一個(gè)余弦波放大?αα?倍的過程和結(jié)果。其中，黑色的曲線表示原信號(hào)?f(x)f(x)?，藍(lán)色的曲線表示變化后的信號(hào)?f(x+δ)f(x+δ)?，青色的曲線表示對(duì)這個(gè)信號(hào)的泰勒級(jí)數(shù)逼近?f(x)+δ(t)?f(x)?xf(x)+δ(t)?f(x)?x，綠色的曲線表示我們分離出來的變化的部分。我們將這個(gè)部分放大?αα?倍再加回原信號(hào)就得到放大后的信號(hào)，圖中紅色的曲線表示這個(gè)放大后的信號(hào)?f(x)+(1+α)B(x,t))f(x)+(1+α)B(x,t))?。

非理想情況

不過，有些時(shí)候，我們并沒有那么幸運(yùn)——變化信號(hào)?δ(t)δ(t)?的頻率范圍超出了我們所選的頻段范圍。這種情況下，應(yīng)用帶通濾波意味著只是保留了一部分的變化信號(hào)，而其他頻率超出范圍的信號(hào)將會(huì)被減弱。因此，我們用?γk(t)γk(t)?來表示在?tt?時(shí)刻變化第?kk?個(gè)變化信號(hào)減弱的倍數(shù)（0≤γk≤10≤γk≤1）3?3它的值和所選的帶通濾波器有關(guān)。實(shí)際上論文原作者考慮這個(gè)情況僅僅是出于論證上的嚴(yán)密，在實(shí)現(xiàn)時(shí)我們不需要計(jì)算它，只需要得到 B(x,t) ，而這個(gè)值是上一步的帶通結(jié)果。，則有：

eq: 8 ?

B(x,t)=∑kγkδk(t)?f(x)?xB(x,t)=∑kγkδk(t)?f(x)?x

之后我們又要對(duì)它乘以放大倍數(shù)?αα?。既然兩步都是乘以一個(gè)倍數(shù)，為了方便起見，我們干脆把這兩個(gè)線性變化合為一步，即讓放大倍數(shù)?αk=γkααk=γkα?，則：

eq: 9 ?

I~(x,t)≈f(x+∑k(1+αk)δk(t))I~(x,t)≈f(x+∑k(1+αk)δk(t))

放大倍數(shù)限制

線性的 EVM 方法會(huì)在放大動(dòng)作變化的同時(shí)放大噪聲，為了避免造成太大的失真，可以設(shè)置一個(gè)合理的放大倍數(shù)限制。假定信號(hào)的空間波長(zhǎng)為?λ=2πωλ=2πω?，這個(gè)限制可以用公式 9 來表示：

eq: 10 ?

(1+α)δ(t)<λ8(1+α)δ(t)<λ8

當(dāng)超出這個(gè)邊界的時(shí)候，我們可以讓?αα?維持在一個(gè)邊界值，如下圖所示：

不過，如果要放大的是顏色的變化，那么從視覺上看并不會(huì)很受影響（或者說這種顏色的失真就是我們想要的），這時(shí)候就可以不用對(duì)?αα?進(jìn)行限制，或者使用一個(gè)更小的空間波長(zhǎng)下限值。

算法實(shí)現(xiàn)

這一節(jié)將介紹我使用 OpenCV 實(shí)現(xiàn)線性歐拉影像放大算法的心得。

本文的源碼遵守 LGPL v3 協(xié)議，但這個(gè)技術(shù)已由原作者申請(qǐng)了專利，因此請(qǐng)勿直接用于商業(yè)用途。
動(dòng)作放大的模塊引用了 Yusuke Tomoto 的?ofxEvm?項(xiàng)目的相關(guān)代碼。這個(gè)項(xiàng)目基于 EVM 實(shí)現(xiàn)了對(duì)動(dòng)作的放大。對(duì)于初學(xué)者，我建議先閱讀他的代碼。
另外特別感謝?Daniel Ron?和?Aalessandro Gentilini?的大力援助，沒有他們我無法完成這個(gè)程序。

顏色空間轉(zhuǎn)換

顏色空間轉(zhuǎn)換在原論文中只是一筆帶過，但這一步對(duì)于放大動(dòng)作還是非常有用的。在進(jìn)入整個(gè) framework 之前，作者建議先將圖像的顏色空間由 RGB 轉(zhuǎn)換到?YIQ?。YIQ 是 NTSC 電視機(jī)系統(tǒng)所采用的顏色空間，Y?是提供黑白電視及彩色電視的亮度信號(hào)，I?代表 In-phase，色彩從橙色到青色，Q?代表 Quadrature-phase，色彩從紫色到黃綠色。

采用這種顏色空間可以方便在后期使用一個(gè)衰減因子來減少噪聲：對(duì)于只想放大動(dòng)作變化的情況，顏色就應(yīng)該不會(huì)發(fā)生太大變化，所以我們可以用這個(gè)衰減因子來減小放大后的信號(hào)的 I 和 Q 兩個(gè)分量的值。然后再轉(zhuǎn)回 RGB 顏色空間44類似的顏色空間還有 Lab，經(jīng)過我的測(cè)試，使用 Lab 顏色空間也有效。。兩個(gè)轉(zhuǎn)換函數(shù)實(shí)現(xiàn)如下：

void rgb2ntsc(const Mat_<Vec3f>& src, Mat_<Vec3f>& dst)

{

Mat ret = src.clone();

Mat T = (Mat_<float>(3,3) << 1, 1, 1, 0.956, -0.272, -1.106, 0.621, -0.647, 1.703);

T = T.inv(); //here inverse!

for (int j=0; j<src.rows; j++) {

for (int i=0; i<src.cols; i++) {

ret.at<Vec3f>(j,i)(0) = src.at<Vec3f>(j,i)(0) * T.at<float>(0.0)

+ src.at<Vec3f>(j,i)(1) * T.at<float>(0,1)

+ src.at<Vec3f>(j,i)(2) * T.at<float>(0,2);

ret.at<Vec3f>(j,i)(1) = src.at<Vec3f>(j,i)(0) * T.at<float>(1.0)

+ src.at<Vec3f>(j,i)(1) * T.at<float>(1,1)

+ src.at<Vec3f>(j,i)(2) * T.at<float>(1,2);

ret.at<Vec3f>(j,i)(2) = src.at<Vec3f>(j,i)(0) * T.at<float>(2.0)

+ src.at<Vec3f>(j,i)(1) * T.at<float>(2,1)

+ src.at<Vec3f>(j,i)(2) * T.at<float>(2,2);

}

dst = ret;

}

void ntsc2rgb(const Mat_<Vec3f>& src, Mat_<Vec3f>& dst)

{

Mat ret = src.clone();

Mat T = (Mat_<float>(3,3) << 1.0, 0.956, 0.621, 1.0, -0.272, -0.647, 1.0, -1.106, 1.703);

T = T.t(); //here transpose!

for (int j=0; j<src.rows; j++) {

for (int i=0; i<src.cols; i++) {

ret.at<Vec3f>(j,i)(0) = src.at<Vec3f>(j,i)(0) * T.at<float>(0.0)

+ src.at<Vec3f>(j,i)(1) * T.at<float>(0,1)

+ src.at<Vec3f>(j,i)(2) * T.at<float>(0,2);

ret.at<Vec3f>(j,i)(1) = src.at<Vec3f>(j,i)(0) * T.at<float>(1.0)

+ src.at<Vec3f>(j,i)(1) * T.at<float>(1,1)

+ src.at<Vec3f>(j,i)(2) * T.at<float>(1,2);

ret.at<Vec3f>(j,i)(2) = src.at<Vec3f>(j,i)(0) * T.at<float>(2.0)

+ src.at<Vec3f>(j,i)(1) * T.at<float>(2,1)

+ src.at<Vec3f>(j,i)(2) * T.at<float>(2,2);

}

dst = ret;

}

空間濾波

如前面所述，EVM 算法可以使用拉普拉斯金字塔和高斯金字塔來進(jìn)行空間濾波。使用哪個(gè)金字塔得根據(jù)具體需求而定。如果要放大的是動(dòng)作的變化，那么可以選擇拉普拉斯金字塔，構(gòu)造多個(gè)不同空間頻率的基帶；如果要放大的是顏色的變化，不同基帶的 SNR 應(yīng)該比較接近，因此可以選擇高斯金字塔，只取最頂層下采樣和低通濾波的結(jié)果。這兩個(gè)金字塔可以很容易地利用 OpenCV 的?cv::PyDown()?和?cv::PyUp()?兩個(gè)函數(shù)來構(gòu)造：

/**

* buildLaplacianPyramid - construct a laplacian pyramid from given image

* @param img - source image

* @param levels - levels of the destinate pyramids

* @param pyramid - destinate image

* @return true if success

bool buildLaplacianPyramid(const cv::Mat &img, const int levels,

std::vector<cv::Mat_<cv::Vec3f> > &pyramid)

{

if (levels < 1){

perror("Levels should be larger than 1");

return false;

}

pyramid.clear();

cv::Mat currentImg = img;

for (int l=0; l<levels; l++) {

cv::Mat down,up;

pyrDown(currentImg, down);

pyrUp(down, up, currentImg.size());

cv::Mat lap = currentImg - up;

pyramid.push_back(lap);

currentImg = down;

}

pyramid.push_back(currentImg);

return true;

}

/**

* buildGaussianPyramid - construct a gaussian pyramid from a given image

* @param img - source image

* @param levels - levels of the destinate pyramids

* @param pyramid - destinate image

* @return true if success

bool buildGaussianPyramid(const cv::Mat &img,

const int levels,

std::vector<cv::Mat_<cv::Vec3f> > &pyramid)

{

if (levels < 1){

perror("Levels should be larger than 1");

return false;

}

pyramid.clear();

cv::Mat currentImg = img;

for (int l=0; l<levels; l++) {

cv::Mat down;

cv::pyrDown(currentImg, down);

pyramid.push_back(down);

currentImg = down;

}

return true;

}

時(shí)域?yàn)V波

同樣，時(shí)域?yàn)V波可以根據(jù)不同的需求選擇不同的帶通濾波器。如果需要對(duì)放大結(jié)果進(jìn)行后續(xù)的時(shí)頻分析，則可以選擇理想帶通濾波器；如果不需要對(duì)放大結(jié)果進(jìn)行時(shí)頻分析，可以選擇寬通帶的濾波器，如 Butterworth 帶通濾波器，二級(jí)IIR 濾波器等。這里分別實(shí)現(xiàn)了二階 IIR 帶通濾波器和理想帶通濾波器：

/**

* temporalIIRFilter - temporal IIR filtering an image

* (thanks to Yusuke Tomoto)

* @param pyramid - source image

* @param filtered - filtered result

void VideoProcessor::temporalIIRFilter(const cv::Mat &src,

cv::Mat &dst)

{

cv::Mat temp1 = (1-fh)*lowpass1[curLevel] + fh*src;

cv::Mat temp2 = (1-fl)*lowpass2[curLevel] + fl*src;

lowpass1[curLevel] = temp1;

lowpass2[curLevel] = temp2;

dst = lowpass1[curLevel] - lowpass2[curLevel];

}

/**

* temporalIdalFilter - temporal IIR filtering an image pyramid of concat-frames

* (Thanks to Daniel Ron & Alessandro Gentilini)

* @param pyramid - source pyramid of concatenate frames

* @param filtered - concatenate filtered result

void VideoProcessor::temporalIdealFilter(const cv::Mat &src,

cv::Mat &dst)

{

cv::Mat channels[3];

// split into 3 channels

cv::split(src, channels);

for (int i = 0; i < 3; ++i){

cv::Mat current = channels[i]; // current channel

cv::Mat tempImg;

int width = cv::getOptimalDFTSize(current.cols);

int height = cv::getOptimalDFTSize(current.rows);

cv::copyMakeBorder(current, tempImg,

0, height - current.rows,

0, width - current.cols,

cv::BORDER_CONSTANT, cv::Scalar::all(0));

// do the DFT

cv::dft(tempImg, tempImg, cv::DFT_ROWS | cv::DFT_SCALE, tempImg.rows);

// construct the filter

cv::Mat filter = tempImg.clone();

createIdealBandpassFilter(filter, fl, fh, rate);

// apply filter

cv::mulSpectrums(tempImg, filter, tempImg, cv::DFT_ROWS);

// do the inverse DFT on filtered image

cv::idft(tempImg, tempImg, cv::DFT_ROWS | cv::DFT_SCALE, tempImg.rows);

// copy back to the current channel

tempImg(cv::Rect(0, 0, current.cols, current.rows)).copyTo(channels[i]);

}

// merge channels

cv::merge(channels, 3, dst);

// normalize the filtered image

cv::normalize(dst, dst, 0, 1, CV_MINMAX);

}

/**

* createIdealBandpassFilter - create a 1D ideal band-pass filter

* @param filter - destinate filter

* @param fl - low cut-off

* @param fh - high cut-off

* @param rate - sampling rate(i.e. video frame rate)

void VideoProcessor::createIdealBandpassFilter(cv::Mat &filter, double fl, double fh, double rate)

{

int width = filter.cols;

int height = filter.rows;

fl = 2 * fl * width / rate;

fh = 2 * fh * width / rate;

double response;

for (int i = 0; i < height; ++i) {

for (int j = 0; j < width; ++j) {

// filter response

if (j >= fl && j <= fh)

response = 1.0f;

else

response = 0.0f;

filter.at<float>(i, j) = response;

}

放大變化

根據(jù)前面的公式 5 和公式 8，可以設(shè)計(jì)如下的放大函數(shù)：

/**

* amplify - ampilfy the motion

* @param filtered - motion image

void VideoProcessor::amplify(const cv::Mat &src, cv::Mat &dst)

{

float curAlpha;

switch (spatialType) {

case LAPLACIAN: // for motion magnification

//compute modified alpha for this level

curAlpha = lambda/delta/8 - 1;

curAlpha *= exaggeration_factor;

if (curLevel==levels || curLevel==0) // ignore the highest and lowest frequency band

dst = src * 0;

else

dst = src * cv::min(alpha, curAlpha);

break;

case GAUSSIAN: // for color magnification

dst = src * alpha;

break;

default:

break;

}

對(duì)于動(dòng)作信號(hào)的放大，調(diào)用這個(gè)函數(shù)前需要先算出每一層基帶的?λλ?以及?δ(t)δ(t)?的值，以便于根據(jù)公式 10 來計(jì)算當(dāng)前合理的放大倍數(shù)：

delta = lambda_c/8.0/(1.0+alpha);

// the factor to boost alpha above the bound

// (for better visualization)

exaggeration_factor = 2.0;

// compute the representative wavelength lambda

// for the lowest spatial frequency band of Laplacian pyramid

lambda = sqrt(w*w + h*h)/3; // 3 is experimental constant

注意：這里的?λcλc?就是前面所提的空間波長(zhǎng)的下限值。對(duì)于?λ<λcλ<λc?的基帶，αα?值將被減弱。另外，還加了一個(gè)?exaggeration_factor?參數(shù)，它是一個(gè)魔數(shù)（magic number），用來將符合?λ>λcλ>λc?的基帶加倍放大。

合成圖像

先合成變化信號(hào)的圖像，再與原圖進(jìn)行疊加。根據(jù)使用金字塔的類型，編寫對(duì)應(yīng)的合成方法：

/**

* reconImgFromLaplacianPyramid - reconstruct image from given laplacian pyramid

* @param pyramid - source laplacian pyramid

* @param levels - levels of the pyramid

* @param dst - destinate image

void reconImgFromLaplacianPyramid(const std::vector<cv::Mat_<cv::Vec3f> > &pyramid,

const int levels,

cv::Mat_<cv::Vec3f> &dst)

{

cv::Mat currentImg = pyramid[levels];

for (int l=levels-1; l>=0; l--) {

cv::Mat up;

cv::pyrUp(currentImg, up, pyramid[l].size());

currentImg = up + pyramid[l];

}

dst = currentImg.clone();

}

/**

* upsamplingFromGaussianPyramid - up-sampling an image from gaussian pyramid

* @param src - source image

* @param levels - levels of the pyramid

* @param dst - destinate image

void upsamplingFromGaussianPyramid(const cv::Mat &src,

const int levels,

cv::Mat_<cv::Vec3f> &dst)

{

cv::Mat currentLevel = src.clone();

for (int i = 0; i < levels; ++i) {

cv::Mat up;

cv::pyrUp(currentLevel, up);

currentLevel = up;

}

currentLevel.copyTo(dst);

}

衰減 I、Q 通道

對(duì)于動(dòng)作信號(hào)的放大，可以在后期引入一個(gè)衰減因子減弱 I、Q 兩個(gè)通道的變化幅度，最后才轉(zhuǎn)回 RGB 顏色空間。

/**

* attenuate - attenuate I, Q channels

* @param src - source image

* @param dst - destinate image

void VideoProcessor::attenuate(cv::Mat &src, cv::Mat &dst)

{

cv::Mat planes[3];

cv::split(src, planes);

planes[1] = planes[1] * chromAttenuation;

planes[2] = planes[2] * chromAttenuation;

cv::merge(planes, 3, dst);

}

結(jié)果

下面演示使用我的程序，對(duì)論文提供的 face 案例進(jìn)行處理的結(jié)果。

所選參數(shù)如下：

?ααλcλc頻段衰減因子

動(dòng)作變化放大	10	80	0.05~0.4	0.1
顏色變化放大	50	-	0.83~1.0	-

源碼和程序

?QtEVM 的源碼
?QtEVM 的 Win32 可執(zhí)行程序

總結(jié)

以上是生活随笔為你收集整理的Eulerian Video Magnification的全部?jī)?nèi)容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網(wǎng)站內(nèi)容還不錯(cuò)，歡迎將生活随笔推薦給好友。

上一篇： openCV——图像金字塔
下一篇： CAD可以用数位板吗？矢量图有必要用数位