Bug when using single precsion

[bkaradzic]

I ran into issue where calling function convhull_3d_build multiple times with the same input causes different results. The bug exist only when used with CONVHULL_3D_USE_SINGLE_PRECISION.

Repro:

#include <stdlib.h>

#define CONVHULL_3D_ENABLE
#define CONVHULL_3D_USE_SINGLE_PRECISION /* optional */
#include "../convhull_3d.h"

int main(int argc, const char * argv[])
{
	ch_vec3 vert[] =
	{
		{  4.243597,  2.029179, -0.815667 },
		{  4.417669,  0.039982, -0.928770 },
		{ -2.428492,  1.198822,  3.519505 },
		{ -2.254420, -0.790374,  3.406402 },
		{  3.153918,  2.029179, -2.492747 },
		{  3.327990,  0.039982, -2.605850 },
		{ -3.518171,  1.198822,  1.842424 },
		{ -3.344099, -0.790374,  1.729321 },
	};

	for (int ii = 0; ii < 10; ++ii)
	{
		printf("%d\n", ii);
		int* indices = NULL;
		int numFaces;
		convhull_3d_build(vert, 8, &indices, &numFaces);
		free(indices);
	}

	return 0;
}

Program hits this assert:

convhull_3d/convhull_3d.h

Line 923 in 88eca3c

assert(detA>0.0);

> gcc bug.c -lm -o bug; and ./bug
0
1
2
bug: ../convhull_3d.h:923: convhull_3d_build: Assertion `detA>0.0' failed.

[bkaradzic]

I replaced rand with new random function and the issue is gone for this particular input.

CH_FLOAT rnd(CH_FLOAT x, CH_FLOAT y)
{
    // Reference(s):
    //
    // - Improvements to the canonical one-liner GLSL rand() for OpenGL ES 2.0
    //   http://byteblacksmith.com/improvements-to-the-canonical-one-liner-glsl-rand-for-opengl-es-2-0/
    //
    CH_FLOAT a  = 12.9898;
    CH_FLOAT b  = 78.233;
    CH_FLOAT c  = 43758.5453;
    CH_FLOAT dt = x*a + y*b;
#ifdef CONVHULL_3D_USE_SINGLE_PRECISION
    float sn = fmodf(dt, 3.14);
    float intpart;
    return modff(sin(sn) * c, &intpart);
#else
    double sn = fmod(dt, 3.14);
    double intpart;
    return modf(sin(sn) * c, &intpart);
#endif // CONVHULL_3D_USE_SINGLE_PRECISION
}

[bkaradzic]

Input that fails with Assertion 'detA>0.0' failed.:

	ch_vec3 vert[] =
	{
		{  2.531722,  6.652832, -0.009874 },
		{  3.408637,  5.542159, -1.423177 },
		{  0.189403, -0.000172,  3.765184 },
		{  1.066318, -1.110846,  2.351881 },
		{  0.832276,  6.652832, -1.064333 },
		{  1.709191,  5.542159, -2.477636 },
		{ -1.510042, -0.000172,  2.710725 },
		{ -0.633127, -1.110846,  1.297422 },
	};

Input in these tests is just box that rotates.

[leomccormack]

Hi! Thank you for letting me know, and for investigating it further. This looks like a problem of the numerical precision (or rather: in-precision) when trying to figure out the orientation of a face. This assert is basically to stop the loop from going on for infinity. Indeed, the introduction of a bit noise is inelegant, but I've not found a better way yet. So I am very open to suggestions. Note that the issue is significantly more rare when using double precision, which is why this is the default configuration. Double precision is also the default for MATLAB, so the algorithm in the original implementation (which convhull_3d is derived from) also used double precision:

https://se.mathworks.com/matlabcentral/fileexchange/48509-computational-geometry-toolbox?focused=3851286&tab=function

[bkaradzic]

the introduction of a bit noise is inelegant, but I've not found a better way yet

rnd function I proposed above does solve this issue with the same input not producing the same output. I can send PR with it as a fix.

Note that the issue is significantly more rare when using double precision

One thing I was thinking it to introduce another define. Right now define single or double precision affects input/output, and processing data. It would be ideal if input/output data could be single precission while processing happens as double precission.

Another idea (I didn't actually check much code, so it might already happen) option is to scale inputs into homogenous space do all calculation, and then restore it back into model scale.

[leomccormack]

@bkaradzic good suggestions! If rnd minimises the problem, then please feel free to send a PR.

And true, the interfacing with the function could be single precision while having the internals using double precision. I actually wrote this wrapper for another project, but it could be integrated into convhull_3d_build here.
Actually, there could also be e.g.: convhull_3d_build_float and convhull_3d_build_double functions instead of this define switch case, if that would be better?

And do you mean to scale-up the vertices values prior to building the convexhull? If it helps, then I say go for it : )

[bkaradzic]

Ok, I fixed random issue: #8

I tested moving vertices into origin into homogeneous -1, 1 space, but while it improved situation, it didn't fully solve the issue. I actually need to step thru to completely understand why this is happening.

Wrapper function does solve the issue on simple models I've been using. Probably good enough, but I'm really curious to dig into it, to figure out what actually is causing it. I feel my simple box, with 8 points should be error free even with single-precision, floating point range is pretty low, and it has no overlapping points.