vertex multigrid

src/main.c

@@ -25,6 +25,7 @@
 #include"multigrid_correct_shader.h"
 #include"subtract_pressure_shader.h"
 #include"splat_shader.h"
+#include"blit_residue_fragment_shader.h"
 
 void* memcpy(void* dst, void* src, unsigned long n)
 {
@@ -134,6 +135,16 @@ typedef struct jacobi_program
 
 } jacobi_program;
 
+typedef struct blit_residue_program
+{
+	GLuint prog;
+
+	GLint pos;
+	GLint mvp;
+	GLint xTex;
+	GLint bTex;
+} blit_residue_program;
+
 typedef struct multigrid_restrict_residual_program
 {
 	GLuint prog;
@@ -204,6 +215,7 @@ subtract_program subtractProgram;
 splat_program splatProgram;
 blit_program blitProgram;
 blit_program blitDivProgram;
+blit_residue_program blitResidueProgram;
 
 frame_buffer colorBuffer;
 frame_buffer velocityBuffer;
@@ -325,6 +337,16 @@ void init_blit_div(blit_program* program)
 	program->tex = glGetUniformLocation(program->prog, "tex");
 }
 
+void init_blit_residue(blit_residue_program* program)
+{
+	console_log_fmt("compiling blit residue...");
+	program->prog = compile_shader(blit_vertex_shader_src, blit_residue_fragment_shader_src);
+	program->pos = glGetAttribLocation(program->prog, "pos");
+	program->mvp = glGetUniformLocation(program->prog, "mvp");
+	program->xTex = glGetUniformLocation(program->prog, "xTex");
+	program->bTex = glGetUniformLocation(program->prog, "bTex");
+}
+
 
 GLuint create_texture(int width, int height, GLenum internalFormat, GLenum format, GLenum type, char* initData)
 {
@@ -384,12 +406,9 @@ void frame_buffer_swap(frame_buffer* buffer)
 //----------------------------------------------------------------
 //NOTE(martin): entry point
 //----------------------------------------------------------------
-/*
-#define texWidth 512
-#define texHeight 512
-*/
-#define texWidth (128)
-#define texHeight (128)
+
+#define texWidth (256)
+#define texHeight (256)
 
 float colorInitData[texWidth][texHeight][4] = {0};
 float velocityInitData[texWidth][texHeight][4] = {0};
@@ -408,11 +427,11 @@ void reset_texture(GLuint texture, float width, float height, char* initData)
 	glTexImage2D(GL_TEXTURE_2D, 0, TEX_INTERNAL_FORMAT, width, height, 0, TEX_FORMAT, TEX_TYPE, initData);
 }
 
-static bool resetCmd = false;
+//static bool resetCmd = false;
 
 void reset()
 {
-	resetCmd = true;
+//	resetCmd = true;
 	console_log_fmt("reset");
 
 	reset_texture(colorBuffer.textures[0], texWidth, texHeight, (char*)colorInitData);
@@ -490,6 +509,7 @@ void init_velocity_vortex()
 			float y = 2*i/(float)texWidth - 1;
 			velocityInitData[i][j][0] = sinf(2*PI*y);
 			velocityInitData[i][j][1] = sinf(2*PI*x);
+
 		}
 	}
 }
@@ -504,8 +524,8 @@ int init(float canvasSize)
 
 	WAJS_SetupCanvas(canvasSize, canvasSize);
 
-	init_color_checker();
-	init_velocity_vortex();
+//	init_color_checker();
+//	init_velocity_vortex();
 
 	// init programs
 	init_advect(&advectProgram);
@@ -513,6 +533,7 @@ int init(float canvasSize)
 	init_jacobi(&jacobiProgram);
 	init_multigrid_restrict_residual(&multigridRestrictResidualProgram);
 	init_multigrid_correct(&multigridCorrectProgram);
+	init_blit_residue(&blitResidueProgram);
 
 	init_subtract(&subtractProgram);
 	init_splat(&splatProgram);
@@ -565,6 +586,56 @@ int init(float canvasSize)
 	return(0);
 }
 
