利用着色器迭代求解离散泊松方程
一直想整理一下研究生期间做的东西，在这里做一个开端吧。从毕业论文中用到的算法开始。
求解离散泊松方程，源于当时的研究项目，从图像生成三维模型，具体内容此处不做细谈。
为了生成具有一定特征凸起的平滑表面，利用G2连续性曲面的二阶偏微分具有G0连续性的特性，构造特征参数表面，该表面只需G0连续性，通过距离变换等计算，可以很容易获得。通过迭代计算，即可获得所需的具有G2连续性的表面。
构造泊松等式：
d2z(x,y)/dx2 +d2z(x,y)/dy2 =f(x,y)
其中xy为自变量，z为所求的表面高度，f(x,y)为z对于x，y的偏微分参数，即特征表面（保证G0连续）。
将等式写成有限差分形式：
d2z(i,j)/dx2 +d2z(i,j)/dy2 = 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迭代:
zn+1(i,j)=1/4(zn(i-1,j)+zn(i,j-1)+zn(i+1,j)+zn(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 }
