From 890e2bad880a1ba9396fbe4ed0ebf9de87f576ec Mon Sep 17 00:00:00 2001 From: Andreas Grois Date: Fri, 9 Mar 2018 23:00:49 +0100 Subject: Make this thing render. Also: Steal integer hash from stackexchange... --- BuddhaTest/Shaders/BuddhaCompute.glsl | 185 ++++++++++++++++++++++++++++++++-- BuddhaTest/src/BuddhaTest.cpp | 11 +- 2 files changed, 186 insertions(+), 10 deletions(-) diff --git a/BuddhaTest/Shaders/BuddhaCompute.glsl b/BuddhaTest/Shaders/BuddhaCompute.glsl index 49209ae..4ca5141 100644 --- a/BuddhaTest/Shaders/BuddhaCompute.glsl +++ b/BuddhaTest/Shaders/BuddhaCompute.glsl @@ -9,7 +9,9 @@ layout(std430, binding=2) buffer renderedData uniform uint iterationCount; -void addToColorAt(uvec2 cell, uvec3 toAdd) +uniform uint orbitLength; + +void addToColorOfCell(uvec2 cell, uvec3 toAdd) { uint firstIndex = 3*(cell.x + cell.y * width); atomicAdd(counts_SSBO[firstIndex],toAdd.x); @@ -17,11 +19,182 @@ void addToColorAt(uvec2 cell, uvec3 toAdd) atomicAdd(counts_SSBO[firstIndex+2],toAdd.z); } +uvec2 getCell(vec2 complex) +{ + vec2 uv = clamp(vec2((complex.x+2.5)/3.5, 0.5*(complex.y + 1.0)),vec2(0.0),vec2(1.0)); + return uvec2(width * uv.x, height * uv.y); +} + +void addToColorAt(vec2 complex, uvec3 toAdd) +{ + uvec2 cell = getCell(complex); + addToColorOfCell(cell,toAdd); +} + +/* +vec2 hash2( uint n, out uint hash ) +{ + // integer hash copied from Hugo Elias + n = (n << 13U) ^ n; + n = n * (n * n * 15731U + 789221U) + 1376312589U; + hash = n*n; + uvec2 k = n * uvec2(n,n*16807U); + return vec2( k & uvec2(0x7fffffffU))/float(0x7fffffff); +} +*/ + +uint intHash(uint x) { + x = ((x >> 16) ^ x) * 0x45d9f3bU; + x = ((x >> 16) ^ x) * 0x45d9f3bU; + x = (x >> 16) ^ x; + return x; +} +vec2 hash2(uint n, out uint hash) +{ + uint ih =intHash(n); + hash = intHash(ih); + uvec2 k = uvec2(ih,hash); + return vec2(k & uvec2(0x7fffffffU))/float(0x7fffffff); +} + +vec2 compSqr(in vec2 v) +{ + return vec2(v.x*v.x-v.y*v.y, 2.0*v.x*v.y); +} + +float isInMainCardioid(vec2 v) +{ + /* + The condition that a point c is in the main cardioid is that its orbit has + an attracting fixed point. In other words, it must fulfill + z**2 -z + v = 0 (z**2+v has a fixed point at z) + and + d/dz(z**2+v) < 1 (fixed point at z is attractive) + Solving these equations yields + v = u/2*(1-u/2), where u is a complex number inside the unit circle + + Sadly we only know v, not u, and inverting this formula leads to a complex square root. + This is the old code that uses the compSqrt method, which is rather slow... + + vec2 toRoot = (vec2(1.0,0.0)-vec2(4.0)*v); + vec2 root = compSqrt(toRoot); + vec2 t1,t2; + t1=vec2(1,0)-root; + t2=vec2(1,0)+root; + float retval = step(dot(t1,t1),0.99999); + retval += step(dot(t2,t2),0.99999); + retval = min(retval,1.0); + return retval; + + On several websites (and several mandelbrot-related shaders on ShaderToy) one can + find various faster formulas that check the same inequality. + What however is hard to find is the actual derivation of those formulas, + and they are not as trivial as one might think at first. + That's why I'm writing this lengthy comment, to preserve the scrap notes I made + for future reuse by myself and others... + + Now the actual derivation looks like this: + We start with the following line + u = 1 +- sqrt(1-4*v) + don't mind the +-, that's just me being too lazy to check which solution is the right one + We know that |u| < 1, and we immediately square the whole beast to get rid of the root + Some definitions: r := sqrt(1-4*v), z = 1-4*v -layout (local_size_x = 1, local_size_y = 1) in; + 1 > Re(u)**2+Im(u)**2 = (1 +- Re(r))**2+Im(r)**2 + 1 > 1 +- 2*Re(r) + Re(r)**2+Im(r)**2 + 1 > 1 +- 2*Re(r) + |r|**2 + For complex values the square root of the norm is the same as the norm of the square root + 1 > 1 +- 2*Re(r) + |z| + +-2*Re(r) > |z| + 4*Re(r)**2 > |z|**2 + This step is now a bit arcane. If one solves the two coupled equations + (a+i*b) = (x+i*y)*(x+i*y) component-wise for x and y, + one can see that x**2 = (|a+i*b|+a)/2 + With this follows + 2*(|z|+Re(z)) > |z|**2 + |z| > 1/2*|z|**2-Re(z) + |z|**2 > (1/2*|z|**2-Re(z))**2 + + And long story short, the result is the following few operations. + */ + vec2 z = vec2(1.0,0.0)-4.0*v; + float zNormSqr = dot(z,z); + float rhsSqrt = 0.5*zNormSqr - z.x; + return step(rhsSqrt*rhsSqrt,zNormSqr); +} + +float isInMainBulb(vec2 v) +{ + //The condition for this is that f(f(z)) = z + //where f(z) = z*z+v + //This is an equation of fourth order, but two solutions are the same as + //for f(z) = z, which has been solved for the cardioid above. + //Factoring those out, you end up with a second order equation. + //Again, the solutions need to be attractive (|d/dz f(f(z))| < 1) + //Well, after a bit of magic one finds out it's a circle around -1, radius 1/4. + vec2 shifted = v + vec2(1,0); + float sqrRadius = dot(shifted,shifted); + return step(sqrRadius,0.062499999); +} + +vec2 getStartValue(uint seed) +{ + uint hash = seed; + + vec2 retval = vec2(0); + for(int i = 0; i < 5; ++i) + { + vec2 random = hash2(hash,hash); + vec2 point = vec2(random.x * 3.5-2.5,random.y*3.1-1.55); + float useThisPoint =1.0-(isInMainBulb(point) + isInMainCardioid(point)); + retval = mix(retval,point,useThisPoint); + } + return retval; +} + +bool isGoingToBeDrawn(vec2 offset) +{ + vec2 val = vec2(0); + for(int i = 0; i < orbitLength;++i) + { + val = compSqr(val) + offset; + if(dot(val,val) > 4.0) + { + return true; + } + } + return false; +} + +void drawOrbit(vec2 offset) +{ + vec2 val = vec2(0); + for(int i = 0; i < orbitLength;++i) + { + val = compSqr(val) + offset; + if(dot(val,val) > 20.0) + { + return; + } + if(val.x > -2.5 && val.x < 1.0 && val.y > -1.0 && val.y < 1.0) + addToColorAt(val,uvec3(1,1,1)); + } +} + +layout (local_size_x = 1024) in; //to be safe, we limit our local work group size to 1024. That's the minimum a GL 4.3 capable driver must support. void main() { - if(iterationCount < 1024) - addToColorAt(uvec2(iterationCount%width,iterationCount/width),uvec3(1,0,0)); - else - addToColorAt(uvec2(iterationCount%width,iterationCount/width),uvec3(0,1,0)); + //we need to know how many total work groups are running this iteration + + uvec3 totalWorkersPerDimension = gl_WorkGroupSize * gl_NumWorkGroups; + uint totalWorkers = totalWorkersPerDimension.x*totalWorkersPerDimension.y*totalWorkersPerDimension.z; + + //TODO: Check this once I've had some sleep. Anyhow, I'm using 1D, so y and z components globalInfocationID should be zero anyhow. + uint seed = iterationCount * totalWorkers + gl_GlobalInvocationID.x + gl_GlobalInvocationID.y*totalWorkersPerDimension.x + gl_GlobalInvocationID.z*(totalWorkersPerDimension.x + totalWorkersPerDimension.y); + + vec2 offset = getStartValue(seed); + + if(!isGoingToBeDrawn(offset)) + return; + + drawOrbit(offset); } diff --git a/BuddhaTest/src/BuddhaTest.cpp b/BuddhaTest/src/BuddhaTest.cpp index fcbf3bf..790b0b2 100644 --- a/BuddhaTest/src/BuddhaTest.cpp +++ b/BuddhaTest/src/BuddhaTest.cpp @@ -11,9 +11,10 @@ void error_callback(int error, const char* description) int main() { - unsigned int bufferWidth = 640; - unsigned int bufferHeight = 480; + unsigned int bufferWidth = 1600; + unsigned int bufferHeight = 900; + unsigned int orbitLength = 1000; GLFWwindow* window; @@ -28,7 +29,7 @@ int main() glfwWindowHint(GLFW_OPENGL_PROFILE, GLFW_OPENGL_CORE_PROFILE); /* Create a windowed mode window and its OpenGL context */ - window = glfwCreateWindow(640, 480, "Buddhabrot", NULL, NULL); + window = glfwCreateWindow(1600, 900, "Buddhabrot", NULL, NULL); if (!window) { std::cerr << "Failed to create OpenGL 4.3 core context. We do not support compatibility contexts." << std::endl; @@ -79,6 +80,8 @@ int main() uint32_t iterationCount{0}; glUseProgram(ComputeShader); GLint iterationCountUniformHandle = glGetUniformLocation(ComputeShader, "iterationCount"); + GLint orbitLengthUniformHandle = glGetUniformLocation(ComputeShader, "orbitLength"); + glUniform1ui(orbitLengthUniformHandle,orbitLength); glClearColor(0.0f, 0.0f, 0.4f, 0.0f); @@ -89,7 +92,7 @@ int main() glUseProgram(ComputeShader); //increase iterationCount, which is used for pseudo random generation glUniform1ui(iterationCountUniformHandle,++iterationCount); - glDispatchCompute(1920, 1080, 1); + glDispatchCompute(128, 1, 1); //before reading the values in the ssbo, we need a memory barrier: glMemoryBarrier(GL_SHADER_STORAGE_BARRIER_BIT); //I hope this is the correct (and only required) bit -- cgit v1.2.3