图像水波纹特效原理分析和实现

 　　前段时间注意到一些软件上有图像的水波纹特效，似乎很炫，想深入了解下该效果的具体原理与实现方式，上网搜了不少些资料，都讲得不清不楚，没办法只能靠自己了。花了一整个下午先去复习了高中物理的波的知识，试着自己来推导原理并实现了下。下面的推导是我根据一些资料以及自己分析出的，如有错误，望请指出。上张效果图先：

基本原理 

---

　　水波效果反映到图像上，则是像素点的偏移。因此对图像的处理就如同图像缩放一样，对于输出图像的每个点，计算其对应于原始输入图像的像素点，通过插值的方法即可计算其颜色值。 

　　先说下关于水波的基础知识。波的传播方向与质点振动方向垂直的为横波，相同则为纵波，水波是横波和纵波的叠加，因此水波在传播时，每个质点存在水平运动和上下运动，合起来就是一个椭圆的圆周运动，如下图的蓝色椭圆。 

　　如图，红轴表示水面，对于图像来讲，只有在垂直于成像视线的方向（图中蓝色线段）产生的运动才会引起像素的偏移，而平行于实现方向则可以忽略。那么对于某个点x，其运动轨迹上的每个点都可以分解为与视线平行和垂直的2个方向上的位移，我们只要计算出垂直视线的位移y即可计算出水波导致像素点的偏移位移。 

 　　假定水波是平面简谐波，则能得到其波动方程

 　　式中，A是振幅，v是波速，t是时间，lamda是波长， x为距离波源的距离，phi为初始相位。但是随着时间的推移和波的传播，其能量会衰减，即其振幅会随着时间和位置衰减，则上面的公式需修改为：

 　　函数表征了水波能量的衰减情况，因此该函数必须是关于x和t的单调递减函数，实际水波的递减是什么规律我不清楚，但这里自行设定一个即可，一般的衰减函数有线性衰减和指数衰减，具体好坏只要实验结果才具有说服力。 

　　一个波源可以由波长，波速，初始相位和振幅唯一确定，因此对于某个时刻的某个点，我们可以计算出其位移。另一方面，水波是向四面八分同时传播的，对图像坐标系来看，每个点的位移都可以分解为水平与垂直位移，从而就可以确定像素点的精确偏移了。当然，对于多个波源，可以分别考虑，每个点的位移是所有波源引起的位移和。

水波纹特效具体实现 

---

 下面就是具体实现了，这个效果实现应该比较简单的。首先有一个波源对象。

typedef unsigned char BYTE;
#define PI 3.1415926
//一个波源
class CWaveSource
{
public:
    CWaveSource(float f32Cx = 100.0f, float f32Cy = 100.0f, 
                float f32MaxAmp = 20.0f, float f32Phase = -PI/2,
                float f32WaveLen = 10.0f, float f32V = 1.0f);
    ~CWaveSource();

    //该波源在当前时刻，x，y位置处的振幅
    void GetAmplitude(int x, int y, float &f32ax, float &f32ay);

    //波传播到下一时刻，需修改相应的波的状态
    bool Propagate();

private:
    float m_f32Cx;        //波源中心点坐标
    float m_f32Cy;
    float m_f32MaxAmp;    //t时刻中心点最大振幅，t=0最大
    float m_f32InitPhase; //初始相位
    float m_f32WaveLen;   //波长
    float m_f32Velocity;  //波速
    int   m_nTime;        //时间

    bool  m_bDisappear;   //该波源是否可以消失，随着时间变化，其能量衰减到很小的时候可以认为其消失了
    float m_f32CurRadius; //当前时刻波的半径,用于节约计算量的变量
};

这里面注释比较清楚，就不多说了。该对象对应的实现为：

CWaveSource::CWaveSource(float f32Cx, float f32Cy, 
                         float f32MaxAmp, float f32Phase, 
                         float f32WaveLen, float f32V)
{
    m_f32Cx = f32Cx;
    m_f32Cy = f32Cy;
    m_f32MaxAmp = f32MaxAmp;
    m_f32InitPhase = f32Phase;
    m_f32WaveLen = f32WaveLen;
    m_f32Velocity = f32V;
    m_nTime = 0;
    m_bDisappear = false;
    m_f32CurRadius = 1e-12;
}
CWaveSource::~CWaveSource()
{

}

