圖像水波紋特效原理分析和實現

 　　前段時間注意到一些軟件上有圖像的水波紋特效，彷佛很炫，想深刻了解下該效果的具體原理與實現方式，上網搜了很多些資料，都講得不清不楚，沒辦法只能靠本身了。花了一整個下午先去複習了高中物理的波的知識，試着本身來推導原理並實現了下。下面的推導是我根據一些資料以及本身分析出的，若有錯誤，望請指出。上張效果圖先：

 

基本原理 

---

　　水波效果反映到圖像上，則是像素點的偏移。所以對圖像的處理就如同圖像縮放同樣，對於輸出圖像的每一個點，計算其對應於原始輸入圖像的像素點，經過插值的方法便可計算其顏色值。 

　　先說下關於水波的基礎知識。波的傳播方向與質點振動方向垂直的爲橫波，相同則爲縱波，水波是橫波和縱波的疊加，所以水波在傳播時，每一個質點存在水平運動和上下運動，合起來就是一個橢圓的圓周運動，以下圖的藍色橢圓。 

　　如圖，紅軸表示水面，對於圖像來說，只有在垂直於成像視線的方向（圖中藍色線段）產生的運動纔會引發像素的偏移，而平行於實現方向則能夠忽略。那麼對於某個點x，其運動軌跡上的每一個點均可以分解爲與視線平行和垂直的2個方向上的位移，咱們只要計算出垂直視線的位移y便可計算出水波致使像素點的偏移位移。 

 　　假定水波是平面簡諧波，則能獲得其波動方程

 　　式中，A是振幅，v是波速，t是時間，lamda是波長， x爲距離波源的距離，phi爲初始相位。可是隨着時間的推移和波的傳播，其能量會衰減，即其振幅會隨着時間和位置衰減，則上面的公式需修改成：

 　　函數表徵了水波能量的衰減狀況，所以該函數必須是關於x和t的單調遞減函數，實際水波的遞減是什麼規律我不清楚，但這裏自行設定一個便可，通常的衰減函數有線性衰減和指數衰減，具體好壞只要實驗結果才具備說服力。 

　　一個波源能夠由波長，波速，初始相位和振幅惟一肯定，所以對於某個時刻的某個點，咱們能夠計算出其位移。另外一方面，水波是向四面八分同時傳播的，對圖像座標系來看，每一個點的位移均可以分解爲水平與垂直位移，從而就能夠肯定像素點的精確偏移了。固然，對於多個波源，能夠分別考慮，每一個點的位移是全部波源引發的位移和。

 

水波紋特效具體實現 

---

 下面就是具體實現了，這個效果實現應該比較簡單的。首先有一個波源對象。

 1 typedef unsigned char BYTE;
 2 #define PI 3.1415926
 3 //一個波源
 4 class CWaveSource
 5 {
 6 public:
 7     CWaveSource(float f32Cx = 100.0f, float f32Cy = 100.0f, 
 8                 float f32MaxAmp = 20.0f, float f32Phase = -PI/2,
 9                 float f32WaveLen = 10.0f, float f32V = 1.0f);
10     ~CWaveSource();
11 
12     //該波源在當前時刻，x，y位置處的振幅
13     void GetAmplitude(int x, int y, float &f32ax, float &f32ay);
14 
15     //波傳播到下一時刻，需修改相應的波的狀態
16     bool Propagate();
17 
18 private:
19     float m_f32Cx;        //波源中心點座標
20     float m_f32Cy;
21     float m_f32MaxAmp;    //t時刻中心點最大振幅，t=0最大
22     float m_f32InitPhase; //初始相位
23     float m_f32WaveLen;   //波長
24     float m_f32Velocity;  //波速
25     int   m_nTime;        //時間
26 
27     bool  m_bDisappear;   //該波源是否能夠消失，隨着時間變化，其能量衰減到很小的時候能夠認爲其消失了
28     float m_f32CurRadius; //當前時刻波的半徑,用於節約計算量的變量
29 };

CWaveSource