+void apply_splat(float splatPosX, float splatPosY, float splatVelX, float splatVelY, float r, float g, float b, bool randomize)
+{
+	glUseProgram(splatProgram.prog);
+
+	if(randomize)
+	{
+		glUniform1f(splatProgram.randomize, 1.);
+	}
+	else
+	{
+		glUniform1f(splatProgram.randomize, 0.);
+	}
+
+	// force
+	glBindFramebuffer(GL_FRAMEBUFFER, velocityBuffer.fbos[1]);
+
+	glActiveTexture(GL_TEXTURE0);
+	glBindTexture(GL_TEXTURE_2D, velocityBuffer.textures[0]);
+	glUniform1i(splatProgram.src, 0);
+
+	glUniform2f(splatProgram.splatPos, splatPosX, splatPosY);
+	glUniform3f(splatProgram.splatColor, splatVelX, splatVelY, 0);
+	glUniform1f(splatProgram.additive, 1);
+	glUniform1f(splatProgram.blending, 0);
+
+	glUniform1f(splatProgram.radius, 0.01);
+
+	glDrawArrays(GL_TRIANGLES, 0, 6);
+
+	frame_buffer_swap(&velocityBuffer);
+
+	// dye
+	glBindFramebuffer(GL_FRAMEBUFFER, colorBuffer.fbos[1]);
+
+	glActiveTexture(GL_TEXTURE0);
+	glBindTexture(GL_TEXTURE_2D, colorBuffer.textures[0]);
+	glUniform1i(splatProgram.src, 0);
+
+	glUniform2f(splatProgram.splatPos, splatPosX, splatPosY);
+	glUniform3f(splatProgram.splatColor, r, g, b);
+	glUniform1f(splatProgram.additive, 0);
+	glUniform1f(splatProgram.blending, 1);
+	glUniform1f(splatProgram.radius, 0.01);
+
+	glDrawArrays(GL_TRIANGLES, 0, 6);
+
+	frame_buffer_swap(&colorBuffer);
+
+}
+
 void jacobi_solve(frame_buffer* x, frame_buffer* b, float invGridSize, int iterationCount)
 {
 	glUseProgram(jacobiProgram.prog);
@@ -587,6 +658,79 @@ void jacobi_solve(frame_buffer* x, frame_buffer* b, float invGridSize, int itera
 	}
 }
 