bool CWaveSource::Propagate()
{
    //下一个时刻
    m_nTime++;
    m_f32CurRadius += m_f32Velocity;
    //考虑最大振幅随时间的衰减,这里的衰减函数可以自己设定
    //只要振幅随着时间递减即可，最好能做到比较自然
    m_f32MaxAmp /= 1.01f;

    //能量衰减到一定程度后，可以认为该波源已消失，通过返回让外面把它删掉
    if(m_f32MaxAmp < 1e-3)
    {
        m_bDisappear = true;
    }
    return m_bDisappear;
}

void CWaveSource::GetAmplitude(int x, int y, float &f32ax, float &f32ay)
{
    f32ax = f32ay = 0.0f;
    //波源已消失，防止上层未删掉，耗计算量
    if(m_bDisappear) return;
    //注：这里可以建一个表，图像区域每个点相对波源中心点位置是固定的，即是f32Dist固定
    //那么位置对应的衰减系数就应该是固定的
    //这样对于一个波源只需计算一次f32Dist和衰减系数r了
    //此处为了代码的可读性，暂时不那样实现    
    float f32Radius = m_f32CurRadius;
    float f32Dx = (m_f32Cx - x);
    float f32Dy = (m_f32Cy - y);

    float f32Dist = f32Dx * f32Dx + f32Dy * f32Dy;
    if(f32Dist > f32Radius * f32Radius) return;

    f32Dist = sqrt(f32Dist);
    //考虑当前点最大振幅随位置的衰减,这里的衰减函数也可以自己设定
    float r = 1 - f32Dist / 1000.0f;//exp(-f32Dist/10);
    float f32MaxAmp = m_f32MaxAmp * r;

    //计算当前点的振幅,每个点的相位在不停变化的 y = Acos(2*pi*(vt-x)/wavelen + phase)
    float f32Temp = 2 * PI * (m_f32Velocity * m_nTime - f32Dist) / m_f32WaveLen + m_f32InitPhase;
    float f32CurAmp = f32MaxAmp * sin(f32Temp);

    f32ax = f32CurAmp * ABS(f32Dx)/f32Dist;
    f32ax = f32CurAmp * ABS(f32Dy)/f32Dist;
}

好了，这样一个波源就已经构造好了。对图像的水波特效处理对象如下：

class CRippling
{
public:
    CRippling();
    ~CRippling();
    
    //增加一个波源，该函数最好能重载一个直接输波源参数的，这里暂时忽略
    void AddSource(CWaveSource &WaveSource);
    
    //当前已有的波源对图像进行处理
    void Process(BYTE *pu8Dst, BYTE *pu8Src, int nWidth, int nHeight, int nChannels = 1);
    
    //所有波源传播
    void Propagate();

private:
    //获取当前时刻所有波源作用下，像素点x，y对应于原始图像的像素位置，返回到x，y中
    void GetPos(float &x, float &y);

    //获取图像亚像素值，采用双线性插值
    void GetSubPixel(BYTE *pu8Dst, BYTE *pu8Src, int nWidth, int nHeight, int nChannels, float x, float y);

    //波源相关信息,因为需要插入和删除故采用list管理，也可以使用数组或者其他数据结构
    list<CWaveSource> m_lWaveSrc;
};

CRippling::CRippling()
{
}

CRippling::~CRippling()
{
}

void CRippling::AddSource(CWaveSource &WaveSource)
{
    m_lWaveSrc.push_back(WaveSource);
}

void CRippling::Process(BYTE *pu8Dst, BYTE *pu8Src, int nWidth, int nHeight, int nChannels)
{
    int nRow, nCol;
    float x,y;
    for(nRow = 0; nRow < nHeight; nRow++)
    {
        for(nCol = 0; nCol < nWidth; nCol++)
        {
            x = nCol;
            y = nRow;
            GetPos(x, y);
            GetSubPixel(pu8Dst, pu8Src, nWidth, nHeight, nChannels, x, y);
            pu8Dst += nChannels;
        }
    }
}

