利用着色器迭代求解离散泊松方程

一直想整理一下研究生期间做的东西，在这里做一个开端吧。从毕业论文中用到的算法开始。

求解离散泊松方程，源于当时的研究项目，从图像生成三维模型，具体内容此处不做细谈。

为了生成具有一定特征凸起的平滑表面，利用G2连续性曲面的二阶偏微分具有G0连续性的特性，构造特征参数表面，该表面只需G0连续性，通过距离变换等计算，可以很容易获得。通过迭代计算，即可获得所需的具有G2连续性的表面。

构造泊松等式：

d²z(x,y)/dx² +d²z(x,y)/dy² =f(x,y)

其中xy为自变量，z为所求的表面高度，f(x,y)为z对于x，y的偏微分参数，即特征表面（保证G0连续）。

将等式写成有限差分形式：

d²z(i,j)/dx² +d²z(i,j)/dy² = z(i-1,j)+z(i,j-1)+z(i+1,j)+z(i,j+1)-4z(i,j)

利用等式的右侧可以得到：

z(i,j)=1/4(z(i-1,j)+z(i,j-1)+z(i+1,j)+z(i,j+1)+f(i,j))

构造Jacobi迭代:

zⁿ⁺¹(i,j)=1/4(zⁿ(i-1,j)+zⁿ(i,j-1)+zⁿ(i+1,j)+zⁿ(i,j+1)+f(i,j))

迭代计算z(i,j)，直到最终收敛。

这样一个等式的计算，可以很容易的写出来，但是对于1000*1000的点阵可能就会耗时数小时。

为了加快计算速度，采用了片段着色器绘制来模拟计算过程,让输出的颜色值表示高度，计算过程可以加快数十倍。

使用两张FrameBufferTexture，轮换作为输入和输出来模拟迭代计算的过程，可以非常迅速的求解。

由于整个工程浩大，这里仅贴上GLSL部分的核心计算代码。其中constraintImg是边界条件对应的图像，红色代表dirichlet边界(边界处z值保持不变)，蓝色代表neuman边界(边界处z值的梯度为0)，绿色代表边界处的高度随输入渐变。（ps:前两个边界条件是常见泊松方程边界条件，最后一个是为了满足特殊需求新定义的边界条件）

#version 400
#extension GL_ARB_texture_rectangle : enable


//original image which will not be changed
uniform sampler2DRect inputImg;

//red as dirichlet boundary,blue as neumann boundary
uniform sampler2DRect constraintImg;

uniform sampler2DRect tmpImg;

uniform vec2 size;


void main()
{
    vec2 coord = gl_FragCoord.xy;
    vec4 inputColor=texture2DRect(inputImg, coord);
    vec4 tmpColor = texture2DRect(tmpImg, coord);
    vec4 constraintColor= texture2DRect(constraintImg, coord);
    
    if(inputColor.w<0.01)//input 4channel as mask
    {
        //outside mask,should be zero
        gl_FragColor = vec4(0,0,0,1);
    }
    else if(constraintColor.x>0.99)//red
    {
        //dirichlet boundary,remain the original color
        gl_FragColor = inputColor;
    }
    else if(constraintColor.x<0.01 && constraintColor.y>0.01)//green
    {
        if(constraintColor.y>0.99)
            gl_FragColor = inputColor;
        else
        {
            vec4 f = vec4(constraintColor.w,constraintColor.w,constraintColor.w,constraintColor.w);
            vec4 a = texture2DRect(tmpImg, coord + vec2(-1, 0));
            vec4 b = texture2DRect(tmpImg, coord + vec2(1, 0));
            vec4 c = texture2DRect(tmpImg, coord + vec2(0, -1));
            vec4 d = texture2DRect(tmpImg, coord + vec2(0, 1));
            vec4 res= 0.25*(a+b+c+d+f);
            res=(1 - constraintColor.y) * res + constraintColor.y * inputColor;
            gl_FragColor = res;
        }
    }
    else if(constraintColor.x<0.01 && constraintColor.z>0.99)//blue
    {
        //neumann boundary,derivative is zero,should be inside mask
        //so we decide to use the average color of its inside region neighbors as the approximation
        int count=0;
        vec4 res = vec4(0,0,0,0);
        vec4 ma = texture2DRect(inputImg, coord+vec2(-1,0));
        if(ma.w>=0.9)//inside mask region
        {
            vec4 a = texture2DRect(tmpImg, coord+vec2(-1,0));
            res+=a;
            count+=1;
            
        }

        vec4 mb = texture2DRect(inputImg, coord+vec2(1,0));
        if(mb.w>=0.9)
        {
            vec4 b = texture2DRect(tmpImg, coord + vec2(1, 0));
            res+=b;
            count+=1;
        }

        vec4 mc = texture2DRect(inputImg, coord+vec2(0,-1));
        if(mc.w>=0.9)
        {
            vec4 c = texture2DRect(tmpImg, coord + vec2(0, -1));
            res+=c;
            count+=1;
        }

        vec4 md = texture2DRect(inputImg, coord+vec2(0,1));
        if(md.w>=0.9)
        {        
            vec4 d = texture2DRect(tmpImg, coord + vec2(0, 1));
            res+=d;
            count+=1;
        }
        
        if(count>0)
            res/=count;

        gl_FragColor=res;
    }
    else
    {
        //inside poisson,jacobi iteration
        vec4 f = vec4(constraintColor.w,constraintColor.w,constraintColor.w,constraintColor.w);
        vec4 a = texture2DRect(tmpImg, coord + vec2(-1, 0));
        vec4 b = texture2DRect(tmpImg, coord + vec2(1, 0));
        vec4 c = texture2DRect(tmpImg, coord + vec2(0, -1));
        vec4 d = texture2DRect(tmpImg, coord + vec2(0, 1));
        vec4 res = 0.25*(a+b+c+d+f);

        gl_FragColor=res;
    }
    
}