GLSL Matrix/Inverse multiplication precision

Question

I am trying to do some work with the GPU for a cloth simulation and I am having some problems with different hardware. I am using threejs as the framework, but I believe that isn't relevant for the problem I am having.

Basically what I do is I upload a matrix and an inverse of that matrix in order to transform points from local coords to world, do some math in world coords (like collision detection), and then transform them back to local. This works great on my laptop when I am using floating point textures, however I noticed on my phone there are some bizarre artifacts:

correct:

incorrect:

After doing some debugging I have narrowed it down to two problems. Both of them are related to decimal precision. The collapsing of the vertices due to constraints (and precision problems during the constraints) and losing precision when using the matrix multiplication and the inverse.

The reason I believe the problem is related to precision is because if I use a floating point texture it works on my computer, but if I use a half float I have the same problems. My phone supports floating point textures which is one reason I am confused about why this would be happening on my phone. I narrowed the problem down so all of the cloth simulation is disabled and if I run the application with half float textures on my computer, without any gravity but with the transformation and inverse the plane sort of flickers in weird ways

while if the transformation and inverse are disabled then it looks normal.

I am out of ideas for how to handle this problem though or if I am even going down the right path. I believe that half float textures have limited decimal precision but I don't understand why this would be causing my problems since it should only affect the output of the shader, not the math going on in the shader.

The code for the shader looks like the following:

    '   vec2 cellSize  = 1.0 / res;',
    '   vec4 pos = texture2D(vertexPositions, vuv.xy );',


    '   vec2 newUV;',
    '   if(type == 0.0){',
        '   float px = floor(vuv.x * res.x );',
        '   float spacingx = px- (2.0 * floor(px/2.0));',
        '   float py = floor(vuv.y * res.y );',
        '   float spacingy = py- (2.0 * floor(py/2.0));',
        '   float total = spacingx + spacingy;',
        '   total = total- (2.0 * floor(total/2.0));',

        '   if(total == 0.0){',
        '       newUV = vuv + (direction * cellSize);',
        '   }',
        '   else{',
        '       newUV = vuv - (direction * cellSize);',
        '   }',
    '   }',
    '   if(type == 1.0){',
        '   float px = floor(vuv.x * res.x );',
        '   float spacingx = px- (2.0 * floor(px/2.0));',

        '   float total = spacingx;',


        '   if(total == 0.0){',
        '       newUV = vuv + (direction * cellSize);',
        '   }',
        '   else{',
        '       newUV = vuv - (direction * cellSize);',
        '   }',
    '   }',






    '   vec4 totalDisplacement = vec4(0.0);',

    '           if(newUV.x > 0.0 && newUV.x < 1.0 && newUV.y > 0.0 && newUV.y < 1.0){ ',
    '               vec4 posOld = texture2D(vertexPositionsStart, vuv);' ,
    '               vec4 posOld2 = texture2D(vertexPositionsStart, newUV);' ,

    '               float targetDistance = length(posOld - posOld2);',
    '               vec4 newPos =  texture2D(vertexPositions, newUV);',
    '               float dx = pos.x - newPos.x;',
    '               float dy = pos.y - newPos.y;',
    '               float dz = pos.z - newPos.z;',
    '               float distance = sqrt(dx * dx + dy * dy + dz * dz);',
    '               float difference = targetDistance- distance;',
    '               float percent = difference / distance / 2.0;',
    '               float offsetX = dx * percent * rigid;',
    '               float offsetY = dy * percent * rigid;',
    '               float offsetZ = dz * percent * rigid;',
    '               totalDisplacement.x += offsetX;',
    '               totalDisplacement.y += offsetY;',
    '               totalDisplacement.z += offsetZ;',
    '           }',
    '       }',
    '   }',

    '   pos += totalDisplacement;',
    '   if(  vuv.x  > 1.0 - cellSize.x  && topConstrain == 1 ){',
    '       pos =transformation *  texture2D(vertexPositionsStart, vuv.xy );',
    '   }',

    '   if(  vuv.x  < cellSize.x  && bottomConstrain == 1 ){',
    '       pos =transformation *  texture2D(vertexPositionsStart, vuv.xy );',
    '   }',

    '   if(  vuv.y  < cellSize.y  && leftConstrain == 1 ){',
    '       pos =transformation *  texture2D(vertexPositionsStart, vuv.xy );',
    '   }',


    '   if(  vuv.y  > 1.0 - cellSize.y && rightConstrain == 1 ){',
    '       pos =transformation *  texture2D(vertexPositionsStart, vuv.xy );',
    '   }',




    '   gl_FragColor = vec4( pos.xyz , 1.0 );',

Answer 1

To ensure that your shader creates high precision variables for floating point calculations, the following should be added to the beginning of your vertex shader:

precision highp float;
precision highp int;

And in the fragment shader, floating point variable declarations should be declared as follows:

precision highp float;

Floating point errors are magnified if you utilise values in your calculations which are the result of previous calculations stored as floats. Otherwise known as intermediate values.

To minimise these errors, you should limit the number of intermediate calculations performed in the shader. For instance, you can fully expand the calculations for newUV:

newUV = vuv + ( direction * ( 1.0 / res ) );

You can also fully expand the calculation for totalDisplacement, step by step, first substitute offsetX like so:

totalDisplacement.x += ( dx * percent * rigid )

Now substitute each variable dx, and percent into the above:

totalDisplacement.x += ( ( pos.x - newPos.x ) * ( difference / distance / 2.0 ) * rigid )

You can now see that the equation can be expanded even further, substitute in difference like so:

totalDisplacement.x += ( ( pos.x - newPos.x ) * ( ( targetDistance- distance ) / ( distance * 2.0 ) ) * rigid );

At this point you can do some algebra to simplify and cancel out some variables (division by distance). By simplifying the above equation we now get the following:

totalDisplacement.x += ( ( pos.x - newPos.x ) * ( ( targetDistance / ( distance * 2.0 ) - 0.5 ) * rigid );

Finally we can substitute the formula for the targetDistance like so:

totalDisplacement.x += ( ( pos.x - newPos.x ) * ( ( length(posOld - posOld2) / ( distance * 2.0 ) - 0.5 ) * rigid );

And respectively for the other coordinates:

totalDisplacement.y += ( ( pos.y - newPos.y ) * ( ( length(posOld - posOld2) / ( distance * 2.0 ) - 0.5 ) * rigid );
totalDisplacement.z += ( ( pos.z - newPos.z ) * ( ( length(posOld - posOld2) / ( distance * 2.0 ) - 0.5 ) * rigid );

And obviously you can keep going, i.e. by substituting in the value for posOld, posOld2 and newPos.

Note now that up to this point, to get the equation above we have eliminated the need to store 5 of the intermediate values in float variables. Also through the simplification of the equation, (division by distance), the distance variable is only used in the calculation once. Compare this to the initial implementation, distance is used to calculate difference and percent. Combining it all in one equation allows you to simplify and reduce the amount of times the same floating point value is used. Hence also reducing the total floating point error. The trade off here is that the resulting equation is less human readable.

If you are curious, you can also inspect the level of precision for a given shader compiler by calling glGetShaderPrecisionFormat :

int range[2], precision;
glGetShaderPrecisionFormat(GL_FRAGMENT_SHADER, GL_HIGH_FLOAT, range, &precision);

If you inspect the result of this for the desktop and mobile version of your application, you will be able to compare the difference in precision for the two versions.