+void multigrid_coarsen_residual(frame_buffer* output, frame_buffer* x, frame_buffer* b, float invFineGridSize)
+{
+	//NOTE: compute residual and downsample to coarser grid, put result in coarser buffer
+	glUseProgram(multigridRestrictResidualProgram.prog);
+	glBindFramebuffer(GL_FRAMEBUFFER, output->fbos[1]);
+
+	glActiveTexture(GL_TEXTURE0);
+	glBindTexture(GL_TEXTURE_2D, x->textures[0]);
+	glUniform1i(multigridRestrictResidualProgram.xTex, 0);
+
+	glActiveTexture(GL_TEXTURE1);
+	glBindTexture(GL_TEXTURE_2D, b->textures[0]);
+	glUniform1i(multigridRestrictResidualProgram.bTex, 1);
+
+	glDrawArrays(GL_TRIANGLES, 0, 6);
+
+	frame_buffer_swap(output);
+}
+
+void multigrid_prolongate_and_correct(frame_buffer* x, frame_buffer* error, float invFineGridSize)
+{
+	//NOTE: correct finer pressure
+	glUseProgram(multigridCorrectProgram.prog);
+	glBindFramebuffer(GL_FRAMEBUFFER, x->fbos[1]);
+
+	glActiveTexture(GL_TEXTURE0);
+	glBindTexture(GL_TEXTURE_2D, x->textures[0]);
+	glUniform1i(multigridCorrectProgram.src, 0);
+
+	glActiveTexture(GL_TEXTURE1);
+	glBindTexture(GL_TEXTURE_2D, error->textures[0]);
+	glUniform1i(multigridCorrectProgram.error, 1);
+
+	glUniform1f(multigridCorrectProgram.invGridSize, invFineGridSize);
+
+	glDrawArrays(GL_TRIANGLES, 0, 6);
+
+	frame_buffer_swap(x);
+}
+
+void multigrid_clear(frame_buffer* error)
+{
+	glBindFramebuffer(GL_FRAMEBUFFER, error->fbos[0]);
+	glClear(GL_COLOR_BUFFER_BIT);
+}
+
+void input_splat(float t)
+{
+	//NOTE: apply force and dye
+	if(mouseInput.down && (mouseInput.deltaX || mouseInput.deltaY))
+	{
+		float canvasWidth = canvas_width();
+		float canvasHeight = canvas_height();
+
+		float splatPosX = mouseInput.x/canvasWidth;
+		float splatPosY = 1 - mouseInput.y/canvasHeight;
+
+		float splatVelX = 10000.*DELTA*mouseInput.deltaX/canvasWidth;
+		float splatVelY = -10000.*DELTA*mouseInput.deltaY/canvasWidth;
+
+		float intensity = 100*sqrtf(square(mouseInput.deltaX/canvasWidth) + square(mouseInput.deltaY/canvasHeight));
+
+		float r = intensity * (sinf(2*PI*0.1*t) + 1);
+		float g = 0.5*intensity * (cosf(2*PI*0.1/EULER*t + 654) + 1);
+		float b = intensity * (sinf(2*PI*0.1/SQRT2*t + 937) + 1);
+
+		apply_splat(splatPosX, splatPosY, splatVelX, splatVelY, r, g, b, false);
+
+		mouseInput.deltaX = 0;
+		mouseInput.deltaY = 0;
+	}
+}
+
 // This function is called by loader.js every frame
 void WAFNDraw()
 {
@@ -617,6 +761,61 @@ void WAFNDraw()
 
 	frame_buffer_swap(&velocityBuffer);
 
+	/*
+	//DEBUG
+	static bool splatTrig = false;
+	static bool splat = false;
+	static float splatStart = 0;
+	static int splatDir = 0;
+
+	static int frameCount = 0;
+
+	if(resetCmd)
+	{
+		frameCount = 0;
+		splat = true;
+		splatStart = frameT;
+	}
+
+	if(splat)
+	{
+		if(frameT - splatStart >= 0.5)
+		{
+			splat = false;
+			splatDir++;
+			splatDir = splatDir % 3;
+		}
+		float dirX = 0;
+		float dirY = 0;
+		if(splatDir == 0)
+		{
+			dirX = 0;
+			dirY = 0.3;
+		}
+		if(splatDir == 1)
+		{
+			dirX = 0.3;
+			dirY = 0;
+		}
+		if(splatDir == 2)
+		{
+			dirX = 0.2121;
+			dirY = 0.2121;
+		}
+		apply_splat(0.5, 0.5, dirX, dirY, 1.5, 1., 0.1, false);
+	}
+	resetCmd = false;
+
+	if(frameCount>20)
+	{
+		return;
+	}
+	frameCount++;
+*/
+
+	input_splat(t);
+
+
 	//NOTE: compute divergence of advected velocity
 	glUseProgram(divProgram.prog);
 	glBindFramebuffer(GL_FRAMEBUFFER, divBuffer[0].fbos[1]);
@@ -633,7 +832,34 @@ void WAFNDraw()
 	glBindFramebuffer(GL_FRAMEBUFFER, pressureBuffer[0].fbos[1]);
 	glClear(GL_COLOR_BUFFER_BIT);
 
-	jacobi_solve(&pressureBuffer[0], &divBuffer[0], INV_GRID_SIZE, 20000);
+	#if 0
+		multigrid_clear(&pressureBuffer[0]);
+		jacobi_solve(&pressureBuffer[0], &divBuffer[0], INV_GRID_SIZE, texWidth*texHeight);
+	#else
+		multigrid_clear(&pressureBuffer[0]);
+
+		for(int i=0; i<2; i++)
+		{
+			jacobi_solve(&pressureBuffer[0], &divBuffer[0], INV_GRID_SIZE, 2);
+			multigrid_coarsen_residual(&divBuffer[1], &pressureBuffer[0], &divBuffer[0], INV_GRID_SIZE);
+
+			multigrid_clear(&pressureBuffer[1]);
+			jacobi_solve(&pressureBuffer[1], &divBuffer[1], 2*INV_GRID_SIZE, 2);
+			multigrid_coarsen_residual(&divBuffer[2], &pressureBuffer[1], &divBuffer[1], 2*INV_GRID_SIZE);
+
+			multigrid_clear(&pressureBuffer[2]);
+			jacobi_solve(&pressureBuffer[2], &divBuffer[2], 4*INV_GRID_SIZE, 40);
+
+			multigrid_prolongate_and_correct(&pressureBuffer[1], &pressureBuffer[2], 2*INV_GRID_SIZE);
+			jacobi_solve(&pressureBuffer[1], &divBuffer[1], 2*INV_GRID_SIZE, 8);
+
+			multigrid_prolongate_and_correct(&pressureBuffer[0], &pressureBuffer[1], INV_GRID_SIZE);
+			jacobi_solve(&pressureBuffer[0], &divBuffer[0], INV_GRID_SIZE, 4);
+		}
+	#endif
+
+
+
 
 	//NOTE: subtract pressure gradient to advected velocity
 	glUseProgram(subtractProgram.prog);
@@ -667,7 +893,7 @@ void WAFNDraw()
 
 	glUniform1f(advectProgram.delta, DELTA);
 
-	glUniform1f(advectProgram.dissipation, 0.3);
+	glUniform1f(advectProgram.dissipation, 0.001);
 
 	glDrawArrays(GL_TRIANGLES, 0, 6);
 
@@ -675,6 +901,9 @@ void WAFNDraw()
 
 	//NOTE: Blit color texture to screen
 
+	//NOTE: blit residue to screen
+	glViewport(0, 0, canvas_width(), canvas_height());
+
 	float displayMatrix[16] = {
 		1/aspectRatio, 0, 0, 0,
 		0, 1, 0, 0,
@@ -682,8 +911,22 @@ void WAFNDraw()
 		0, 0, 0, 1 };
 
 /*
-	glViewport(0, 0, canvas_width(), canvas_height());
+	glUseProgram(blitResidueProgram.prog);
+	glBindFramebuffer(GL_FRAMEBUFFER, 0);
+
+	glActiveTexture(GL_TEXTURE0);
+	glBindTexture(GL_TEXTURE_2D, pressureBuffer[0].textures[0]);
+	glUniform1i(blitResidueProgram.xTex, 0);
 
+	glActiveTexture(GL_TEXTURE1);
+	glBindTexture(GL_TEXTURE_2D, divBuffer[0].textures[0]);
+	glUniform1i(blitResidueProgram.bTex, 1);
+
+	glUniformMatrix4fv(blitResidueProgram.mvp, 1, GL_FALSE, displayMatrix);
+
+	glDrawArrays(GL_TRIANGLES, 0, 6);
+//*/
+//*
 	glUseProgram(blitProgram.prog);
 	glBindFramebuffer(GL_FRAMEBUFFER, 0);
 
@@ -695,6 +938,7 @@ void WAFNDraw()
 
 	glDrawArrays(GL_TRIANGLES, 0, 6);
 /*/
+
 	//NOTE: recompute divergence of (corrected) velocity
 	glUseProgram(divProgram.prog);
 	glBindFramebuffer(GL_FRAMEBUFFER, divBuffer[0].fbos[1]);
@@ -719,5 +963,6 @@ void WAFNDraw()
 	glUniformMatrix4fv(blitDivProgram.mvp, 1, GL_FALSE, displayMatrix);
 
 	glDrawArrays(GL_TRIANGLES, 0, 6);
+
 //*/
 }

