{"id":1190,"date":"2015-07-15T13:47:08","date_gmt":"2015-07-15T18:47:08","guid":{"rendered":"https:\/\/blogs.mathworks.com\/loren\/?p=1190"},"modified":"2017-09-24T19:34:17","modified_gmt":"2017-09-25T00:34:17","slug":"incremental-delaunay-construction","status":"publish","type":"post","link":"https:\/\/blogs.mathworks.com\/loren\/2015\/07\/15\/incremental-delaunay-construction\/","title":{"rendered":"Incremental Delaunay Construction"},"content":{"rendered":"<div class=\"content\"><!--introduction--><p>I'm happy to welcome back Damian Sheehy as guest blogger. Last time Damian wrote about how Natural Neighbor interpolation addresses FAQs in scattered data interpolation. In this blog he will answer a FAQ on adaptively editing a Delaunay triangulation.<\/p><!--\/introduction--><h3>Contents<\/h3><div><ul><li><a href=\"#1cacaed1-f06b-4c25-9174-610bedc44a94\">Is <tt>delaunayTriangulation<\/tt> More Efficient than <tt>delaunay<\/tt>?<\/a><\/li><li><a href=\"#46f99a36-fd9f-42de-8ea0-badc301144c4\">When is Incremental Delaunay Important?<\/a><\/li><li><a href=\"#c06cefa9-4e4a-4fdb-a61f-e73d41e2f270\">Performance Example of Incremental Delaunay Construction<\/a><\/li><li><a href=\"#6a4ed5b2-13cd-4789-abf5-92af268621e8\">Your Need for Geometric Tools?<\/a><\/li><\/ul><\/div><h4>Is <tt>delaunayTriangulation<\/tt> More Efficient than <tt>delaunay<\/tt>?<a name=\"1cacaed1-f06b-4c25-9174-610bedc44a94\"><\/a><\/h4><p>A technical support question that occasionally crops up asks about the best and most efficient way to construct a Delaunay triangulation by adding points to an existing triangulation. We call this operation incremental editing as the goal is to add or remove points in an incremental manner as opposed to recreating from scratch; for example calling <a href=\"https:\/\/www.mathworks.com\/help\/matlab\/ref\/delaunay.html\"><tt>delaunay<\/tt><\/a> multiple times. The documentation for the <tt>delaunayTriangulation<\/tt> class provides examples that show the syntax that allows you to edit a Delaunay triangulation by adding or removing points. Let's look at a simple 2-D example to highlight the concept.<\/p><p>In this example we will load <tt>trimesh2d.mat<\/tt> which ships with MATLAB. The file contains points with coordinates (x, y) and triangulation edge constraints that define the boundary of a domain. Let's triangulate the data and take a look at that.<\/p><pre class=\"codeinput\">load <span class=\"string\">trimesh2d<\/span>\r\ndt = delaunayTriangulation(x,y,Constraints);\r\ntin = find(dt.isInterior());\r\ntriplot(dt(tin,:),dt.Points(:,1),dt.Points(:,2))\r\naxis <span class=\"string\">equal<\/span>\r\n<\/pre><img decoding=\"async\" vspace=\"5\" hspace=\"5\" src=\"https:\/\/blogs.mathworks.com\/images\/loren\/2015\/IncrementalDelaunayBlog2_01.png\" alt=\"\"> <p>Suppose we want to add the circumcenter points to this triangulation. The circumcenter of a triangle is the center of a circumscribed circle that passes through the vertices of the triangle. The <tt>delaunayTriangulation<\/tt> class provides a method to compute them - <a href=\"https:\/\/www.mathworks.com\/help\/matlab\/ref\/triangulation.circumcenter.html\"><tt>delaunayTriangulation\/circumcenter<\/tt><\/a>. Compute them using this method and add them to the plot. Note, some triangles may share the same circumcenter, so I will call the <a href=\"https:\/\/www.mathworks.com\/help\/matlab\/ref\/uniquetol.html\"><tt>uniquetol<\/tt><\/a> function to eliminate the near-duplicates.<\/p><pre class=\"codeinput\">cc = dt.circumcenter(tin);\r\ncc = uniquetol(cc, <span class=\"string\">'ByRows'<\/span>,true);\r\nhold <span class=\"string\">on<\/span>\r\nplot(cc(:,1),cc(:,2),<span class=\"string\">'.r'<\/span>)\r\nhold <span class=\"string\">off<\/span>\r\n<\/pre><img decoding=\"async\" vspace=\"5\" hspace=\"5\" src=\"https:\/\/blogs.mathworks.com\/images\/loren\/2015\/IncrementalDelaunayBlog2_02.png\" alt=\"\"> <p>Now, I will add the circumcenter points to the triangulation and plot the result in a new figure.<\/p><pre class=\"codeinput\">numcc = size(cc,1);\r\ndt.Points(end+(1:numcc),:) = cc;\r\ntin = find(dt.isInterior());\r\nfigure\r\ntriplot(dt(tin,:),dt.Points(:,1),dt.Points(:,2))\r\naxis <span class=\"string\">equal<\/span>\r\n<\/pre><img decoding=\"async\" vspace=\"5\" hspace=\"5\" src=\"https:\/\/blogs.mathworks.com\/images\/loren\/2015\/IncrementalDelaunayBlog2_03.png\" alt=\"\"> <p>So that's basically an incremental change we made to the triangulation. When prototyping applications like this at the command line we may find adding a few points to a large triangulation can take almost as long as creating the full triangulation of all the points.  Why is that, surely the operation should be more efficient?<\/p><h4>When is Incremental Delaunay Important?<a name=\"46f99a36-fd9f-42de-8ea0-badc301144c4\"><\/a><\/h4><p>There are some practical applications that rely on this expected level of efficiency. For example, an incremental algorithm for 2-D Mesh Generation using <a href=\"http:\/\/en.wikipedia.org\/wiki\/Ruppert%27s_algorithm\">Ruppert's Algorithm<\/a>. Delaunay-based algorithms for reconstructing surfaces from point clouds may also be incremental; triangulating an initial set of points and adaptively adding more points to recover the surface. In fact the additive points in these algorithms are often circumcenters and that's why I chose them in the example. But the question remains, shouldn't an incremental addition of a few points to a large triangulation be more efficient than a complete triangulation of all points. Absolutely, this behavior is honored in the scenario where algorithms written in MATLAB are designed to run most efficiently. If we write the algorithm in a function in a file, the incremental change will be performed efficiently. If we prototype at the command line the performance we get may underperform the efficiency we get from the function-in-a-file format.<\/p><h4>Performance Example of Incremental Delaunay Construction<a name=\"c06cefa9-4e4a-4fdb-a61f-e73d41e2f270\"><\/a><\/h4><p>Let's run a little example to test this; the following code creates a 3-D Delaunay triangulation of a half-million points and subsequently adds 40K points in 4 incremental updates. Here's the output when the code runs at the command line:<\/p><p><tt>timeToCreateDelaunay =<\/tt><\/p><pre>   9.7799<\/pre><p><tt>timeToIncrementallyAddPoints =<\/tt><\/p><pre>  10.1267<\/pre><p>When I put the code in a function in a file, execution gives the following:<\/p><pre class=\"language-matlab\">\r\n<span class=\"keyword\">function<\/span> DelaunayIncrementalTest()\r\n    tic; \r\n    dt = delaunayTriangulation(rand(500000,3)); \r\n    timeToCreateDelaunay = toc\r\n\r\n    tic;\r\n    dt.Points(end+(1:10000),:) = rand(10000,3);\r\n    dt.Points(end+(1:10000),:) = rand(10000,3);\r\n    dt.Points(end+(1:10000),:) = rand(10000,3);\r\n    dt.Points(end+(1:10000),:) = rand(10000,3);\r\n    timeToIncrementallyAddPoints = toc\r\n<span class=\"keyword\">end<\/span>\r\n\r\n<\/pre><pre class=\"codeinput\">DelaunayIncrementalTest()\r\n<\/pre><pre class=\"codeoutput\">timeToCreateDelaunay =\r\n       9.9275\r\ntimeToIncrementallyAddPoints =\r\n       1.1396\r\n<\/pre><p>Why the significant performance improvement? MATLAB can execute code more efficiently when it is in the function-in-file format. Loren has a past blog on <a href=\"https:\/\/blogs.mathworks.com\/loren\/2007\/03\/22\/in-place-operations-on-data\/?s_tid=Blog_Loren_Category\">In-place Operations<\/a> that highlights the behavior that improves the efficiency here. So to gauge the performance of your MATLAB code it's good to structure your code into functions that reside in files. Then run the performance analyzer and eliminate the bottlenecks.<\/p><h4>Your Need for Geometric Tools?<a name=\"6a4ed5b2-13cd-4789-abf5-92af268621e8\"><\/a><\/h4><p>Does your work involve triangulations or geometric computing? Does your application area require geometric tools or features that are not well supported in core MATLAB? Close the loop and share your experience <a href=\"https:\/\/blogs.mathworks.com\/loren\/?p=1190#respond\">here<\/a>; hearing what users' need is the compass that charts our feature enhancements!<\/p><script language=\"JavaScript\"> <!-- \r\n    function grabCode_d43b137209db4553b86d5a34b1135f5c() {\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='d43b137209db4553b86d5a34b1135f5c ' + '##### ' + 'SOURCE BEGIN' + ' #####';\r\n        t2='##### ' + 'SOURCE END' + ' #####' + ' d43b137209db4553b86d5a34b1135f5c';\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        copyright = 'Copyright 2015 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 copyright line at the bottom if specified.\r\n        if (copyright.length > 0) {\r\n            d.writeln('');\r\n            d.writeln('%%');\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     --> <\/script><p style=\"text-align: right; font-size: xx-small; font-weight:lighter;   font-style: italic; color: gray\"><br><a href=\"javascript:grabCode_d43b137209db4553b86d5a34b1135f5c()\"><span style=\"font-size: x-small;        font-style: italic;\">Get \r\n      the MATLAB code <noscript>(requires JavaScript)<\/noscript><\/span><\/a><br><br>\r\n      Published with MATLAB&reg; R2015a<br><\/p><\/div><!--\r\nd43b137209db4553b86d5a34b1135f5c ##### SOURCE BEGIN #####\r\n%% Incremental Delaunay Construction\r\n% I'm happy to welcome back Damian Sheehy as guest blogger. Last time\r\n% Damian wrote about how Natural Neighbor interpolation addresses FAQs in\r\n% scattered data interpolation. In this blog he will answer a FAQ on\r\n% adaptively editing a Delaunay triangulation.\r\n% \r\n%% Is |delaunayTriangulation| More Efficient than |delaunay|?\r\n% A technical support question that occasionally crops up asks about the\r\n% best and most efficient way to construct a Delaunay triangulation by\r\n% adding points to an existing triangulation. We call this operation\r\n% incremental editing as the goal is to add or remove points in an\r\n% incremental manner as opposed to recreating from scratch; for example\r\n% calling <https:\/\/www.mathworks.com\/help\/matlab\/ref\/delaunay.html\r\n% |delaunay|> multiple times. The documentation for the\r\n% <https:\/\/www.mathworks.com\/help\/matlab\/ref\/delaunaytriangulation-class.html\r\n% |delaunayTriangulation|> class provides examples that show the syntax\r\n% that allows you to edit a Delaunay triangulation by adding or removing\r\n% points. Let's look at a simple 2-D example to highlight the concept.\r\n% \r\n% In this example we will load |trimesh2d.mat| which ships with MATLAB. The\r\n% file contains points with coordinates (x, y) and triangulation edge\r\n% constraints that define the boundary of a domain. Let's triangulate the\r\n% data and take a look at that.\r\nload trimesh2d\r\ndt = delaunayTriangulation(x,y,Constraints);\r\ntin = find(dt.isInterior());\r\ntriplot(dt(tin,:),dt.Points(:,1),dt.Points(:,2))\r\naxis equal\r\n\r\n%% \r\n% Suppose we want to add the circumcenter points to this triangulation. The\r\n% circumcenter of a triangle is the center of a circumscribed circle that\r\n% passes through the vertices of the triangle. The |delaunayTriangulation|\r\n% class provides a method to compute them -\r\n% <https:\/\/www.mathworks.com\/help\/matlab\/ref\/triangulation.circumcenter.html\r\n% |delaunayTriangulation\/circumcenter|>. Compute them using this method and\r\n% add them to the plot. Note, some triangles may share the same\r\n% circumcenter, so I will call the\r\n% <https:\/\/www.mathworks.com\/help\/matlab\/ref\/uniquetol.html |uniquetol|>\r\n% function to eliminate the near-duplicates.\r\n\r\ncc = dt.circumcenter(tin);\r\ncc = uniquetol(cc, 'ByRows',true);\r\nhold on\r\nplot(cc(:,1),cc(:,2),'.r')\r\nhold off\r\n\r\n%% \r\n% Now, I will add the circumcenter points to the triangulation and plot the\r\n% result in a new figure.\r\nnumcc = size(cc,1);\r\ndt.Points(end+(1:numcc),:) = cc;\r\ntin = find(dt.isInterior());\r\nfigure\r\ntriplot(dt(tin,:),dt.Points(:,1),dt.Points(:,2))\r\naxis equal\r\n\r\n%% \r\n% So that's basically an incremental change we made to the triangulation.\r\n% When prototyping applications like this at the command line we may find\r\n% adding a few points to a large triangulation can take almost as long as\r\n% creating the full triangulation of all the points.  Why is that, surely\r\n% the operation should be more efficient?\r\n%\r\n%% When is Incremental Delaunay Important?\r\n% There are some practical applications that rely on this expected level of\r\n% efficiency. For example, an incremental algorithm for 2-D Mesh Generation\r\n% using <http:\/\/en.wikipedia.org\/wiki\/Ruppert%27s_algorithm Ruppert's\r\n% Algorithm>. Delaunay-based algorithms for reconstructing surfaces from\r\n% point clouds may also be incremental; triangulating an initial set of\r\n% points and adaptively adding more points to recover the surface. In fact\r\n% the additive points in these algorithms are often circumcenters and\r\n% that's why I chose them in the example. But the question remains,\r\n% shouldn't an incremental addition of a few points to a large\r\n% triangulation be more efficient than a complete triangulation of all\r\n% points. Absolutely, this behavior is honored in the scenario where\r\n% algorithms written in MATLAB are designed to run most efficiently. If we\r\n% write the algorithm in a function in a file, the incremental change will\r\n% be performed efficiently. If we prototype at the command line the\r\n% performance we get may underperform the efficiency we get from the\r\n% function-in-a-file format.\r\n%% Performance Example of Incremental Delaunay Construction\r\n% Let's run a little example to test this; the following code creates a 3-D\r\n% Delaunay triangulation of a half-million points and subsequently adds 40K\r\n% points in 4 incremental updates. Here's the output when the code runs at\r\n% the command line:\r\n% \r\n% |timeToCreateDelaunay =|\r\n% \r\n%     9.7799\r\n% \r\n% |timeToIncrementallyAddPoints =|\r\n% \r\n%    10.1267\r\n%\r\n%\r\n% When I put the code in a function in a file, execution gives the\r\n% following:\r\n%\r\n% <include>DelaunayIncrementalTest.m<\/include>\r\n%\r\nDelaunayIncrementalTest()\r\n\r\n%%\r\n% Why the significant performance improvement? MATLAB can execute code more\r\n% efficiently when it is in the function-in-file format. Loren has a past\r\n% blog on\r\n% <https:\/\/blogs.mathworks.com\/loren\/2007\/03\/22\/in-place-operations-on-data\/?s_tid=Blog_Loren_Category\r\n% In-place Operations> that highlights the behavior that improves the\r\n% efficiency here. So to gauge the performance of your MATLAB code it's\r\n% good to structure your code into functions that reside in files. Then run\r\n% the performance analyzer and eliminate the bottlenecks.\r\n\r\n%% Your Need for Geometric Tools?\r\n% Does your work involve triangulations or geometric computing? Does your\r\n% application area require geometric tools or features that are not well\r\n% supported in core MATLAB? Close the loop and share your experience\r\n% <https:\/\/blogs.mathworks.com\/loren\/?p=1190#respond here>; hearing what\r\n% users' need is the compass that charts our feature enhancements!\r\n\r\n\r\n\r\n\r\n\r\n##### SOURCE END ##### d43b137209db4553b86d5a34b1135f5c\r\n-->","protected":false},"excerpt":{"rendered":"<div class=\"overview-image\"><img decoding=\"async\"  class=\"img-responsive\" src=\"https:\/\/blogs.mathworks.com\/images\/loren\/2015\/IncrementalDelaunayBlog2_03.png\" onError=\"this.style.display ='none';\" \/><\/div><!--introduction--><p>I'm happy to welcome back Damian Sheehy as guest blogger. Last time Damian wrote about how Natural Neighbor interpolation addresses FAQs in scattered data interpolation. In this blog he will answer a FAQ on adaptively editing a Delaunay triangulation.... <a class=\"read-more\" href=\"https:\/\/blogs.mathworks.com\/loren\/2015\/07\/15\/incremental-delaunay-construction\/\">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,23],"tags":[],"_links":{"self":[{"href":"https:\/\/blogs.mathworks.com\/loren\/wp-json\/wp\/v2\/posts\/1190"}],"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=1190"}],"version-history":[{"count":4,"href":"https:\/\/blogs.mathworks.com\/loren\/wp-json\/wp\/v2\/posts\/1190\/revisions"}],"predecessor-version":[{"id":2448,"href":"https:\/\/blogs.mathworks.com\/loren\/wp-json\/wp\/v2\/posts\/1190\/revisions\/2448"}],"wp:attachment":[{"href":"https:\/\/blogs.mathworks.com\/loren\/wp-json\/wp\/v2\/media?parent=1190"}],"wp:term":[{"taxonomy":"category","embeddable":true,"href":"https:\/\/blogs.mathworks.com\/loren\/wp-json\/wp\/v2\/categories?post=1190"},{"taxonomy":"post_tag","embeddable":true,"href":"https:\/\/blogs.mathworks.com\/loren\/wp-json\/wp\/v2\/tags?post=1190"}],"curies":[{"name":"wp","href":"https:\/\/api.w.org\/{rel}","templated":true}]}}