Interested about number precision

XeduR · July 30, 2018, 12:05pm

I’ve just finished working on a function for creating shadows for concave and convex shapes, and I ran into an interesting issue with number precision while working on the function.

I have attached two images to this post. In vertices1.png, the vertices’ coordinates are as is. In vertices2.png, all of the vertices’ coordinates have been run through math.floor().

My function detects intersections and then acts accordingly in order to create the desired shape for the shadow. In vertices1.png, where the coordinates are not run through the floor function, the function states that sides 8_9 and 14_1 intersect, but visually it clear as day that the sides are far apart. In fact, printing out the coordinates for vertex 9 and 14 state that the vertices are over 50 pixels apart. However, if the coordinates are first run through the floor function, then the function does not mistakenly believe that the lines intersect.

All I understand is that this has something to do with number precision in lua. Does someone know more about this topic and perhaps enlighten me as to why non-floored or non-rounded numbers may lead to such drastic inaccuracies?

remiduchalard · July 30, 2018, 1:00pm

I think your unit is point 5 to point 4 is equal at 1.

Can you try to use more precise value?

MyValue=math.floor(10\*myNotRoundValue)\*0.1

XeduR · July 30, 2018, 1:33pm

The current method that I already use in my code, for precision as you mentioned, is

-- method #1 (in use) comparisonVertex[#comparisonVertex+1] = mathRound(vertex[#vertex-1]\*10)\*0.1 comparisonVertex[#comparisonVertex+1] = mathRound(vertex[#vertex]\*10)\*0.1

which yields sufficiently accurate values for x and y coordinates. The thing that I am, however, interested in this particular issue is that unless the values are run through the round function (method 1) or the floor function (method 2), then the values are occasionally unusable due to some lua or floating point related issue.