src/shaders/advect.glsl

@@ -13,22 +13,14 @@ uniform float dissipation;
 
 vec2 u(ivec2 coord)
 {
-	if(  coord.x <= 0
-	  || coord.x >= textureSize(velocity, 0).x
-	  || coord.y <= 0
-	  || coord.y >= textureSize(velocity, 0).y)
-	{
-		return(vec2(0.));
-	}
-
 	return(texelFetch(velocity, coord, 0).xy);
 }
 
 vec4 q(ivec2 coord)
 {
-	if(  coord.x <= 0
+	if(  coord.x < 0
 	  || coord.x >= textureSize(src, 0).x
-	  || coord.y <= 0
+	  || coord.y < 0
 	  || coord.y >= textureSize(src, 0).y)
 	{
 		return(vec4(0.));

src/shaders/blit_residue_fragment.glsl

@@ -0,0 +1,65 @@
+#version 300 es
+
+precision highp float;
+precision highp sampler2D;
+
+in vec2 texCoord;
+out vec4 fragColor;
+
+uniform sampler2D xTex;
+uniform sampler2D bTex;
+
+float x(ivec2 coord)
+{
+	if(  coord.x <= 0
+	  || coord.x >= textureSize(xTex, 0).x
+	  || coord.y <= 0
+	  || coord.y >= textureSize(xTex, 0).y)
+	{
+		return(0.);
+	}
+	return(texelFetch(xTex, coord, 0).x);
+}
+
+float b(ivec2 coord)
+{
+	if(  coord.x <= 0
+	  || coord.x >= textureSize(bTex, 0).x
+	  || coord.y <= 0
+	  || coord.y >= textureSize(bTex, 0).y)
+	{
+		return(0.);
+	}
+	return(texelFetch(bTex, coord, 0).x);
+}
+
+vec3 color_map(float v)
+{
+	float logv = log(abs(v))/log(10.0);
+	float f = floor(logv + 7.0);
+	float i = floor(4.0*(logv + 7.0 - f));
+
+	if(f < 0.0) return vec3(0.0);
+	if(f < 1.0) return mix(vec3(1.0, 0.0, 0.0), vec3(1.0), i/4.0);
+	if(f < 2.0) return mix(vec3(0.0, 1.0, 0.0), vec3(1.0), i/4.0);
+	if(f < 3.0) return mix(vec3(0.0, 0.0, 1.0), vec3(1.0), i/4.0);
+	if(f < 4.0) return mix(vec3(1.0, 1.0, 0.0), vec3(1.0), i/4.0);
+	if(f < 5.0) return mix(vec3(1.0, 0.0, 1.0), vec3(1.0), i/4.0);
+	if(f < 6.0) return mix(vec3(0.0, 1.0, 1.0), vec3(1.0), i/4.0);
+	if(f < 7.0) return mix(vec3(1.0, 0.5, 0.0), vec3(1.0), i/4.0);
+	if(f < 8.0) return mix(vec3(1.0, 1.0, 1.0), vec3(1.0), i/4.0);
+	return vec3(1.0);
+}
+
+void main()
+{
+	ivec2 pixelCoord = ivec2(floor(texCoord.xy * vec2(textureSize(xTex, 0).xy)));
+
+	float tl = x(pixelCoord + ivec2(-1, 1));
+	float tr = x(pixelCoord + ivec2(1, 1));
+	float bl = x(pixelCoord + ivec2(-1, -1));
+	float br = x(pixelCoord + ivec2(1, -1));
+
+	float residue = b(pixelCoord) - (-tl - tr - bl - br + 4.*x(pixelCoord));
+	fragColor = vec4(color_map(residue), 1);
+}

src/shaders/divergence.glsl

@@ -8,28 +8,9 @@ out vec4 fragColor;
 
 uniform sampler2D src;
 
-float ux(ivec2 coord)
+vec2 u(ivec2 coord)
 {
-	if(  coord.x <= 0
-	  || coord.x >= textureSize(src, 0).x
-	  || coord.y <= 0
-	  || coord.y >= textureSize(src, 0).y)
-	{
-		return(0.);
-	}
-	return(texelFetch(src, coord, 0).x);
-}
-
-float uy(ivec2 coord)
-{
-	if(  coord.x <= 0
-	  || coord.x >= textureSize(src, 0).x
-	  || coord.y <= 0
-	  || coord.y >= textureSize(src, 0).y)
-	{
-		return(0.);
-	}
-	return(texelFetch(src, coord, 0).y);
+	return(texelFetch(src, coord, 0).xy);
 }
 
 void main()
@@ -45,16 +26,16 @@ void main()
 	}
 	else
 	{
-		ivec2 vr = pixelCoord + ivec2(1, 0);
-		ivec2 vl = pixelCoord - ivec2(1, 0);
-		ivec2 vt = pixelCoord + ivec2(0, 1);
-		ivec2 vb = pixelCoord - ivec2(0, 1);
+		vec2 tl = u(pixelCoord + ivec2(-1, 0));
+		vec2 tr = u(pixelCoord);
+		vec2 bl = u(pixelCoord + ivec2(-1, -1));
+		vec2 br = u(pixelCoord + ivec2(0, -1));
 
-		float r = ux(vr);
-		float l = ux(vl);
-		float t = uy(vt);
-		float b = uy(vb);
+		float r = (tr.x + br.x)/2.;
+		float l = (tl.x + bl.x)/2.;
+		float t = (tl.y + tr.y)/2.;
+		float b = (bl.y + br.y)/2.;
 
-		fragColor = vec4(-(r - l + t - b)/2., 0, 0, 1);
+		fragColor = vec4(-2.*(r - l + t - b), 0, 0, 1);
 	}
 }

src/shaders/jacobi_step.glsl

@@ -38,23 +38,18 @@ void main()
 	ivec2 pixelCoord = ivec2(floor(gl_FragCoord.xy));
 
 	if(  pixelCoord.x <= 0
-	  || pixelCoord.x >= textureSize(xTex, 0).x
-	  || pixelCoord.y <= 0
-	  || pixelCoord.y >= textureSize(xTex, 0).y)
+	  || pixelCoord.y <= 0)
 	{
 		fragColor = vec4(0, 0, 0, 1);
 	}
 	else
 	{
-		ivec2 vr = pixelCoord + ivec2(1, 0);
-		ivec2 vl = pixelCoord - ivec2(1, 0);
-		ivec2 vt = pixelCoord + ivec2(0, 1);
-		ivec2 vb = pixelCoord - ivec2(0, 1);
+		float tl = x(pixelCoord + ivec2(-1, 1));
+		float tr = x(pixelCoord + ivec2(1, 1));
+		float bl = x(pixelCoord + ivec2(-1, -1));
+		float br = x(pixelCoord + ivec2(1, -1));
 
-		float w = 0.8;
-		float standardJacobi = (x(vl) + x(vr) + x(vt) + x(vb) + b(pixelCoord))/4.;
-		float weightedJacobi = w*standardJacobi + (1.-w)*x(pixelCoord);
-
-		fragColor = vec4(weightedJacobi, 0, 0, 1);
+		float jacobi = (tl + tr + bl + br + b(pixelCoord))/4.;
+		fragColor = vec4(jacobi, 0, 0, 1);
 	}
 }

