truncate float32 to lesser precision (22 bit mantissa)

[xiio...@gmail.com]

ok I have a problem with rounding errors on floats which I think is unavoidable. Just as specific background this is happening for me on tests of vectors (lines) intersecting with bounding boxes (ie oblongs, axis aligned). Two code snippets shows the core of the maths..

front_n = (x_front - s.X) / v.X

ripY = s.Y + front_n*v.Y

ripZ = s.Z + front_n*v.Z

if (front_n > 0) && (ymin <= ripY && ripY <= ymax) && (zmin <= ripZ && ripZ <= zmax) {

...

           etc

and so on eg :

back_n = (x_back - s.X) / v.X

ripY = s.Y + back_n*v.Y

ripZ = s.Z + back_n*v.Z

if (back_n > 0) && (ymin <= ripY && ripY <= ymax) && (zmin <= ripZ && ripZ <= zmax) {

...

          etc 

 

where xmin, ymin etc represent the boundaries of the box, anything .X , .Y etc is one component of a 3d vector, and anything suffix _n is a scalar. All are floats

One of the two boundaries (front or back face) is occasionally skipped .. The issue occurs obviously (?) when the vector intersects and crosses the box close to or on and edge/corner are reached. It's rare, and clearly (?) switching to float64 makes it less rare, but doesn't solve the problem. I considered multiplying out the divides but that greatly complicates the boundary conditions, and I don't think it actually solves the problem.. I think using a 'delta' in the boundaries just kicks the same problem down the road.

(I've got a simple workaround that doesn't solve the core problem..)

I think I need to truncate the accuracy of the float to 22 bits  - so the question is really about guard and round bits etc - we got anything supporting that in go ? It looked like I would have to use pkg unsafe

[Michael Jones]

Any kind of computational geometry code will have subtle issues in just this area. Broadly, the issue is classification of a point or line as IN, OUT, or ON. Of these, ‘ON’ is the difficult one because floating point prevents perfect calculation where real numbers would allow it. This is because floating point is actually integer math with a companion power function, and the integer math part has the problem in cases of division where the numerators and denominators are relatively prime. This leads to differences in edge and point positions and frustrates the ON calculation, and this the in and out calculation near the edge.

 

I suggest that you think of this carefully because there is generally no good solution to be found by considering only floating point precision—the same problem will exist at any non-infinite precision.

 

However, to answer your specific question, yes you can easily truncate float32 to a lower precision. Just force the bits you want to truncate to zero.

 

Hp is high precision

Lp is low precision

 

Lp := (unsafe interpret as float)(Mask & (unsafe interpret as integer(Hp))

 

Mask would be used to discard the lower bits of the word and thus precision.

 Mask := ^uint32(0) << k //discard lowest k bits

 

Nice examples of what this accomplishes are here:

http://www.h-schmidt.net/FloatConverter/IEEE754.html

 

 

--
You received this message because you are subscribed to the Google Groups "golang-nuts" group.
To unsubscribe from this group and stop receiving emails from it, send an email to golang-nuts...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[xiio...@gmail.com]

Yes thanks. 

I decided on catching the missed ON (edge/corner) case at the end of the function with an extra if statement when the numerical work returns EXITS_BOX=true, ENTERS_BOX=false - then I do a relatively simple check for the 'ray' outside the box - if Ray.Origin was outside and EXITS_BOX=true then ENTERS_BOX should have been true too - so a contradiction here lets be catch the miss. (There's another reverse case I don't need to catch).

if s.X < xmin || s.Y < ymin || s.Z < zmin || s.X > xmax || s.Y > ymax || s.Z > zmax {
		//origin is outside box *but* ray intersected with back plane
		//	so ray should have intersected with front plane
		//	probable rounding error at edge/corner intersection
		return true, true, back_n, back_n, back_p, back_p
	}

The float conversion is a little fiddly because of that implicit bit in the mantissa .. if the compiler is working well the catch is only 6 compares, and the rounding conversion adds a bit-ops AND, OR and SHIFT.. so I doubt the unsafe will make a huge difference in time.. though I think it could be preferable.

.. Mainly relieved to here someone say this problem is unavoidable with floats .. I can stop looking for maths work arounds.. Thanks again.

[jnc...@gmail.com]

Michael is of course correct about floating point but it is possible to robustly handle cases like this. I'm not sure on your exact goals but you've not already seen there are several papers about robust ray-box intersection e.g. An Efficient and Robust Ray-Box Intersection Algorithm by Amy Williams et al., and https://tavianator.com/fast-branchless-raybounding-box-intersections/ has some good information about handling precision, infs and NaNs (Vermeer LT implements some of this in Go but mainly in asm unfortunately so it's not much use for borrowing :). Also the chapter in PBRT about floating point error handling is very handy for understanding the issues (there is a preprint PDF floating).

Cheers,

Jamie

[xiio...@gmail.com]

Thanks for those links - I started using a something similar to what is described as Smit's algorithm (Slab method) in the Amy Williams paper. My code us subtly different to that usually given because I check the "back face" intersection first -for an early terminate ; this basically involves swapping x_plane_min and x_plane_max early if v.X<0 .. then looking explicitly for the back plane intersection points first..  as well as being slightly more verbose it seems to have the side effect of avoiding the divide by -0 problem described in the paper. (misses ray parallel to plane, catches adjoining edges) I think my version is sub-optimal but not a priority right now.

I've noticed that the divides in all the "slab" methods could be multiplied out - this which inverts the boundary condition if the divisor was negative, but can be caught using 

if (multiplied_x < undivided_boundary_value) ==(divisor>0)

instead of the vanilla.

if (x < boundary_value) 

not only does this avoid any division (hurray!), but also reduces the cumulative float errors - the number of operations on both sides of the condition of is roughly equal, rather than being unbalanced with most on one side in the un-premultuplied version.. It's a little more complicated, but divides are, afaik, still slow even on modern hardware, which is something I suspect may be a bottleneck within the function.. I need to bench this.. It may mitigate but doesn't solve the float precision problem.

The second link is a nice optimisation for simd too. Optimisation for this function isn't a big priority for me because I'm primarily using it it get the intersection point for a bounding volume for a "3d texture" (or 'non-sparse voxel array' in other words) - my main routine steps through the voxel grid front to back using a different algorithm.

PBRT looks good - wasn't aware of it - found the pre-publication chapter section 3.9.2 "CONSERVATIVE RAY–BOUNDS INTERSECTIONS" - though I didn't proof read it it looks like it covers using calculated expected errors (deltas) well.

[Nigel Tao]

On Sun, Sep 18, 2016 at 1:57 PM, <xiio...@gmail.com> wrote:
> It looked like I would have to use pkg unsafe

BTW, you don't have to use package unsafe. Use package math's
Float32bits and Float32frombits functions.

[Michael Jones]

Good point. I missed that. Also, I did not mean “unsolvabe” just not solvable using slightly a different precision.

Common solutions:

Remove the divides by multiplying them out. Then it just works (at least in FP where integer precision is not such an issue as in graphics hardware)

Avoid the clipping/division tests by doing comparisons alone.

Deal with relative primality by direct means. I think I have a patent (SGI patent / OpenGL) about dealing with this as part of the reference polygon extension. There was an elaboration that resolved the “cracks at T-vertices” problem too, but that died with SGI.

[jnc...@gmail.com]

Sure, as you say for your use case it might be best for you, and you might need different information from the test anyway.   Usually in a fast ray tracer you'd pre-calculate the inverse of the direction components and store them (since you do LOTS of ray-bbox intersections for each ray while traversing the bounding hierarchies) then the slab test is pure multiplication and adds.  The branchless version then replaces all branch tests with min/max operations.  I'd be interested if your div free version is faster and doesn't need the precalculation/storage costs.

Cheers,

Jamie

[xiio...@gmail.com]

On Thursday, 22 September 2016 19:43:35 UTC+1, Jamie Clarkson wrote:

Sure, as you say for your use case it might be best for you, and you might need different information from the test anyway.   Usually in a fast ray tracer you'd pre-calculate the inverse of the direction components and store them (since you do LOTS of ray-bbox intersections for each ray while traversing the bounding hierarchies) then the slab test is pure multiplication and adds.  The branchless version then replaces all branch tests with min/max operations.  I'd be interested if your div free version is faster and doesn't need the precalculation/storage costs.

Cheers,

Jamie

I'm currently fiddling with the "no-division" version - it's not optimised. It doesn't actually avoid the division, unless all that is needed is a true/false (intersects) return - for a function that also returns the scalar along the ray the division is still needed, but it's not in the boundary check 

Easier to show : Original code (in my first post)

front_n = (x_front - s.X) / v.X

ripY = s.Y + front_n*v.Y

ripZ = s.Z + front_n*v.Z

if (front_n > 0) && (ymin <= ripY && ripY <= ymax) && (zmin <= ripZ && ripZ <= zmax) { ..etc ..

becomes  (completely non-optimized)

 

front_n_d = (x_front - s.X) // non-divided version

positive = (v.X > 0)

 

if (front_n_d > ) == positive &&

((ymin-s.Y)*v.X <= front_n_d*v.Y)==positive &&

(front_n_d*v.Y <= (ymax-s.Y)*v.X)==positive &&

((zmin-s.Z)*v.X <= front_n_d*v.Z)==positive &&

(front_n_d*v.Z <= (zmax-s.Z)*v.X))==positive {

..etc..

 

(terms such as ymin-s.Y are reused many times .. I haven't remove them here as it obscures the math.)

 

the obvious alternative is to add an extra if positive .. else .. to simplify all those ==positive out of the condition.. either way I'm not expecting it to be faster, but the numerical errors should be much less, since each side has only 1 multiply and 1 add/subtract, whereas the original form had all the math on one side of the boundary condition - consisting of a subtract followed by a divide ie (x_front - s.X) / v.X (or a multiplication by a reciprocal which adds more potential numerical error), compounded into an additional multiplication and add ie s.Y + front_n*v.Y - so in terms of numerical accuracy the non-multiply version should be much better.

however if the condition is true the divide is still required to get the scalar if needed :

front_n = front_n_d / v.X

So for multiple rays the pre-store of the reciprocal of v.X, v.Y, v.Z still make sense. 

For the usual raytracing bounding boxes a few errors, or slightly oversized bounding boxes are fine, but in my case the ray-box intersection is used to generate indexes for a 3d voxel array - so errors here can throw an index out of range (about 1 in 1 million rays, but still a run time panic for me)