{"id":192,"date":"2009-07-15T13:46:29","date_gmt":"2009-07-15T13:46:29","guid":{"rendered":"https:\/\/blogs.mathworks.com\/loren\/2009\/07\/15\/computational-geometry-in-matlab-r2009a-part-i\/"},"modified":"2016-08-01T10:45:31","modified_gmt":"2016-08-01T15:45:31","slug":"computational-geometry-in-matlab-r2009a-part-i","status":"publish","type":"post","link":"https:\/\/blogs.mathworks.com\/loren\/2009\/07\/15\/computational-geometry-in-matlab-r2009a-part-i\/","title":{"rendered":"Computational Geometry in MATLAB R2009a, Part I"},"content":{"rendered":"<div xmlns:mwsh=\"https:\/\/www.mathworks.com\/namespace\/mcode\/v1\/syntaxhighlight.dtd\" class=\"content\">\r\n   <introduction>\r\n      <p>I'm pleased to introduce Damian Sheehy as this week's guest blogger. Damian is a geometry developer at The MathWorks, his\r\n         background is in the area of geometric modeling and mesh generation.\r\n      <\/p>\r\n   <\/introduction>\r\n   <h3>Contents<\/h3>\r\n   <div>\r\n      <ul>\r\n         <li><a href=\"#1\">A Creature Called Epsilon<\/a><\/li>\r\n         <li><a href=\"#2\">Why Geometric Computing is Susceptible to Robustness Issues<\/a><\/li>\r\n         <li><a href=\"#7\">Robust Delaunay Triangulations<\/a><\/li>\r\n         <li><a href=\"#10\">Collinear Test Revisited<\/a><\/li>\r\n         <li><a href=\"#20\">What Conclusions Can We Draw From This?<\/a><\/li>\r\n         <li><a href=\"#34\">The Bottom Line?<\/a><\/li>\r\n         <li><a href=\"#35\">Coming up Next<\/a><\/li>\r\n      <\/ul>\r\n   <\/div>\r\n   <h3>A Creature Called Epsilon<a name=\"1\"><\/a><\/h3>\r\n   <p>When I set out to write this I reviewed some of Loren's past blogs to refresh my sense of what goes. When I read her blog\r\n      on <a href=\"https:\/\/blogs.mathworks.com\/loren\/2006\/08\/23\/a-glimpse-into-floating-point-accuracy\/\">A Glimpse into Floating-Point Accuracy<\/a>, <a href=\"https:\/\/blogs.mathworks.com\/loren\/2008\/05\/20\/another-lesson-in-floating-point\/\">Another Lesson in Floating Point<\/a>, and <a href=\"https:\/\/blogs.mathworks.com\/loren\/2008\/06\/06\/collinearity\/\">Collinearity<\/a>, I knew I was in familiar territory. The chances are, you too have had a brush with epsilon at some point in your programming\r\n      past.\r\n   <\/p>\r\n   <p>I had my first encounter with <a href=\"http:\/\/en.wikipedia.org\/wiki\/Unit_in_the_last_place\">ulp<\/a> (unit in the last place) and his accomplice <a href=\"https:\/\/www.mathworks.com\/help\/releases\/R2009a\/techdoc\/ref\/eps.html\"><tt>eps<\/tt><\/a> when I was a graduate student in the early 90s. My advisor and I were working on an algorithm that was based on a <a href=\"http:\/\/en.wikipedia.org\/wiki\/Delaunay_triangulation\">Delaunay triangulation<\/a>. The triangulation algorithm would often \"fall over\" due to numerical problems. I wasted effort trying to \"fine tune\" the\r\n      tolerances to get our collection of datasets to triangulate successfully. As you know, this is like trying to push a wrinkle\r\n      out of a fully-fitted carpet. Sure enough, before long my advisor would present a new dataset from a research collaborator\r\n      and I was back in business. The frustrating part was the lack of information on how to deal with these failures, because they\r\n      were rarely addressed in texts or technical publications. When I was finishing up at graduate school and thereafter things\r\n      began to change for the better. Technical journals and conferences dealing with geometric computing began to solicit research\r\n      papers on robustness. The papers that appeared through the years were encouraging and I had a reassuring sense that help was\r\n      on the way.\r\n   <\/p>\r\n   <h3>Why Geometric Computing is Susceptible to Robustness Issues<a name=\"2\"><\/a><\/h3>\r\n   <p>In Loren's blog on <a href=\"https:\/\/blogs.mathworks.com\/loren\/2008\/06\/06\/collinearity\/\">Collinearity<\/a>, we realized the collinearity test for three points could viewed as computing the area of a triangle formed by the points\r\n      and then checking for zero area. This reduced the geometric test to computing a determinant. In fact, this geometric test\r\n      can also be used to tell whether a query point lies to the left or the right of an oriented line defined by two points. Here,\r\n      give it a try . . .\r\n   <\/p><pre style=\"background: #F9F7F3; padding: 10px; border: 1px solid rgb(200,200,200)\">p1 = [1 1];\r\np2 = [5 5];\r\np3 = [2 4];\r\nplotPointLine(p1,p2,p3);\r\naxis <span style=\"color: #A020F0\">equal<\/span>;\r\npointLineOrientation(p1,p2,p3)<\/pre><pre style=\"font-style:oblique\">ans =\r\nLEFT\r\n<\/pre><img decoding=\"async\" vspace=\"5\" hspace=\"5\" src=\"https:\/\/blogs.mathworks.com\/images\/loren\/192\/CgBlogPart1_01.png\"> <pre style=\"background: #F9F7F3; padding: 10px; border: 1px solid rgb(200,200,200)\">q1 = [1 1];\r\nq2 = [5 5];\r\nq3 = [4 2];\r\npointLineOrientation(q1,q2,q3)<\/pre><pre style=\"font-style:oblique\">ans =\r\nRIGHT\r\n<\/pre><pre style=\"background: #F9F7F3; padding: 10px; border: 1px solid rgb(200,200,200)\">dbtype <span style=\"color: #A020F0\">pointLineOrientation<\/span><\/pre><pre style=\"font-style:oblique\">\r\n1     function orient = pointLineOrientation(p1,p2,p3)\r\n2          mat = [p1(1)-p3(1) p1(2)-p3(2); ...\r\n3                 p2(1)-p3(1) p2(2)-p3(2)];\r\n4          detmat = det(mat);   \r\n5          if(detmat == 0)\r\n6              orient = 'ON';\r\n7          else if(detmat &lt; 0)\r\n8                  orient = 'RIGHT';\r\n9          else orient = 'LEFT';\r\n10             end\r\n11         end\r\n12         \r\n\r\n<\/pre><pre style=\"background: #F9F7F3; padding: 10px; border: 1px solid rgb(200,200,200)\">dbtype <span style=\"color: #A020F0\">plotPointLine<\/span><\/pre><pre style=\"font-style:oblique\">\r\n1     function plotPointLine(p1,p2,p3)\r\n2     P = [p1;p2;p3];\r\n3     plot(P(1:2,1),P(1:2,2));\r\n4     hold on;\r\n5     plot(p3(1),p3(2),'*r')\r\n6     ptlabels = arrayfun(@(n) {sprintf('  %d', n)}, (1:3)');\r\n7     Hpl = text(P(:,1),P(:,2),ptlabels, 'FontWeight', 'bold');\r\n\r\n<\/pre><p>If the point lies to the left of the line, the sign of the determinant will be positive &#8211; corresponding to a triangle of positive\r\n      area. Conversely if the point lies on the other side. As we saw before, if the query point is on the line all three points\r\n      are collinear and the determinant is zero. Many problems in geometric computing reduce to evaluating the sign of a determinant.\r\n      The construction of a 2D Delaunay triangulation is based on two geometric tests. The point-line orientation test which we\r\n      just saw, and a point-in-circle test which is used to decide if a query point is within the circumcircle defined by the three\r\n      vertices of a triangle. It is surprising that you only need two simple geometric tests that you learned in high-school and\r\n      some programming knowledge to write a 2D Delaunay triangulation algorithm. However, if you were to try this for real you would\r\n      find things can sometimes go wrong &#8211; very wrong.\r\n   <\/p>\r\n   <h3>Robust Delaunay Triangulations<a name=\"7\"><\/a><\/h3>\r\n   <p>The integrity of the Delaunay triangulation algorithm hinges on the correct evaluation of the determinants used in the geometric\r\n      tests. Problems arise when the determinant is ambiguously close to zero. Inconsistencies creep in due to the presence of numerical\r\n      roundoff, orientation tests can return an incorrect outcome, and the algorithm can subsequently fail.\r\n   <\/p>\r\n   <p>The Qhull computational geometry library, which was developed in the early 1990s, identified these numerical precision problems\r\n      during the computation. Qhull communicated the problems to the user through a warning or error message, and provided some\r\n      helpful interactive options to try to work around the issues.\r\n   <\/p>\r\n   <p>For example, if we compute the Delaunay triangulation of the corner points of a unit square using the Qhull-based delaunayn\r\n      method, a numerical precision warning is issued; this is because all four points are co-circular. Qhull provides an option\r\n      to \"workaround\" these problems; the 'QJ' option adds noise to the points to avoid this degenerate scenario.\r\n   <\/p><pre style=\"background: #F9F7F3; padding: 10px; border: 1px solid rgb(200,200,200)\">X = [0 0; 1 0; 1 1; 0 1]\r\ntri = delaunayn(X, {<span style=\"color: #A020F0\">'QJ'<\/span>})<\/pre><pre style=\"font-style:oblique\">X =\r\n     0     0\r\n     1     0\r\n     1     1\r\n     0     1\r\nWarning: qhull precision warning:\r\nThe initial hull is narrow (cosine of min. angle is 1.0000000000000000).\r\nA coplanar point may lead to a wide facet.  Options 'QbB' (scale to unit box)\r\nor 'Qbb' (scale last coordinate) may remove this warning.  Use 'Pp' to skip\r\nthis warning.  See 'Limitations' in http:\/\/www.qhull.org\/html\/qh-impre.htm\r\n\r\n \r\ntri =\r\n     1     2     4\r\n     4     2     3\r\n<\/pre><p>While this workaround is helpful, is not guaranteed to be robust and failure can still arise in practice.  Fortunately, technology\r\n      has progressed since the early 1990s.\r\n   <\/p>\r\n   <p>In Loren's earlier blog on Collinearity, follow-up comments cited numerical inaccuracies associated with the computation of\r\n      a determinant, and recommended the use of rand\/SVD as being more numerically reliable. But, reducing numerical roundoff does\r\n      not guarantee the correct outcome of the geometric test either. The implementation is still vulnerable to numerical failure.\r\n      In the past decade researchers have focused attention on robustness issues like these and on robustly computing Delaunay triangulations\r\n      in particular. Regarding Delaunay triangulations, the generally accepted solution is based on the notion of Exact Geometric\r\n      Computing (EGC).\r\n   <\/p>\r\n   <p>One reply posted on Loren's <a href=\"https:\/\/blogs.mathworks.com\/loren\/2008\/06\/06\/collinearity\/\">Collinearity<\/a> blog highlighted this solution, but it didn't catch focus. The idea is based on evaluating the determinant using exact arithmetic.\r\n      Since this is computationally expensive, a two-step process is applied; the determinant is computed using floating point arithmetic\r\n      in the usual way and a corresponding error bound is worked out along with this computation. The error bound is used to filter\r\n      out the ambiguous cases, and EGC is then applied to these cases to ensure a correct outcome.\r\n   <\/p>\r\n   <p>In R2009a we adopted 2D and 3D Delaunay triangulations from the Computational Geometry Algorithms Library <a href=\"http:\/\/www.cgal.org\">(CGAL)<\/a> to provide more robust, faster, and memory-efficient solutions in MATLAB. CGAL employs EGC and floating point filters to\r\n      guarantee numerical robustness. Let's see how EGC performs for our collinear test. We will use the new class <a href=\"https:\/\/www.mathworks.com\/help\/releases\/R2009a\/techdoc\/index.html?\/access\/helpdesk\/help\/releases\/R2009a\/techdoc\/ref\/delaunaytri.html\"><tt>DelaunayTri<\/tt><\/a>, which is based on CGAL, to experiment a little.\r\n   <\/p>\r\n   <h3>Collinear Test Revisited<a name=\"10\"><\/a><\/h3>\r\n   <p>Let's take a look at Loren's collinear test again, first choose 3 points we know are collinear. If we attempt to construct\r\n      a Delaunay triangulation using <tt>DelaunayTri<\/tt>, then a triangle should not be formed because the three points are collinear.\r\n   <\/p><pre style=\"background: #F9F7F3; padding: 10px; border: 1px solid rgb(200,200,200)\">format <span style=\"color: #A020F0\">long<\/span>\r\np1 = [1 1]; p2 = [7500 7500]; p3 = [750000 750000];<\/pre><p>Here's the test based on the determinant computation<\/p><pre style=\"background: #F9F7F3; padding: 10px; border: 1px solid rgb(200,200,200)\">collinear = (det([p2-p1 ; p3-p1]) == 0)<\/pre><pre style=\"font-style:oblique\">collinear =\r\n     1\r\n<\/pre><p>and the test based on rank computation.<\/p><pre style=\"background: #F9F7F3; padding: 10px; border: 1px solid rgb(200,200,200)\">collinear = (rank ([p1 ; p2 ; p3]) &lt; 2)<\/pre><pre style=\"font-style:oblique\">collinear =\r\n     1\r\n<\/pre><pre style=\"background: #F9F7F3; padding: 10px; border: 1px solid rgb(200,200,200)\">P = [p1;p2;p3];<\/pre><p>Here's the test based on EGC.<\/p><pre style=\"background: #F9F7F3; padding: 10px; border: 1px solid rgb(200,200,200)\">dt = DelaunayTri(P)<\/pre><pre style=\"font-style:oblique\">dt = \r\n  DelaunayTri\r\n\r\n  Properties:\r\n      Constraints: []\r\n                X: [3x2 double]\r\n    Triangulation: []\r\n<\/pre><p>The Triangulation is empty <tt>[]<\/tt> since the points are collinear, this is what we expect. Now move point <tt>p1<\/tt> by a small amount to make the three points non-collinear.\r\n   <\/p><pre style=\"background: #F9F7F3; padding: 10px; border: 1px solid rgb(200,200,200)\">p1 = [1 1+eps(5)];<\/pre><p>Here's the test based on the determinant computation<\/p><pre style=\"background: #F9F7F3; padding: 10px; border: 1px solid rgb(200,200,200)\">collinear = (det([p2-p1 ; p3-p1]) == 0)\t<span style=\"color: #228B22\">% Gives incorrect result<\/span><\/pre><pre style=\"font-style:oblique\">collinear =\r\n     1\r\n<\/pre><p>and the test based on rank computation.<\/p><pre style=\"background: #F9F7F3; padding: 10px; border: 1px solid rgb(200,200,200)\">collinear = (rank ([p1 ; p2 ; p3]) &lt; 2) <span style=\"color: #228B22\">% Gives incorrect result<\/span><\/pre><pre style=\"font-style:oblique\">collinear =\r\n     1\r\n<\/pre><p>Finally, here's the EGC-based test.<\/p><pre style=\"background: #F9F7F3; padding: 10px; border: 1px solid rgb(200,200,200)\">P = [p1;p2;p3];\r\ndt = DelaunayTri(P);\r\ncollinear = isempty(dt(:,:))            <span style=\"color: #228B22\">% Gives the correct result<\/span><\/pre><pre style=\"font-style:oblique\">collinear =\r\n     0\r\n<\/pre><p>In this instance a triangle was constructed which indicates the three points are not considered to be collinear. We know this\r\n      is the correct assessment.\r\n   <\/p>\r\n   <h3>What Conclusions Can We Draw From This?<a name=\"20\"><\/a><\/h3>\r\n   <p>Well, EGC plays a very important role in the robust implementation of algorithms like Delaunay triangulations, Voronoi diagrams,\r\n      and convex hulls. These algorithms have a long history of being notoriously susceptible to numerical precision issues.\r\n   <\/p>\r\n   <p>This raises an interesting question; should we adopt exact geometric tests for general algorithm development within MATLAB?\r\n      For example; should we be using exact collinear tests and exact point-line orientation tests and the like?\r\n   <\/p>\r\n   <p>In general we should not, because the environment we are working in has finite precision. If we used EGC extensively we would\r\n      realize that while it is numerically correct, it may fail to capture our programming intent. In the flow of an algorithm,\r\n      we take the output from one step of the computation and pass it as input to the next.\r\n   <\/p>\r\n   <p>What happens if we pass inexact input to an exact geometric test? Suppose we begin our computation with three collinear points,\r\n      and next we apply a matrix of rotation to reorient the points in a particular direction. In doing so we introduce inaccuracies;\r\n      will the exact geometric test right those inaccuracies? No, it will not; see for yourself.\r\n   <\/p><pre style=\"background: #F9F7F3; padding: 10px; border: 1px solid rgb(200,200,200)\">p1 = [1 1];\r\np2 = [7500 7500];\r\np3 = [750000 750000];\r\nP = [p1;p2;p3];\r\nT = [cos(pi\/3) sin(pi\/3); -sin(pi\/3) cos(pi\/3)];\r\nP2 = P*T;\r\ndt = DelaunayTri(P2);\r\ncollinear = isempty(dt(:,:))<\/pre><pre style=\"font-style:oblique\">collinear =\r\n     0\r\n<\/pre><p>A triangle was constructed from points {p1, p2, p3}, so in terms of EGC the three points are considered to be noncollinear.<\/p>\r\n   <p>While the outcome is technically correct it's probably not what you expected or hoped for. In the context of Delaunay triangulations\r\n      and their applications, the consequences of imprecise input are relatively benign. EGC is a powerful tool for addressing a\r\n      long-standing robustness problem in this domain.\r\n   <\/p>\r\n   <p>In our programming world of finite precision a judicious choice of tolerance is generally sufficient to capture the correct\r\n      behavior. Now, that just begs the question; what's a judicious choice of tolerance? What would be a suitable tolerance for\r\n      use in the collinearity test, for example? Well, in this context, a good tolerance in one that captures the loss in the computation\r\n      of the determinant in a dynamic or relative sense. It should work when the coordinates of the points are small and it should\r\n      work when they are large. A predefined fixed value will not work over a wide range of coordinate values.  When I compute a\r\n      tolerance-sensitive determinant to be used in a geometric test, I examine the products that are produced during the expansion.\r\n      I then choose the product with the largest magnitude and use it to scale eps.\r\n   <\/p>\r\n   <p>For example, when computing the determinant of matrix <tt>[a b; c d]<\/tt>, I choose the magnitude,\r\n   <\/p><pre>     mag = max(abs(a*d), abs(b*c))<\/pre><p>and then compute the relative tolerance,<\/p><pre>     relativeTol = eps(mag);<\/pre><p>and I apply the tolerance as follows.<\/p><pre>     collinear = (abs(det([p2-p1 ; p3-p1])  &lt;= relativeTol))<\/pre><p>(For simplicity, the products <tt>a*d<\/tt> and <tt>b*c<\/tt> are not reused to compute <tt>det<\/tt>.)\r\n   <\/p>\r\n   <h3>The Bottom Line?<a name=\"34\"><\/a><\/h3>\r\n   <p>So, what's the point of telling you all about Exact Geometric Computing? Ohhh, I am hoping this overview gives you some insight\r\n      into numerical issues in geometric computing. But the important point for you, the user, is you no longer need to worry about\r\n      numerical robustness issues when constructing 2D\/3D Delaunay triangulations, Voronoi diagrams or Convex hulls in MATLAB. Since\r\n      other applications like nearest neighbor queries and scattered data interpolation are built on these features, you can expect\r\n      robstness, performance improvements, and memory efficiency there as well. If you would like to learn more, check out the <tt>video<\/tt> and the Delaunay triangulation demo highlighting the new computational geometry features in R2009a. Follow these links for the reference pages for the new features:\r\n   <\/p>\r\n   <div>\r\n      <ul>\r\n         <li><a href=\"https:\/\/www.mathworks.com\/help\/releases\/R2009a\/techdoc\/index.html?\/access\/helpdesk\/help\/releases\/R2009a\/techdoc\/ref\/trirep.html\"><tt>TriRep<\/tt><\/a>,\r\n         <\/li>\r\n         <li><a href=\"https:\/\/www.mathworks.com\/help\/releases\/R2009a\/techdoc\/index.html?\/access\/helpdesk\/help\/releases\/R2009a\/techdoc\/ref\/delaunaytri.html\"><tt>DelaunayTri<\/tt><\/a>,\r\n         <\/li>\r\n         <li><a href=\"https:\/\/www.mathworks.com\/help\/releases\/R2009a\/techdoc\/index.html?\/access\/helpdesk\/help\/releases\/R2009a\/techdoc\/ref\/triscatteredinterp.html\"><tt>TriScatteredInterp<\/tt><\/a>.\r\n         <\/li>\r\n      <\/ul>\r\n   <\/div>\r\n   <h3>Coming up Next<a name=\"35\"><\/a><\/h3>\r\n   <p>Next week I will talk about an important application area of Delaunay triangulations in MATLAB, and that's Scattered Data\r\n      Interpolation. We introduced new interpolation features in R2009a that improve usability, robustness, and performance, so\r\n      stay tuned for more insight into that. In the meantime, please let me know your views, thoughts, and opinions <a href=\"https:\/\/blogs.mathworks.com\/loren\/?p=192#respond\">here<\/a>.\r\n   <\/p><script language=\"JavaScript\">\r\n<!--\r\n\r\n    function grabCode_46e09184895a4b9db2d4c1ff205ca9b7() {\r\n        \/\/ Remember the title so we can use it in the new page\r\n        title = document.title;\r\n\r\n        \/\/ Break up these strings so that their presence\r\n        \/\/ in the Javascript doesn't mess up the search for\r\n        \/\/ the MATLAB code.\r\n        t1='46e09184895a4b9db2d4c1ff205ca9b7 ' + '##### ' + 'SOURCE BEGIN' + ' #####';\r\n        t2='##### ' + 'SOURCE END' + ' #####' + ' 46e09184895a4b9db2d4c1ff205ca9b7';\r\n    \r\n        b=document.getElementsByTagName('body')[0];\r\n        i1=b.innerHTML.indexOf(t1)+t1.length;\r\n        i2=b.innerHTML.indexOf(t2);\r\n \r\n        code_string = b.innerHTML.substring(i1, i2);\r\n        code_string = code_string.replace(\/REPLACE_WITH_DASH_DASH\/g,'--');\r\n\r\n        \/\/ Use \/x3C\/g instead of the less-than character to avoid errors \r\n        \/\/ in the XML parser.\r\n        \/\/ Use '\\x26#60;' instead of '<' so that the XML parser\r\n        \/\/ doesn't go ahead and substitute the less-than character. \r\n        code_string = code_string.replace(\/\\x3C\/g, '\\x26#60;');\r\n\r\n        author = 'Loren Shure';\r\n        copyright = 'Copyright 2009 The MathWorks, Inc.';\r\n\r\n        w = window.open();\r\n        d = w.document;\r\n        d.write('<pre>\\n');\r\n        d.write(code_string);\r\n\r\n        \/\/ Add author and copyright lines at the bottom if specified.\r\n        if ((author.length > 0) || (copyright.length > 0)) {\r\n            d.writeln('');\r\n            d.writeln('%%');\r\n            if (author.length > 0) {\r\n                d.writeln('% _' + author + '_');\r\n            }\r\n            if (copyright.length > 0) {\r\n                d.writeln('% _' + copyright + '_');\r\n            }\r\n        }\r\n\r\n        d.write('<\/pre>\\n');\r\n      \r\n      d.title = title + ' (MATLAB code)';\r\n      d.close();\r\n      }   \r\n      \r\n-->\r\n<\/script><p style=\"text-align: right; font-size: xx-small; font-weight:lighter;   font-style: italic; color: gray\"><br><a href=\"javascript:grabCode_46e09184895a4b9db2d4c1ff205ca9b7()\"><span style=\"font-size: x-small;        font-style: italic;\">Get \r\n            the MATLAB code \r\n            <noscript>(requires JavaScript)<\/noscript><\/span><\/a><br><br>\r\n      Published with MATLAB&reg; 7.8<br><\/p>\r\n<\/div>\r\n<!--\r\n46e09184895a4b9db2d4c1ff205ca9b7 ##### SOURCE BEGIN #####\r\n%% Computational Geometry in MATLAB R2009a, Part I\r\n% I'm pleased to introduce Damian Sheehy as this week's guest blogger. \r\n% Damian is a geometry developer at The MathWorks, his background is in the \r\n% area of geometric modeling and mesh generation.\r\n\r\n%% A Creature Called Epsilon\r\n% When I set out to write this I reviewed some of Loren's past blogs to refresh \r\n% my sense of what goes. When I read her blog on \r\n% <https:\/\/blogs.mathworks.com\/loren\/2006\/08\/23\/a-glimpse-into-floating-point-accuracy\/ A Glimpse into Floating-Point Accuracy>, \r\n% <https:\/\/blogs.mathworks.com\/loren\/2008\/05\/20\/another-lesson-in-floating-point\/ Another Lesson in Floating Point>, and \r\n% <https:\/\/blogs.mathworks.com\/loren\/2008\/06\/06\/collinearity\/ Collinearity>,\r\n% I knew I was in familiar territory. The chances are, you too have had a \r\n% brush with epsilon at some point in your programming past.\r\n%\r\n% I had my first encounter with \r\n% <http:\/\/en.wikipedia.org\/wiki\/Unit_in_the_last_place ulp> (unit in the last place)\r\n% and his accomplice \r\n% <https:\/\/www.mathworks.com\/help\/releases\/R2009a\/techdoc\/ref\/eps.html |eps|>\r\n% when I was a graduate student in the early 90s. My advisor and I were working \r\n% on an algorithm that was based on a \r\n% <http:\/\/en.wikipedia.org\/wiki\/Delaunay_triangulation Delaunay triangulation>. \r\n% The triangulation algorithm would often \"fall over\" due to numerical problems. \r\n% I wasted effort trying to \"fine tune\" the tolerances to get our collection \r\n% of datasets to triangulate successfully. As you know, this is like trying \r\n% to push a wrinkle out of a fully-fitted carpet. Sure enough, before long \r\n% my advisor would present a new dataset from a research collaborator and I \r\n% was back in business. The frustrating part was the lack of information on \r\n% how to deal with these failures, because they were rarely addressed in \r\n% texts or technical publications. When I was finishing up at graduate school \r\n% and thereafter things began to change for the better. Technical journals \r\n% and conferences dealing with geometric computing began to solicit research \r\n% papers on robustness. The papers that appeared through the years were \r\n% encouraging and I had a reassuring sense that help was on the way.\r\n\r\n%% Why Geometric Computing is Susceptible to Robustness Issues\r\n% In Loren's blog on\r\n% <https:\/\/blogs.mathworks.com\/loren\/2008\/06\/06\/collinearity\/ Collinearity>, \r\n% we realized the collinearity test for three points could viewed as computing \r\n% the area of a triangle formed by the points and then checking for zero \r\n% area. This reduced the geometric test to computing a determinant. In fact, \r\n% this geometric test can also be used to tell whether a query point lies \r\n% to the left or the right of an oriented line defined by two points. Here, \r\n% give it a try . . . \r\np1 = [1 1];\r\np2 = [5 5];\r\np3 = [2 4];\r\nplotPointLine(p1,p2,p3);\r\naxis equal;\r\npointLineOrientation(p1,p2,p3)\r\n%%\r\nq1 = [1 1];\r\nq2 = [5 5];\r\nq3 = [4 2];\r\npointLineOrientation(q1,q2,q3)\r\n%%\r\ndbtype pointLineOrientation\r\n%%\r\ndbtype plotPointLine\r\n%%\r\n% If the point lies to the left of the line, the sign of the determinant will \r\n% be positive \u00e2\u20ac\u201c corresponding to a triangle of positive area. Conversely\r\n% if the point lies on the other side. As we saw before, if the query point \r\n% is on the line all three points are collinear and the determinant is zero. \r\n% Many problems in geometric computing reduce to evaluating the sign of a \r\n% determinant. The construction of a 2D Delaunay triangulation is based on \r\n% two geometric tests. The point-line orientation test which we just saw, \r\n% and a point-in-circle test which is used to decide if a query point is \r\n% within the circumcircle defined by the three vertices of a triangle. \r\n% It is surprising that you only need two simple geometric tests that you \r\n% learned in high-school and some programming knowledge to write a 2D \r\n% Delaunay triangulation algorithm. However, if you were to try this for \r\n% real you would find things can sometimes go wrong \u00e2\u20ac\u201c very wrong.\r\n\r\n%% Robust Delaunay Triangulations\r\n% The integrity of the Delaunay triangulation algorithm hinges on the correct \r\n% evaluation of the determinants used in the geometric tests. Problems arise \r\n% when the determinant is ambiguously close to zero. Inconsistencies creep \r\n% in due to the presence of numerical roundoff, orientation tests can return \r\n% an incorrect outcome, and the algorithm can subsequently fail. \r\n%\r\n% The Qhull computational geometry library, which was developed in the early \r\n% 1990s, identified these numerical precision problems during the computation. \r\n% Qhull communicated the problems to the user through a warning or error message, \r\n% and provided some helpful interactive options to try to work around the \r\n% issues. \r\n%\r\n% For example, if we compute the Delaunay triangulation of the corner \r\n% points of a unit square using the Qhull-based delaunayn method, a numerical \r\n% precision warning is issued; this is because all four points are co-circular. \r\n% Qhull provides an option to \"workaround\" these problems; the 'QJ' option \r\n% adds noise to the points to avoid this degenerate scenario.\r\n\r\nX = [0 0; 1 0; 1 1; 0 1]\r\ntri = delaunayn(X, {'QJ'})\r\n%%\r\n% While this workaround is helpful, is not guaranteed to be robust and \r\n% failure can still arise in practice.  Fortunately, technology has progressed \r\n% since the early 1990s.\r\n%\r\n% In Loren's earlier blog on Collinearity, follow-up comments cited numerical \r\n% inaccuracies associated with the computation of a determinant, and \r\n% recommended the use of rand\/SVD as being more numerically reliable. But, \r\n% reducing numerical roundoff does not guarantee the correct outcome of the \r\n% geometric test either. The implementation is still vulnerable to numerical \r\n% failure. In the past decade researchers have focused attention on robustness \r\n% issues like these and on robustly computing Delaunay triangulations in \r\n% particular. Regarding Delaunay triangulations, the generally accepted \r\n% solution is based on the notion of Exact Geometric Computing (EGC). \r\n%\r\n% One reply posted on Loren's \r\n% <https:\/\/blogs.mathworks.com\/loren\/2008\/06\/06\/collinearity\/ Collinearity>\r\n% blog highlighted this solution, but it didn't catch focus. The idea is based \r\n% on evaluating the determinant using exact arithmetic. Since this is \r\n% computationally expensive, a two-step process is applied; the determinant \r\n% is computed using floating point arithmetic in the usual way and a \r\n% corresponding error bound is worked out along with this computation. The \r\n% error bound is used to filter out the ambiguous cases, and EGC is then \r\n% applied to these cases to ensure a correct outcome. \r\n%%\r\n% In R2009a we adopted 2D and 3D Delaunay triangulations from the Computational \r\n% Geometry Algorithms Library <http:\/\/www.cgal.org (CGAL)> to provide more \r\n% robust, faster, and memory-efficient solutions in MATLAB. CGAL employs \r\n% EGC and floating point filters to guarantee numerical robustness. Let's \r\n% see how EGC performs for our collinear test. \r\n% We will use the new class \r\n% <https:\/\/www.mathworks.com\/help\/releases\/R2009a\/techdoc\/index.html?\/access\/helpdesk\/help\/releases\/R2009a\/techdoc\/ref\/delaunaytri.html |DelaunayTri|>, \r\n% which is based on CGAL, to experiment a little.\r\n\r\n%% Collinear Test Revisited\r\n% Let's take a look at Loren's collinear test again, first choose 3 points \r\n% we know are collinear. If we attempt to construct a Delaunay triangulation \r\n% using |DelaunayTri|, then a triangle should not be formed because the three \r\n% points are collinear.\r\n\r\nformat long\r\np1 = [1 1]; p2 = [7500 7500]; p3 = [750000 750000];\r\n%%\r\n% Here's the test based on the determinant computation\r\ncollinear = (det([p2-p1 ; p3-p1]) == 0)\r\n%%\r\n% and the test based on rank computation.\r\ncollinear = (rank ([p1 ; p2 ; p3]) < 2)\r\n%%\r\nP = [p1;p2;p3];\r\n%%\r\n% Here's the test based on EGC.\r\ndt = DelaunayTri(P)\r\n\r\n%%\r\n% The Triangulation is empty |[]| since the points are collinear, this is what \r\n% we expect. Now move point |p1| by a small amount to make the three points \r\n% non-collinear.\r\np1 = [1 1+eps(5)];\r\n%%\r\n% Here's the test based on the determinant computation\r\ncollinear = (det([p2-p1 ; p3-p1]) == 0)\t% Gives incorrect result\r\n%%\r\n% and the test based on rank computation.\r\ncollinear = (rank ([p1 ; p2 ; p3]) < 2) % Gives incorrect result\r\n%%\r\n% Finally, here's the EGC-based test.\r\nP = [p1;p2;p3];\r\ndt = DelaunayTri(P);\r\ncollinear = isempty(dt(:,:))            % Gives the correct result\r\n%%\r\n% In this instance a triangle was constructed which indicates the three \r\n% points are not considered to be collinear. We know this is the correct \r\n% assessment.\r\n\r\n%% What Conclusions Can We Draw From This?\r\n% Well, EGC plays a very important role in the robust implementation of \r\n% algorithms like Delaunay triangulations, Voronoi diagrams, and convex hulls. \r\n% These algorithms have a long history of being notoriously susceptible to \r\n% numerical precision issues. \r\n%%\r\n% This raises an interesting question; should \r\n% we adopt exact geometric tests for general algorithm development within \r\n% MATLAB? For example; should we be using exact collinear tests and exact \r\n% point-line orientation tests and the like? \r\n%%\r\n% In general we should not, because the environment we are working in has \r\n% finite precision. If we used EGC extensively we would realize that while \r\n% it is numerically correct, it may fail to capture our programming intent. \r\n% In the flow of an algorithm, we take the output from one step of the \r\n% computation and pass it as input to the next. \r\n%%\r\n% What happens if we pass \r\n% inexact input to an exact geometric test? Suppose we begin our computation \r\n% with three collinear points, and next we apply a matrix of rotation to reorient \r\n% the points in a particular direction. In doing so we introduce inaccuracies; \r\n% will the exact geometric test right those inaccuracies? No, it will not;  \r\n% see for yourself.\r\np1 = [1 1];\r\np2 = [7500 7500];\r\np3 = [750000 750000];\r\nP = [p1;p2;p3];\r\nT = [cos(pi\/3) sin(pi\/3); -sin(pi\/3) cos(pi\/3)];\r\nP2 = P*T;\r\ndt = DelaunayTri(P2);\r\ncollinear = isempty(dt(:,:))\r\n%%\r\n% A triangle was constructed from points {p1, p2, p3}, so in terms of EGC\r\n% the three points are considered to be noncollinear.\r\n%%\r\n% While the outcome is technically correct it's probably not what you expected \r\n% or hoped for. \r\n% In the context of Delaunay triangulations and their applications, the \r\n% consequences of imprecise input are relatively benign. EGC is a powerful \r\n% tool for addressing a long-standing robustness problem in this domain.\r\n%%\r\n% In our programming world of finite precision a judicious choice \r\n% of tolerance is generally sufficient to capture the correct behavior. Now, \r\n% that just begs the question; what's a judicious choice of tolerance? \r\n% What would be a suitable tolerance for use in the collinearity test, \r\n% for example? Well, in this context, a good tolerance in one that captures \r\n% the loss in the computation of the determinant in a dynamic or relative \r\n% sense. It should work when the coordinates of the points are small and it \r\n% should work when they are large. A predefined fixed value will not work \r\n% over a wide range of coordinate values.  When I compute a tolerance-sensitive \r\n% determinant to be used in a geometric test, I examine the products that \r\n% are produced during the expansion. I then choose the product with the \r\n% largest magnitude and use it to scale eps. \r\n%%\r\n% For example, when computing the determinant of matrix |[a b; c d]|, I choose \r\n% the magnitude, \r\n%%\r\n%       mag = max(abs(a*d), abs(b*c)) \r\n%%\r\n% and then compute the relative tolerance,\r\n%%\r\n%       relativeTol = eps(mag);\r\n%%\r\n% and I apply the tolerance as follows.\r\n%%\r\n%       collinear = (abs(det([p2-p1 ; p3-p1])  <= relativeTol))\r\n%%\r\n% (For simplicity, the products |a*d| and |b*c| are not reused to compute |det|.)\r\n%% The Bottom Line?\r\n% So, what's the point of telling you all about Exact Geometric Computing?\r\n% Ohhh, I am hoping this overview gives you some insight into numerical issues\r\n% in geometric computing. But the important point for you, the user, is you \r\n% no longer need to worry about numerical robustness issues when constructing \r\n% 2D\/3D Delaunay triangulations, Voronoi diagrams or Convex hulls in MATLAB.\r\n% Since other applications like nearest neighbor queries and scattered \r\n% data interpolation are built on these features, you can expect robstness, \r\n% performance improvements, and memory efficiency there as well.\r\n% If you would like to learn more, check out the \r\n% <https:\/\/www.mathworks.com\/products.htmldemos\/shipping\/matlab\/New-MATLAB-Mathematics-Features-in-R2009a.html |video|>\r\n% and the Delaunay triangulation <https:\/\/www.mathworks.com\/products.htmldemos\/shipping\/matlab\/demoDelaunayTri.html demo> \r\n% highlighting the new computational geometry features in R2009a.\r\n% Follow these links for the reference pages for the new features:\r\n%\r\n% * <https:\/\/www.mathworks.com\/help\/releases\/R2009a\/techdoc\/index.html?\/access\/helpdesk\/help\/releases\/R2009a\/techdoc\/ref\/trirep.html |TriRep|>,\r\n% * <https:\/\/www.mathworks.com\/help\/releases\/R2009a\/techdoc\/index.html?\/access\/helpdesk\/help\/releases\/R2009a\/techdoc\/ref\/delaunaytri.html |DelaunayTri|>,\r\n% * <https:\/\/www.mathworks.com\/help\/releases\/R2009a\/techdoc\/index.html?\/access\/helpdesk\/help\/releases\/R2009a\/techdoc\/ref\/triscatteredinterp.html |TriScatteredInterp|>.\r\n\r\n%% Coming up Next\r\n% Next week I will talk about an important application area of Delaunay \r\n% triangulations in MATLAB, and that's Scattered Data Interpolation. \r\n% We introduced new interpolation features in R2009a that improve usability, \r\n% robustness, and performance, so stay tuned for more insight into that. \r\n% In the meantime, please let me know your views, thoughts, and opinions\r\n% <https:\/\/blogs.mathworks.com\/loren\/?p=192#respond here>.\r\n\r\n##### SOURCE END ##### 46e09184895a4b9db2d4c1ff205ca9b7\r\n-->","protected":false},"excerpt":{"rendered":"<p>\r\n   \r\n      I'm pleased to introduce Damian Sheehy as this week's guest blogger. Damian is a geometry developer at The MathWorks, his\r\n         background is in the area of geometric modeling and... <a class=\"read-more\" href=\"https:\/\/blogs.mathworks.com\/loren\/2009\/07\/15\/computational-geometry-in-matlab-r2009a-part-i\/\">read more >><\/a><\/p>","protected":false},"author":39,"featured_media":0,"comment_status":"open","ping_status":"closed","sticky":false,"template":"","format":"standard","meta":[],"categories":[32,6,11],"tags":[],"_links":{"self":[{"href":"https:\/\/blogs.mathworks.com\/loren\/wp-json\/wp\/v2\/posts\/192"}],"collection":[{"href":"https:\/\/blogs.mathworks.com\/loren\/wp-json\/wp\/v2\/posts"}],"about":[{"href":"https:\/\/blogs.mathworks.com\/loren\/wp-json\/wp\/v2\/types\/post"}],"author":[{"embeddable":true,"href":"https:\/\/blogs.mathworks.com\/loren\/wp-json\/wp\/v2\/users\/39"}],"replies":[{"embeddable":true,"href":"https:\/\/blogs.mathworks.com\/loren\/wp-json\/wp\/v2\/comments?post=192"}],"version-history":[{"count":2,"href":"https:\/\/blogs.mathworks.com\/loren\/wp-json\/wp\/v2\/posts\/192\/revisions"}],"predecessor-version":[{"id":1871,"href":"https:\/\/blogs.mathworks.com\/loren\/wp-json\/wp\/v2\/posts\/192\/revisions\/1871"}],"wp:attachment":[{"href":"https:\/\/blogs.mathworks.com\/loren\/wp-json\/wp\/v2\/media?parent=192"}],"wp:term":[{"taxonomy":"category","embeddable":true,"href":"https:\/\/blogs.mathworks.com\/loren\/wp-json\/wp\/v2\/categories?post=192"},{"taxonomy":"post_tag","embeddable":true,"href":"https:\/\/blogs.mathworks.com\/loren\/wp-json\/wp\/v2\/tags?post=192"}],"curies":[{"name":"wp","href":"https:\/\/api.w.org\/{rel}","templated":true}]}}