void CRippling::Propagate()
{
    //每个波源都需要传播
    list<CWaveSource>::iterator It = m_lWaveSrc.begin();
    while(It != m_lWaveSrc.end())
    {
        bool bDisappear = It->Propagate();
        //波源可以消失，则删掉
        if(bDisappear) It = m_lWaveSrc.erase(It); 
        else It++;
    }
}

void CRippling::GetPos(float &x, float &y)
{
    //根据平面简谐波的特性，某个点的位移是所有波位移和
    float f32AmpX = 0;
    float f32AmpY = 0;
    list<CWaveSource>::iterator It = m_lWaveSrc.begin();
    while(It != m_lWaveSrc.end())
    {
        It->GetAmplitude(x, y, f32AmpX, f32AmpY);
        x += f32AmpX;
        y += f32AmpY;
        It++;
    }
}

void CRippling::GetSubPixel(BYTE *pu8Dst, BYTE *pu8Src, int nWidth, int nHeight, int nChannels, float x, float y)
{
    //采用双线性插值，就实际效果来看，最近邻插值也可以的
    int x0 = x;
    int y0 = y;
    int x1 = x0 + 1;
    int y1 = y0 + 1;

    float wx1 = x - x0;
    float wy1 = y - y0;
    float wx0 = 1.0f - wx1;
    float wy0 = 1.0f - wy1;

    x0 = MAX(x0, 0);
    x1 = MAX(x1, 0);
    y0 = MAX(y0, 0);
    y1 = MAX(y1, 0);

    x0 = MIN(x0, nWidth-1);
    x1 = MIN(x1, nWidth-1);
    y0 = MIN(y0, nHeight-1);
    y1 = MIN(y1, nHeight-1);
    for(int i = 0; i < nChannels; i++)
    {
        u8 * pu8Y0 = pu8Src + y0 * nWidth * nChannels;
        u8 * pu8Y1 = pu8Src + y1 * nWidth * nChannels;

        float sum_y0 = (pu8Y0[x0 * nChannels + i] * wx0 + pu8Y0[x1 * nChannels + i] * wx1);
        float sum_y1 = (pu8Y1[x0 * nChannels + i] * wx0 + pu8Y1[x1 * nChannels + i] * wx1);
        pu8Dst[i] = (BYTE)(sum_y0 * wy0 + sum_y1 * wy1);
    }
}

测试代码采用隔一定时间增加一个波源，可以产生动态效果。其中调用了opencv读写和显示图像，相关的头文件和库自己加上即可。

void test_rippling()
{    
    Mat mImg = imread("C:/test2.jpg");
    imshow("org", mImg);
    CRippling rippling;
    for(int t = 0; t < 10000; t++)
    {
        if(t % 20 == 0)
        {
            CWaveSource wavesource(50 + rand() % 800, 50 + rand() % 300);
            rippling.AddSource(wavesource);
        }

        printf("
times: %d", t);
        Mat mSrc = mImg.clone();
        Mat mDst = mImg.clone();

        rippling.Process(mDst.data, mSrc.data, mSrc.cols, mSrc.rows, 3);
        rippling.Propagate();
        imshow("pro", mDst);
        waitKey(1);
    }
    waitKey(0);
}

这样基本就完成了，可以实现最开始那张图的效果了。

总结 

---

　　上面的分析与实现只是按照原始的产生与传播过程来做的，可以想象，随着波源数的增加，算法的耗时会不断增加，当然在代码注释中也在有的地方说明了可以优化的点。另外，如果我们仅仅是实现这样一个效果，很多波源参数给固定下来，那么可以在算法上进一步推导，从而简化运算。比如，固定波速与波长，让波的周期是4个或2个单位之间，那么应该可以推导出后一时刻位移为前面几个时刻位移的关系，这样就不用每次计算sin函数了。

　　当然上面的推导只考虑了简单情况，特别地，成像一般会近大远小，因此在图像上下方波的扩张速度和质点偏移都是不一样的，对于完整的模型这些都是应该要考虑的因素。

---