這裏面註釋比較清楚，就很少說了。該對象對應的實現爲：

 1 CWaveSource::CWaveSource(float f32Cx, float f32Cy, 
 2                          float f32MaxAmp, float f32Phase, 
 3                          float f32WaveLen, float f32V)
 4 {
 5     m_f32Cx = f32Cx;
 6     m_f32Cy = f32Cy;
 7     m_f32MaxAmp = f32MaxAmp;
 8     m_f32InitPhase = f32Phase;
 9     m_f32WaveLen = f32WaveLen;
10     m_f32Velocity = f32V;
11     m_nTime = 0;
12     m_bDisappear = false;
13     m_f32CurRadius = 1e-12;
14 }
15 CWaveSource::~CWaveSource()
16 {
17 
18 }
19 
20 bool CWaveSource::Propagate()
21 {
22     //下一個時刻
23     m_nTime++;
24     m_f32CurRadius += m_f32Velocity;
25     //考慮最大振幅隨時間的衰減,這裏的衰減函數能夠本身設定
26     //只要振幅隨着時間遞減便可，最好能作到比較天然
27     m_f32MaxAmp /= 1.01f;
28 
29     //能量衰減到必定程度後，能夠認爲該波源已消失，經過返回讓外面把它刪掉
30     if(m_f32MaxAmp < 1e-3)
31     {
32         m_bDisappear = true;
33     }
34     return m_bDisappear;
35 }
36 
37 void CWaveSource::GetAmplitude(int x, int y, float &f32ax, float &f32ay)
38 {
39     f32ax = f32ay = 0.0f;
40     //波源已消失，防止上層未刪掉，耗計算量
41     if(m_bDisappear) return;
42     //注：這裏能夠建一個表，圖像區域每一個點相對波源中心點位置是固定的，便是f32Dist固定
43     //那麼位置對應的衰減係數就應該是固定的
44     //這樣對於一個波源只需計算一次f32Dist和衰減係數r了
45     //此處爲了代碼的可讀性，暫時不那樣實現    
46     float f32Radius = m_f32CurRadius;
47     float f32Dx = (m_f32Cx - x);
48     float f32Dy = (m_f32Cy - y);
49 
50     float f32Dist = f32Dx * f32Dx + f32Dy * f32Dy;
51     if(f32Dist > f32Radius * f32Radius) return;
52 
53     f32Dist = sqrt(f32Dist);
54     //考慮當前點最大振幅隨位置的衰減,這裏的衰減函數也能夠本身設定
55     float r = 1 - f32Dist / 1000.0f;//exp(-f32Dist/10);
56     float f32MaxAmp = m_f32MaxAmp * r;
57 
58     //計算當前點的振幅,每一個點的相位在不停變化的 y = Acos(2*pi*(vt-x)/wavelen + phase)
59     float f32Temp = 2 * PI * (m_f32Velocity * m_nTime - f32Dist) / m_f32WaveLen + m_f32InitPhase;
60     float f32CurAmp = f32MaxAmp * sin(f32Temp);
61 
62     f32ax = f32CurAmp * ABS(f32Dx)/f32Dist;
63     f32ax = f32CurAmp * ABS(f32Dy)/f32Dist;
64 }

View Code

