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

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

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

为了生成具有一定特征凸起的平滑表面，利用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:前两个边界条件是常见泊松方程边界条件，最后一个是为了满足特殊需求新定义的边界条件）

  1 #version 400
  2 #extension GL_ARB_texture_rectangle : enable
  3 
  4 
  5 //original image which will not be changed
  6 uniform sampler2DRect inputImg;
  7 
  8 //red as dirichlet boundary,blue as neumann boundary
  9 uniform sampler2DRect constraintImg;
 10 
 11 uniform sampler2DRect tmpImg;
 12 
 13 uniform vec2 size;
 14 
 15 
 16 void main()
 17 {
 18     vec2 coord = gl_FragCoord.xy;
 19     vec4 inputColor=texture2DRect(inputImg, coord);
 20     vec4 tmpColor = texture2DRect(tmpImg, coord);
 21     vec4 constraintColor= texture2DRect(constraintImg, coord);
 22     
 23     if(inputColor.w<0.01)//input 4channel as mask
 24     {
 25         //outside mask,should be zero
 26         gl_FragColor = vec4(0,0,0,1);
 27     }
 28     else if(constraintColor.x>0.99)//red
 29     {
 30         //dirichlet boundary,remain the original color
 31         gl_FragColor = inputColor;
 32     }
 33     else if(constraintColor.x<0.01 && constraintColor.y>0.01)//green
 34     {
 35         if(constraintColor.y>0.99)
 36             gl_FragColor = inputColor;
 37         else
 38         {
 39             vec4 f = vec4(constraintColor.w,constraintColor.w,constraintColor.w,constraintColor.w);
 40             vec4 a = texture2DRect(tmpImg, coord + vec2(-1, 0));
 41             vec4 b = texture2DRect(tmpImg, coord + vec2(1, 0));
 42             vec4 c = texture2DRect(tmpImg, coord + vec2(0, -1));
 43             vec4 d = texture2DRect(tmpImg, coord + vec2(0, 1));
 44             vec4 res= 0.25*(a+b+c+d+f);
 45             res=(1 - constraintColor.y) * res + constraintColor.y * inputColor;
 46             gl_FragColor = res;
 47         }
 48     }
 49     else if(constraintColor.x<0.01 && constraintColor.z>0.99)//blue
 50     {
 51         //neumann boundary,derivative is zero,should be inside mask
 52         //so we decide to use the average color of its inside region neighbors as the approximation
 53         int count=0;
 54         vec4 res = vec4(0,0,0,0);
 55         vec4 ma = texture2DRect(inputImg, coord+vec2(-1,0));
 56         if(ma.w>=0.9)//inside mask region
 57         {
 58             vec4 a = texture2DRect(tmpImg, coord+vec2(-1,0));
 59             res+=a;
 60             count+=1;
 61             
 62         }
 63 
 64         vec4 mb = texture2DRect(inputImg, coord+vec2(1,0));
 65         if(mb.w>=0.9)
 66         {
 67             vec4 b = texture2DRect(tmpImg, coord + vec2(1, 0));
 68             res+=b;
 69             count+=1;
 70         }
 71 
 72         vec4 mc = texture2DRect(inputImg, coord+vec2(0,-1));
 73         if(mc.w>=0.9)
 74         {
 75             vec4 c = texture2DRect(tmpImg, coord + vec2(0, -1));
 76             res+=c;
 77             count+=1;
 78         }
 79 
 80         vec4 md = texture2DRect(inputImg, coord+vec2(0,1));
 81         if(md.w>=0.9)
 82         {        
 83             vec4 d = texture2DRect(tmpImg, coord + vec2(0, 1));
 84             res+=d;
 85             count+=1;
 86         }
 87         
 88         if(count>0)
 89             res/=count;
 90 
 91         gl_FragColor=res;
 92     }
 93     else
 94     {
 95         //inside poisson,jacobi iteration
 96         vec4 f = vec4(constraintColor.w,constraintColor.w,constraintColor.w,constraintColor.w);
 97         vec4 a = texture2DRect(tmpImg, coord + vec2(-1, 0));
 98         vec4 b = texture2DRect(tmpImg, coord + vec2(1, 0));
 99         vec4 c = texture2DRect(tmpImg, coord + vec2(0, -1));
100         vec4 d = texture2DRect(tmpImg, coord + vec2(0, 1));
101         vec4 res = 0.25*(a+b+c+d+f);
102 
103         gl_FragColor=res;
104     }
105     
106 }