{"id":652,"date":"2012-08-22T16:10:33","date_gmt":"2012-08-22T20:10:33","guid":{"rendered":"https:\/\/blogs.mathworks.com\/steve\/?p=652"},"modified":"2019-10-31T14:54:10","modified_gmt":"2019-10-31T18:54:10","slug":"visualizing-the-floating-point-behavior-of-the-point-on-line-test","status":"publish","type":"post","link":"https:\/\/blogs.mathworks.com\/steve\/2012\/08\/22\/visualizing-the-floating-point-behavior-of-the-point-on-line-test\/","title":{"rendered":"Visualizing the floating-point behavior of the point-on-line test"},"content":{"rendered":"

Today's post is about this picture, which is from the paper \"Classroom Examples of Robustness Problems in Geometric Computations,\" by L. Kettner, K. Mehlhorn, S. Pion, S. Schirra, and C. Yap, Computational Geometry<\/i>, vol. 40, n. 1, May 2008, pp. 61-78.<\/p>

\"\" <\/p>

Yes, I know, I was supposed to do another follow-up on my Cody problem of eliminating unnecessary polygon vertices<\/a>. But I really am following up on the Cody problem, if a bit indirectly.<\/p>

I was looking at solution 108897<\/a> by Cris<\/a>, and I paused over these lines:<\/p>

v1 = p2-p1;\r\nv2 = p3-p2;\r\nif<\/span> ( v1(1)*v2(2) == v1(2)*v2(1) ) && any(v1.*v2>0)\r\n<\/pre>
I recognized the expression ( v1(1)*v2(2) == v1(2)*v2(1) )<\/tt> as a way to see if the point p2<\/tt> is on the line segment formed by p1<\/tt> and p3<\/tt>. It reminded me of Loren's 06-Jun-2008 blog post<\/a> on methods for determining collinearity, and it also reminded me of the paper above.<\/p>
I decided to see if I could reproduce the figure in MATLAB.<\/p>
The figure concerns the behavior of floating-point methods for computing the planar orientation predicate<\/i>. Given three points $p$, $q$, and $r$ in a plane, does the point $r$ lie to the left of the line through $p$ and $q$, to the right of that line, or on the line? One way to compute this is via the equation<\/p>
$$O(p,q,r) = \\mbox{sign}((q_x - p_x)(r_y - p_y) - (q_y - p_y)(r_x - p_x))$$<\/p>
The figure in the paper represents an experiment using $q=(12,12)$, $r=(24,24)$, and $256^2$ different locations for the point $p$, all of which are extremely<\/i> close to the point $p=(0.5,0.5)$. More specifically, the figure is a visualization of $O((p_x+xu, p_y+yu),q,r)$ for $0 \\leq x,y \\leq 255$ and $u=2^{-53}$. This value of $u$ is the increment between adjacent double-precision floating-point numbers in the interval $0.5 < x \\leq 1$.<\/p>
Here's a straightforward implementation of the orientation equation above.<\/p>
function<\/span> out = floatOrient(p,q,r)\r\n<\/pre>
out = sign((q(1) - p(1))*(r(2) - p(2)) - ...<\/span>\r\n    (q(2) - p(2))*(r(1) - p(1)));\r\n<\/pre>
Let's see what we get when we compute the orientation for all those different point locations $p'$ that are extremely close to $(0.5,0.5)$.<\/p>
[m,n] = meshgrid(0:255);\r\nPx = 0.5 + n*eps(0.5);\r\nPy = 0.5 + m*eps(0.5);\r\n\r\nq = [12 12];\r\nr = [24 24];\r\n\r\nout = zeros(size(Px));\r\nfor<\/span> u = 1:256\r\n    for<\/span> v = 1:256\r\n        out(u,v) = floatOrient([Px(u,v), Py(u,v)], q, r);\r\n    end<\/span>\r\nend<\/span>\r\n<\/pre>
(Yes, I know ... floatOrient<\/tt> and the calling code above can all be vectorized.)<\/p>
I'll show in yellow the locations $p'$ where $r$ has been determined to be on the line containing $p'$ and $q$. I'll use blue and red to show where $r$ has been determined to lie on one side or the other of the line.<\/p>
image(out + 2);\r\naxis equal<\/span>;\r\nset(gca,'ydir'<\/span>,'normal'<\/span>)\r\naxis off<\/span>\r\ncolormap([0 0 1; 1 1 0; 1 0 0])\r\nhold on<\/span>\r\nplot([.5 256.5],[.5 256.6],'k'<\/span>,'LineWidth'<\/span>,3)\r\nhold off<\/span>\r\n<\/pre>\"\" 
How do we interpret this? First, I want to emphasize that this matrix represents the results for $p'$ that are all extremely<\/i> close to $(0.5,0.5)$. The lower left corner corresponds exactly to the $(0.5,0.5)$, while the upper right corner corresponds to approximately $(0.500000000000028,0.500000000000028)$. Ideally, only the results lying along the black line should<\/i> be yellow; these are the locations where $p'$, $q$, and $r$ are collinear. Our floating-point computation method is showing many other locations in yellow. Notice also that a few pairs of locations $p'$ are actually being computed reversed<\/i>. A left turn has been determined to be a right turn, and vice versa.<\/p>
function<\/span> out = rankCollinear(p,q,r)\r\nout = rank([p-q ; q-r]) < 2;\r\n<\/pre>
Let's see how this test does on our collection of points.<\/p>
for<\/span> u = 1:256\r\n    for<\/span> v = 1:256\r\n        out(u,v) = rankCollinear([Px(u,v), Py(u,v)], q, r);\r\n    end<\/span>\r\nend<\/span>\r\n\r\nimage(out + 1);\r\naxis equal<\/span>;\r\nset(gca,'ydir'<\/span>,'normal'<\/span>)\r\naxis off<\/span>\r\ncolormap([0 0 0; 1 1 0])\r\nhold on<\/span>\r\nplot([.5 256.5],[.5 256.6],'k'<\/span>,'LineWidth'<\/span>,3)\r\nhold off<\/span>\r\n<\/pre>\"\" 
We only see yellow pixels, not red and blue, because rankCollinear<\/tt> only checks for collinearity, not the left or right turning of noncollinear points. The width of the section of yellow pixels depends on the numerical tolerance used by the rank<\/tt> function.<\/p>
I am fascinated by this particular experiment and visualization. It reminds me once again that seemingly straightforward geometric ideas can be devilishly difficult to implement with a computer program.<\/p>
In the comment thread on Loren's collinearity post, Tim Davis suggested<\/a> using a rank test to determine collinearity. It looks something like this:<\/p>