-- method 2 (not in use) comparisonVertex[#comparisonVertex+1] = mathFloor(vertex[#vertex-1]) comparisonVertex[#comparisonVertex+1] = mathFloor(vertex[#vertex])

The function itself works great and as intended, but I am puzzled by this number precision issue.

For instance, in the vertices1 image above, if I leave the vertices’ coordinates as is and do not round them or floor them, and then I run my intersect function for sides 8_9 and 14_1, then the function states that they intersect. However, if I print out the values and then copy the printed values from the console and into the intersect function, then the function states that they no longer intersect. This most likely points to some floating point number issue, but I am no expert on the matter and I am interested in knowing how can inaccuracies like this cause such problems? I do understand that floating points can be somewhat inaccurate, but confusing coordinates over 50 pixels apart seems absurd.

davebollinger · July 30, 2018, 3:41pm

post the source of your “line intersect” routine, along with your rounding function and the (full, unrounded) coordinates of two line segments that exhibit the problem (fail when unrounded, succeed when rounded). structure it as a stand-alone demo so that it will run unmodified as is, for ex something like:

-- (nonsense values, give real ones) x1, y1 = 1234.5678, 1234.5678 x2, y2 = 1234.5678, 1234.5678 x3, y3 = 1234.5678, 1234.5678 x4, y4 = 1234.5678, 1234.5678 -- expect failure: unrounded: print(intersect(x1,y1,x2,y2,x3,y3,x4,y4)) -- expect success, rounded: print(intersect(round(x1),round(y1),round(x2),round(y2),round(x3),round(y3),round(x4),round(y4))

odds are that it’s not a true “precision” issue, unless testing ==0 instead of abs()<epsilon fe, but some “weakness” in your intersect routine, perhaps when given colinear segments (or some other problematic case). but there’s no way anyone will be able to diagnose that for you without an actual example that replicates the problem.

XeduR · July 30, 2018, 6:06pm

I’ll be posting some code here by tomorrow so that you can replicate the issue as well.

XeduR · July 31, 2018, 11:11am

Here’s the code for the specific scenario. It’ll run as is and has no dependencies.

-- create new vertices for a shape after it is rotated local function rotatePolygon(vertices,r) -- matrixRotate returns rotation matrix for angle theta local function matrixRotate(theta) local t = setmetatable({},mt) t[1] = {math.cos(theta),math.sin(theta)} t[2] = {-math.sin(theta),math.cos(theta)} return t end -- multiply two matrices and return the resulting matrix local function matrixMultiply( m1, m2 ) -- matrixMultiply source: https://github.com/davidm/lua-matrix/blob/master/lua/matrix.lua if #m1[1] ~= #m2 then return nil end local res = {} for i = 1, #m1 do res[i] = {} for j = 1, #m2[1] do res[i][j] = 0 for k = 1, #m2 do res[i][j] = res[i][j] + m1[i][k] \* m2[k][j] end end end return res end local matrix = {} -- turn the vertices table into a matrix for i = 1, #vertices\*0.5 do matrix[i] = {} matrix[i][1] = vertices[i\*2-1] matrix[i][2] = vertices[i\*2] end -- create a rotation matrix for the given rotation local matrixRotation = matrixRotate(math.rad(r)) -- multiply the vertices matrix with the rotation metrix local matrixFinal = matrixMultiply(matrix,matrixRotation) -- update the original vertices and turn the matrix back into a table for i = 1, #matrixFinal do vertices[i\*2-1] = matrixFinal[i][1] vertices[i\*2] = matrixFinal[i][2] end return vertices end -- doLinesIntersect source: https://gist.github.com/HoraceBury/9431861 local function doLinesIntersect( a, b, c, d ) -- parameter conversion local L1 = {X1=a.x,Y1=a.y,X2=b.x,Y2=b.y} local L2 = {X1=c.x,Y1=c.y,X2=d.x,Y2=d.y} -- Denominator for ua and ub are the same, so store this calculation local \_d = (L2.Y2 - L2.Y1) \* (L1.X2 - L1.X1) - (L2.X2 - L2.X1) \* (L1.Y2 - L1.Y1) -- Make sure there is not a division by zero - this also indicates that the lines are parallel. -- If n\_a and n\_b were both equal to zero the lines would be on top of each -- other (coincidental). This check is not done because it is not -- necessary for this implementation (the parallel check accounts for this). if (\_d == 0) then return false end -- n\_a and n\_b are calculated as seperate values for readability local n\_a = (L2.X2 - L2.X1) \* (L1.Y1 - L2.Y1) - (L2.Y2 - L2.Y1) \* (L1.X1 - L2.X1) local n\_b = (L1.X2 - L1.X1) \* (L1.Y1 - L2.Y1) - (L1.Y2 - L1.Y1) \* (L1.X1 - L2.X1) -- Calculate the intermediate fractional point that the lines potentially intersect. local ua = n\_a / \_d local ub = n\_b / \_d -- The fractional point will be between 0 and 1 inclusive if the lines -- intersect. If the fractional calculation is larger than 1 or smaller -- than 0 the lines would need to be longer to intersect. if (ua \>= 0 and ua \<= 1 and ub \>= 0 and ub \<= 1) then local x = L1.X1 + (ua \* (L1.X2 - L1.X1)) local y = L1.Y1 + (ua \* (L1.Y2 - L1.Y1)) return {x=x, y=y} end return false end local light = display.newCircle(0,0,4) light:setFillColor(1,1,0) light.brightness = 1 local shadow = display.newPolygon( 222, 186, {-86,-25,-65,-46,-8,-49,50,6,21,37,-41,102,-63,81,-21,38,-36,23,-64,39,-85,18,-64,-4} ) shadow:setFillColor(1,0,0,0.5) local vertices = {-40,-30,-20,-30,-20,10,0,10,0,-10,20,-10,20,10,40,10,40,30,-40,30} local object = display.newPolygon( 240, 160, vertices ) object.rotation = 224 -- light coordinates 1: five decimals, which ultimately lead to intersect light.x = 346.34442 light.y = 104.40447 -- light coordinates 2: four decimals, which no longer lead to intersect -- light.x = 346.3444 -- light.y = 104.4044 local rotatedVertices = rotatePolygon(vertices,object.rotation) local function newVertex(vX,vY) local shadowLength = math.sqrt((vX-light.x)^2+(vY-light.y)^2)\*(0.5\*light.brightness) local shadowAngle = math.deg(math.atan2(vY-light.y,vX-light.x)) local dx = math.sin(math.rad(shadowAngle-90))\*shadowLength local dy = math.cos(math.rad(shadowAngle-90))\*shadowLength vX = vX-dx vY = vY+dy return vX, vY end local vertex = {} vertex[#vertex+1] = {} vertex[#vertex].x, vertex[#vertex].y = newVertex(rotatedVertices[5]+object.x, rotatedVertices[6]+object.y) vertex[#vertex+1] = {} vertex[#vertex].x, vertex[#vertex].y = newVertex(rotatedVertices[7]+object.x, rotatedVertices[8]+object.y) vertex[#vertex+1] = {} vertex[#vertex].x, vertex[#vertex].y = newVertex(rotatedVertices[13]+object.x, rotatedVertices[14]+object.y) vertex[#vertex+1] = {} vertex[#vertex].x, vertex[#vertex].y = newVertex(rotatedVertices[15]+object.x, rotatedVertices[16]+object.y) local line12 = display.newLine(vertex[1].x,vertex[1].y,vertex[2].x,vertex[2].y) local line34 = display.newLine(vertex[3].x,vertex[3].y,vertex[4].x,vertex[4].y) local number = {} for k = 1, 4 do number[k] = display.newText( k, vertex[k].x,vertex[k].y, native.systemFont, 12 ) number[k]:setFillColor(1,1,0) end -- method 1: using the actual, non-rounded numbers --------------------------------------------------- local intersect = doLinesIntersect( vertex[1], vertex[2], vertex[3], vertex[4] ) -- method 2: print the values and use them directly --------------------------------------------------- -- print(vertex[1].x,vertex[1].y,vertex[2].x,vertex[2].y,vertex[3].x,vertex[3].y,vertex[4].x,vertex[4].y) -- vertex[1].x = 218.82785887369 -- vertex[1].y = 197.84741793315 -- vertex[2].x = 197.24766486353 -- vertex[2].y = 177.00766681938 -- vertex[3].x = 175.66747085337 -- vertex[3].y = 156.16791570561 -- vertex[4].x = 154.08727684321 -- vertex[4].y = 135.32816459184 -- method 3: use floored (or rounded) values --------------------------------------------------- -- for i = 1, 4 do -- vertex[i].x = math.floor(vertex[i].x) -- vertex[i].y = math.floor(vertex[i].y) -- end local intersect = doLinesIntersect( vertex[1], vertex[2], vertex[3], vertex[4] ) print("---") if intersect == false then print("The lines do not intersect.") else print("The lines intersect at x: "..intersect.x..", y: "..intersect.y) end local distance = math.sqrt((vertex[2].x-vertex[3].x)^2+(vertex[2].y-vertex[3].y)^2) print("The distance between vertex 2 and 3 is "..distance.."px.")

I replaced my own intersect function with Horace Bury’s that I found online just to ensure that his function returns the same results as mine, which it does.

I reiterate, the reason why I believe that this has to do with number precision in lua is because,

if the light’s coordinates are given with four decimal points instead of five, then the intersect no longer occurs.
if the four vertices are not altered, then the intersect occurs (as long as light uses 5 decimals), but
if the four vertices are rounded, floored or printed out and inserted with 11 decimal point accuracy, then the intersect no longer occurs.

I am interested in why this kind behaviour of happening. This isn’t really a debugging issue.

davebollinger · July 31, 2018, 2:12pm

as i suspected (re-read my prior post for the “clues” :D)…

what you have here is a classic problem with floating-point equality testing (google it). given the trig involved in coordinate creation, you have very “detailed” numbers, so expecting any two numbers to ever “exactly equal” each other is problematic.

the issue is line 70 where the determinant (aka “denominator” in your code) is compared to exactly zero. it could be addressed by substituting something like the following:

-- if (\_d == 0) then local epsilon = 1e-12 if (math.abs(\_d) \< epsilon) then

this will have a similar effect to rounding the coordinates themselves, but targets the actual location where the problem manifests, rather than brute-forcing all inputs to lower precision. (so you’ll be able to use a much smaller epsilon here)

while it’s true that it’s a precision issue at the core, it’s really a usage issue that causes the problem. i wouldn’t fault your implementation - probably 90% of the line intersection code out there on the internet has the same issue, but they’re all flawed. even “octuple-precision” floats (if it existed) wouldn’t be enough to resolve it, you’d literally need infinite precision - so it’s in knowing how to work with the limited precision.

it’s sort of like blaming a cheeseburger for being overweight, when it’s really the eating of the cheeseburger that needs to be addressed.

hope that helps

XeduR · July 31, 2018, 2:22pm

This definitely helped. Thank you for taking the time to explain the issue.

remiduchalard · July 30, 2018, 1:00pm

I think your unit is point 5 to point 4 is equal at 1.

Can you try to use more precise value?

MyValue=math.floor(10\*myNotRoundValue)\*0.1

XeduR · July 30, 2018, 1:33pm

The current method that I already use in my code, for precision as you mentioned, is

-- method #1 (in use) comparisonVertex[#comparisonVertex+1] = mathRound(vertex[#vertex-1]\*10)\*0.1 comparisonVertex[#comparisonVertex+1] = mathRound(vertex[#vertex]\*10)\*0.1

which yields sufficiently accurate values for x and y coordinates. The thing that I am, however, interested in this particular issue is that unless the values are run through the round function (method 1) or the floor function (method 2), then the values are occasionally unusable due to some lua or floating point related issue.

-- method 2 (not in use) comparisonVertex[#comparisonVertex+1] = mathFloor(vertex[#vertex-1]) comparisonVertex[#comparisonVertex+1] = mathFloor(vertex[#vertex])

The function itself works great and as intended, but I am puzzled by this number precision issue.

For instance, in the vertices1 image above, if I leave the vertices’ coordinates as is and do not round them or floor them, and then I run my intersect function for sides 8_9 and 14_1, then the function states that they intersect. However, if I print out the values and then copy the printed values from the console and into the intersect function, then the function states that they no longer intersect. This most likely points to some floating point number issue, but I am no expert on the matter and I am interested in knowing how can inaccuracies like this cause such problems? I do understand that floating points can be somewhat inaccurate, but confusing coordinates over 50 pixels apart seems absurd.

davebollinger · July 30, 2018, 3:41pm

post the source of your “line intersect” routine, along with your rounding function and the (full, unrounded) coordinates of two line segments that exhibit the problem (fail when unrounded, succeed when rounded). structure it as a stand-alone demo so that it will run unmodified as is, for ex something like:

-- (nonsense values, give real ones) x1, y1 = 1234.5678, 1234.5678 x2, y2 = 1234.5678, 1234.5678 x3, y3 = 1234.5678, 1234.5678 x4, y4 = 1234.5678, 1234.5678 -- expect failure: unrounded: print(intersect(x1,y1,x2,y2,x3,y3,x4,y4)) -- expect success, rounded: print(intersect(round(x1),round(y1),round(x2),round(y2),round(x3),round(y3),round(x4),round(y4))

odds are that it’s not a true “precision” issue, unless testing ==0 instead of abs()<epsilon fe, but some “weakness” in your intersect routine, perhaps when given colinear segments (or some other problematic case). but there’s no way anyone will be able to diagnose that for you without an actual example that replicates the problem.

XeduR · July 30, 2018, 6:06pm

I’ll be posting some code here by tomorrow so that you can replicate the issue as well.

XeduR · July 31, 2018, 11:11am

Here’s the code for the specific scenario. It’ll run as is and has no dependencies.

-- create new vertices for a shape after it is rotated local function rotatePolygon(vertices,r) -- matrixRotate returns rotation matrix for angle theta local function matrixRotate(theta) local t = setmetatable({},mt) t[1] = {math.cos(theta),math.sin(theta)} t[2] = {-math.sin(theta),math.cos(theta)} return t end -- multiply two matrices and return the resulting matrix local function matrixMultiply( m1, m2 ) -- matrixMultiply source: https://github.com/davidm/lua-matrix/blob/master/lua/matrix.lua if #m1[1] ~= #m2 then return nil end local res = {} for i = 1, #m1 do res[i] = {} for j = 1, #m2[1] do res[i][j] = 0 for k = 1, #m2 do res[i][j] = res[i][j] + m1[i][k] \* m2[k][j] end end end return res end local matrix = {} -- turn the vertices table into a matrix for i = 1, #vertices\*0.5 do matrix[i] = {} matrix[i][1] = vertices[i\*2-1] matrix[i][2] = vertices[i\*2] end -- create a rotation matrix for the given rotation local matrixRotation = matrixRotate(math.rad(r)) -- multiply the vertices matrix with the rotation metrix local matrixFinal = matrixMultiply(matrix,matrixRotation) -- update the original vertices and turn the matrix back into a table for i = 1, #matrixFinal do vertices[i\*2-1] = matrixFinal[i][1] vertices[i\*2] = matrixFinal[i][2] end return vertices end -- doLinesIntersect source: https://gist.github.com/HoraceBury/9431861 local function doLinesIntersect( a, b, c, d ) -- parameter conversion local L1 = {X1=a.x,Y1=a.y,X2=b.x,Y2=b.y} local L2 = {X1=c.x,Y1=c.y,X2=d.x,Y2=d.y} -- Denominator for ua and ub are the same, so store this calculation local \_d = (L2.Y2 - L2.Y1) \* (L1.X2 - L1.X1) - (L2.X2 - L2.X1) \* (L1.Y2 - L1.Y1) -- Make sure there is not a division by zero - this also indicates that the lines are parallel. -- If n\_a and n\_b were both equal to zero the lines would be on top of each -- other (coincidental). This check is not done because it is not -- necessary for this implementation (the parallel check accounts for this). if (\_d == 0) then return false end -- n\_a and n\_b are calculated as seperate values for readability local n\_a = (L2.X2 - L2.X1) \* (L1.Y1 - L2.Y1) - (L2.Y2 - L2.Y1) \* (L1.X1 - L2.X1) local n\_b = (L1.X2 - L1.X1) \* (L1.Y1 - L2.Y1) - (L1.Y2 - L1.Y1) \* (L1.X1 - L2.X1) -- Calculate the intermediate fractional point that the lines potentially intersect. local ua = n\_a / \_d local ub = n\_b / \_d -- The fractional point will be between 0 and 1 inclusive if the lines -- intersect. If the fractional calculation is larger than 1 or smaller -- than 0 the lines would need to be longer to intersect. if (ua \>= 0 and ua \<= 1 and ub \>= 0 and ub \<= 1) then local x = L1.X1 + (ua \* (L1.X2 - L1.X1)) local y = L1.Y1 + (ua \* (L1.Y2 - L1.Y1)) return {x=x, y=y} end return false end local light = display.newCircle(0,0,4) light:setFillColor(1,1,0) light.brightness = 1 local shadow = display.newPolygon( 222, 186, {-86,-25,-65,-46,-8,-49,50,6,21,37,-41,102,-63,81,-21,38,-36,23,-64,39,-85,18,-64,-4} ) shadow:setFillColor(1,0,0,0.5) local vertices = {-40,-30,-20,-30,-20,10,0,10,0,-10,20,-10,20,10,40,10,40,30,-40,30} local object = display.newPolygon( 240, 160, vertices ) object.rotation = 224 -- light coordinates 1: five decimals, which ultimately lead to intersect light.x = 346.34442 light.y = 104.40447 -- light coordinates 2: four decimals, which no longer lead to intersect -- light.x = 346.3444 -- light.y = 104.4044 local rotatedVertices = rotatePolygon(vertices,object.rotation) local function newVertex(vX,vY) local shadowLength = math.sqrt((vX-light.x)^2+(vY-light.y)^2)\*(0.5\*light.brightness) local shadowAngle = math.deg(math.atan2(vY-light.y,vX-light.x)) local dx = math.sin(math.rad(shadowAngle-90))\*shadowLength local dy = math.cos(math.rad(shadowAngle-90))\*shadowLength vX = vX-dx vY = vY+dy return vX, vY end local vertex = {} vertex[#vertex+1] = {} vertex[#vertex].x, vertex[#vertex].y = newVertex(rotatedVertices[5]+object.x, rotatedVertices[6]+object.y) vertex[#vertex+1] = {} vertex[#vertex].x, vertex[#vertex].y = newVertex(rotatedVertices[7]+object.x, rotatedVertices[8]+object.y) vertex[#vertex+1] = {} vertex[#vertex].x, vertex[#vertex].y = newVertex(rotatedVertices[13]+object.x, rotatedVertices[14]+object.y) vertex[#vertex+1] = {} vertex[#vertex].x, vertex[#vertex].y = newVertex(rotatedVertices[15]+object.x, rotatedVertices[16]+object.y) local line12 = display.newLine(vertex[1].x,vertex[1].y,vertex[2].x,vertex[2].y) local line34 = display.newLine(vertex[3].x,vertex[3].y,vertex[4].x,vertex[4].y) local number = {} for k = 1, 4 do number[k] = display.newText( k, vertex[k].x,vertex[k].y, native.systemFont, 12 ) number[k]:setFillColor(1,1,0) end -- method 1: using the actual, non-rounded numbers --------------------------------------------------- local intersect = doLinesIntersect( vertex[1], vertex[2], vertex[3], vertex[4] ) -- method 2: print the values and use them directly --------------------------------------------------- -- print(vertex[1].x,vertex[1].y,vertex[2].x,vertex[2].y,vertex[3].x,vertex[3].y,vertex[4].x,vertex[4].y) -- vertex[1].x = 218.82785887369 -- vertex[1].y = 197.84741793315 -- vertex[2].x = 197.24766486353 -- vertex[2].y = 177.00766681938 -- vertex[3].x = 175.66747085337 -- vertex[3].y = 156.16791570561 -- vertex[4].x = 154.08727684321 -- vertex[4].y = 135.32816459184 -- method 3: use floored (or rounded) values --------------------------------------------------- -- for i = 1, 4 do -- vertex[i].x = math.floor(vertex[i].x) -- vertex[i].y = math.floor(vertex[i].y) -- end local intersect = doLinesIntersect( vertex[1], vertex[2], vertex[3], vertex[4] ) print("---") if intersect == false then print("The lines do not intersect.") else print("The lines intersect at x: "..intersect.x..", y: "..intersect.y) end local distance = math.sqrt((vertex[2].x-vertex[3].x)^2+(vertex[2].y-vertex[3].y)^2) print("The distance between vertex 2 and 3 is "..distance.."px.")

I replaced my own intersect function with Horace Bury’s that I found online just to ensure that his function returns the same results as mine, which it does.

I reiterate, the reason why I believe that this has to do with number precision in lua is because,

if the light’s coordinates are given with four decimal points instead of five, then the intersect no longer occurs.
if the four vertices are not altered, then the intersect occurs (as long as light uses 5 decimals), but
if the four vertices are rounded, floored or printed out and inserted with 11 decimal point accuracy, then the intersect no longer occurs.

I am interested in why this kind behaviour of happening. This isn’t really a debugging issue.

davebollinger · July 31, 2018, 2:12pm

as i suspected (re-read my prior post for the “clues” :D)…

what you have here is a classic problem with floating-point equality testing (google it). given the trig involved in coordinate creation, you have very “detailed” numbers, so expecting any two numbers to ever “exactly equal” each other is problematic.

the issue is line 70 where the determinant (aka “denominator” in your code) is compared to exactly zero. it could be addressed by substituting something like the following:

-- if (\_d == 0) then local epsilon = 1e-12 if (math.abs(\_d) \< epsilon) then

this will have a similar effect to rounding the coordinates themselves, but targets the actual location where the problem manifests, rather than brute-forcing all inputs to lower precision. (so you’ll be able to use a much smaller epsilon here)

while it’s true that it’s a precision issue at the core, it’s really a usage issue that causes the problem. i wouldn’t fault your implementation - probably 90% of the line intersection code out there on the internet has the same issue, but they’re all flawed. even “octuple-precision” floats (if it existed) wouldn’t be enough to resolve it, you’d literally need infinite precision - so it’s in knowing how to work with the limited precision.

it’s sort of like blaming a cheeseburger for being overweight, when it’s really the eating of the cheeseburger that needs to be addressed.

hope that helps

XeduR · July 31, 2018, 2:22pm

This definitely helped. Thank you for taking the time to explain the issue.