src/shaders/splat.glsl

@@ -21,29 +21,8 @@ void main()
 	float intensity = exp(-10.*d2/radius);
 	vec2 force = splatColor.xy;
 
-	if(randomize != 0.0)
-	{
-		float n = 119.*gl_FragCoord.x + 107.*gl_FragCoord.y + 113.;
-		n = fract(n*fract(n/2.7182818));
-		intensity += intensity*abs(cos(2.*3.14159*n));
-
-		//NOTE: the force is added in the additive part.
-
-		force.x += 0.2*sin(2.*3.14159*fract(n/1.144213));
-		force.y += 0.2*cos(2.*3.14159*fract(n/1.144213));
-
-		{
-			vec2 u = force/length(force);
-
-			vec2 v = texCoord - splatPos;
-			vec2 norm = v - dot(v, u)*u;
-			float dFromAxis = length(norm);
-
-			force += 30.*norm;
- 		}
-	}
 	vec3 u = texture(src, texCoord).xyz;
-	vec3 uAdd = u + intensity*vec3(force, 0);
+	vec3 uAdd = u + intensity*splatColor.xyz;
 	vec3 uBlend = u*(1.-intensity) + intensity * splatColor;
 
 	fragColor = vec4(uAdd*additive + uBlend*blending, 1);