好了，這樣一個波源就已經構造好了。對圖像的水波特效處理對象以下：

  1 class CRippling
  2 {
  3 public:
  4     CRippling();
  5     ~CRippling();
  6     
  7     //增長一個波源，該函數最好能重載一個直接輸波源參數的，這裏暫時忽略
  8     void AddSource(CWaveSource &WaveSource);
  9     
 10     //當前已有的波源對圖像進行處理
 11     void Process(BYTE *pu8Dst, BYTE *pu8Src, int nWidth, int nHeight, int nChannels = 1);
 12     
 13     //全部波源傳播
 14     void Propagate();
 15 
 16 private:
 17     //獲取當前時刻全部波源做用下，像素點x，y對應於原始圖像的像素位置，返回到x，y中
 18     void GetPos(float &x, float &y);
 19 
 20     //獲取圖像亞像素值，採用雙線性插值
 21     void GetSubPixel(BYTE *pu8Dst, BYTE *pu8Src, int nWidth, int nHeight, int nChannels, float x, float y);
 22 
 23     //波源相關信息,由於須要插入和刪除故採用list管理，也可使用數組或者其餘數據結構
 24     list<CWaveSource> m_lWaveSrc;
 25 };
 26 
 27 CRippling::CRippling()
 28 {
 29 }
 30 
 31 CRippling::~CRippling()
 32 {
 33 }
 34 
 35 void CRippling::AddSource(CWaveSource &WaveSource)
 36 {
 37     m_lWaveSrc.push_back(WaveSource);
 38 }
 39 
 40 void CRippling::Process(BYTE *pu8Dst, BYTE *pu8Src, int nWidth, int nHeight, int nChannels)
 41 {
 42     int nRow, nCol;
 43     float x,y;
 44     for(nRow = 0; nRow < nHeight; nRow++)
 45     {
 46         for(nCol = 0; nCol < nWidth; nCol++)
 47         {
 48             x = nCol;
 49             y = nRow;
 50             GetPos(x, y);
 51             GetSubPixel(pu8Dst, pu8Src, nWidth, nHeight, nChannels, x, y);
 52             pu8Dst += nChannels;
 53         }
 54     }
 55 }
 56 
 57 void CRippling::Propagate()
 58 {
 59     //每一個波源都須要傳播
 60     list<CWaveSource>::iterator It = m_lWaveSrc.begin();
 61     while(It != m_lWaveSrc.end())
 62     {
 63         bool bDisappear = It->Propagate();
 64         //波源能夠消失，則刪掉
 65         if(bDisappear) It = m_lWaveSrc.erase(It); 
 66         else It++;
 67     }
 68 }
 69 
 70 void CRippling::GetPos(float &x, float &y)
 71 {
 72     //根據平面簡諧波的特性，某個點的位移是全部波位移和
 73     float f32AmpX = 0;
 74     float f32AmpY = 0;
 75     list<CWaveSource>::iterator It = m_lWaveSrc.begin();
 76     while(It != m_lWaveSrc.end())
 77     {
 78         It->GetAmplitude(x, y, f32AmpX, f32AmpY);
 79         x += f32AmpX;
 80         y += f32AmpY;
 81         It++;
 82     }
 83 }
 84 
 85 void CRippling::GetSubPixel(BYTE *pu8Dst, BYTE *pu8Src, int nWidth, int nHeight, int nChannels, float x, float y)
 86 {
 87     //採用雙線性插值，就實際效果來看，最近鄰插值也能夠的
 88     int x0 = x;
 89     int y0 = y;
 90     int x1 = x0 + 1;
 91     int y1 = y0 + 1;
 92 
 93     float wx1 = x - x0;
 94     float wy1 = y - y0;
 95     float wx0 = 1.0f - wx1;
 96     float wy0 = 1.0f - wy1;
 97 
 98     x0 = MAX(x0, 0);
 99     x1 = MAX(x1, 0);
100     y0 = MAX(y0, 0);
101     y1 = MAX(y1, 0);
102 
103     x0 = MIN(x0, nWidth-1);
104     x1 = MIN(x1, nWidth-1);
105     y0 = MIN(y0, nHeight-1);
106     y1 = MIN(y1, nHeight-1);
107     for(int i = 0; i < nChannels; i++)
108     {
109         u8 * pu8Y0 = pu8Src + y0 * nWidth * nChannels;
110         u8 * pu8Y1 = pu8Src + y1 * nWidth * nChannels;
111 
112         float sum_y0 = (pu8Y0[x0 * nChannels + i] * wx0 + pu8Y0[x1 * nChannels + i] * wx1);
113         float sum_y1 = (pu8Y1[x0 * nChannels + i] * wx0 + pu8Y1[x1 * nChannels + i] * wx1);
114         pu8Dst[i] = (BYTE)(sum_y0 * wy0 + sum_y1 * wy1);
115     }
116 }

View Code

測試代碼採用隔必定時間增長一個波源，能夠產生動態效果。其中調用了opencv讀寫和顯示圖像，相關的頭文件和庫本身加上便可。

 1 void test_rippling()
 2 {    
 3     Mat mImg = imread("C:/test2.jpg");
 4     imshow("org", mImg);
 5     CRippling rippling;
 6     for(int t = 0; t < 10000; t++)
 7     {
 8         if(t % 20 == 0)
 9         {
10             CWaveSource wavesource(50 + rand() % 800, 50 + rand() % 300);
11             rippling.AddSource(wavesource);
12         }
13 
14         printf("\rtimes: %d", t);
15         Mat mSrc = mImg.clone();
16         Mat mDst = mImg.clone();
17 
18         rippling.Process(mDst.data, mSrc.data, mSrc.cols, mSrc.rows, 3);
19         rippling.Propagate();
20         imshow("pro", mDst);
21         waitKey(1);
22     }
23     waitKey(0);
24 }

View Code

這樣基本就完成了，能夠實現最開始那張圖的效果了。

 

總結 

---

　　上面的分析與實現只是按照原始的產生與傳播過程來作的，能夠想象，隨着波源數的增長，算法的耗時會不斷增長，固然在代碼註釋中也在有的地方說明了能夠優化的點。另外，若是咱們僅僅是實現這樣一個效果，不少波源參數給固定下來，那麼能夠在算法上進一步推導，從而簡化運算。好比，固定波速與波長，讓波的週期是4個或2個單位之間，那麼應該能夠推導出後一時刻位移爲前面幾個時刻位移的關係，這樣就不用每次計算sin函數了。

　　固然上面的推導只考慮了簡單狀況，特別地，成像通常會近大遠小，所以在圖像上下方波的擴張速度和質點偏移都是不同的，對於完整的模型這些都是應該要考慮的因素。