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, vol. 40, n. 1, May 2008, pp. 61-78.

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

I was looking at solution 108897 by Cris, and I paused over these lines:

v1 = p2-p1;
v2 = p3-p2;
if ( v1(1)*v2(2) == v1(2)*v2(1) ) && any(v1.*v2>0)

I recognized the expression ( v1(1)*v2(2) == v1(2)*v2(1) ) as a way to see if the point p2 is on the line segment formed by p1 and p3. It reminded me of Loren's 06-Jun-2008 blog post on methods for determining collinearity, and it also reminded me of the paper above.

I decided to see if I could reproduce the figure in MATLAB.

The figure concerns the behavior of floating-point methods for computing the planar orientation predicate. 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

$$O(p,q,r) = \mbox{sign}((q_x - p_x)(r_y - p_y) - (q_y - p_y)(r_x - p_x))$$

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 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$.

Here's a straightforward implementation of the orientation equation above.

function out = floatOrient(p,q,r)

out = sign((q(1) - p(1))*(r(2) - p(2)) - ...
    (q(2) - p(2))*(r(1) - p(1)));

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)$.

[m,n] = meshgrid(0:255);
Px = 0.5 + n*eps(0.5);
Py = 0.5 + m*eps(0.5);

q = [12 12];
r = [24 24];

out = zeros(size(Px));
for u = 1:256
    for v = 1:256
        out(u,v) = floatOrient([Px(u,v), Py(u,v)], q, r);
    end
end

(Yes, I know ... floatOrient and the calling code above can all be vectorized.)

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.

image(out + 2);
axis equal;
set(gca,'ydir','normal')
axis off
colormap([0 0 1; 1 1 0; 1 0 0])
hold on
plot([.5 256.5],[.5 256.6],'k','LineWidth',3)
hold off

How do we interpret this? First, I want to emphasize that this matrix represents the results for $p'$ that are all extremely 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 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. A left turn has been determined to be a right turn, and vice versa.

In the comment thread on Loren's collinearity post, Tim Davis suggested using a rank test to determine collinearity. It looks something like this:

function out = rankCollinear(p,q,r)
out = rank([p-q ; q-r]) < 2;

Let's see how this test does on our collection of points.

for u = 1:256
    for v = 1:256
        out(u,v) = rankCollinear([Px(u,v), Py(u,v)], q, r);
    end
end

image(out + 1);
axis equal;
set(gca,'ydir','normal')
axis off
colormap([0 0 0; 1 1 0])
hold on
plot([.5 256.5],[.5 256.6],'k','LineWidth',3)
hold off

We only see yellow pixels, not red and blue, because rankCollinear 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 function.

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.

Get the MATLAB code

Published with MATLAB® 7.14

Published with MATLAB® 7.14