src/shaders/subtract_pressure.glsl

@@ -12,13 +12,6 @@ uniform float invGridSize;
 
 vec2 u(ivec2 coord)
 {
-	if(  coord.x <= 0
-	  || coord.x >= textureSize(src, 0).x
-	  || coord.y <= 0
-	  || coord.y >= textureSize(src, 0).y)
-	{
-		return(vec2(0.));
-	}
 	return(texelFetch(src, coord, 0).xy);
 }
 
@@ -38,22 +31,17 @@ void main()
 {
 	ivec2 pixelCoord = ivec2(floor(gl_FragCoord.xy));
 
-	if(  pixelCoord.x <= 0
-	  || pixelCoord.x >= textureSize(src, 0).x
-	  || pixelCoord.y <= 0
-	  || pixelCoord.y >= textureSize(src, 0).y)
-	{
-		fragColor = vec4(0, 0, 0, 1);
-	}
-	else
-	{
-		ivec2 vr = pixelCoord + ivec2(1, 0);
-		ivec2 vl = pixelCoord - ivec2(1, 0);
-		ivec2 vt = pixelCoord + ivec2(0, 1);
-		ivec2 vb = pixelCoord - ivec2(0, 1);
+	float tl = p(pixelCoord + ivec2(0, 1));
+	float tr = p(pixelCoord + ivec2(1, 1));
+	float bl = p(pixelCoord);
+	float br = p(pixelCoord + ivec2(1, 0));
 
-		vec2 gradP = vec2(p(vr) - p(vl), p(vt) - p(vb))/2.;
+	float r = (tr + br)/2.;
+	float l = (tl + bl)/2.;
+	float t = (tl + tr)/2.;
+	float b = (bl + br)/2.;
 
-		fragColor = vec4(u(pixelCoord) - gradP, 0, 1);
-	}
+	vec2 gradP = vec2(r - l, t - b);
+
+	fragColor = vec4(u(pixelCoord) - gradP, 0, 1);
 }