{"id":10926,"date":"2019-08-23T09:00:45","date_gmt":"2019-08-23T13:00:45","guid":{"rendered":"https:\/\/blogs.mathworks.com\/pick\/?p=10926"},"modified":"2019-08-20T10:50:45","modified_gmt":"2019-08-20T14:50:45","slug":"polygon2voxel","status":"publish","type":"post","link":"https:\/\/blogs.mathworks.com\/pick\/2019\/08\/23\/polygon2voxel\/","title":{"rendered":"Polygon2Voxel"},"content":{"rendered":"<div class=\"content\">\n<p><a href=\"http:\/\/www.mathworks.com\/matlabcentral\/profile\/authors\/3208495\">Sean<\/a>&#8216;s pick this week is <a href=\"http:\/\/www.mathworks.com\/matlabcentral\/fileexchange\/24086\"><tt>polygon2voxel<\/tt><\/a> by <a href=\"http:\/\/www.mathworks.com\/matlabcentral\/profile\/authors\/1097878\">Dirk-Jan Kroon<\/a>.<\/p>\n<p>I was recently asked: &#8220;How can you find the intersecting volume of two polyhedra?&#8221;<\/p>\n<p>R2017b introduced <a href=\"https:\/\/www.mathworks.com\/help\/releases\/R2019a\/matlab\/ref\/polyshape.html\"><tt>polyshape<\/tt><\/a> to MATLAB to make working with two dimensional polygons easy. However, it does not support 3D.<\/p>\n<p>One way to get an approximation to the volume would be to convert the polyhedra to a three dimensional image and find the overlapping voxels. Dirk-Jan&#8217;s <tt>polygon2voxel<\/tt> helps with the first step of this.<\/p>\n<p>First, let&#8217;s create two 3D shapes, one convex, and one not.<\/p>\n<pre style=\"background: #F9F7F3; padding: 10px; border: 1px solid #c8c8c8;\">xyz = gallery(<span style=\"color: #a020f0;\">'uniformdata'<\/span>, [30 3], 0)*100;\r\nDT = delaunayTriangulation(xyz);\r\nfigure\r\ntetramesh(DT)<\/pre>\n<p><img decoding=\"async\" src=\"http:\/\/blogs.mathworks.com\/images\/pick\/Sean\/mainPolygon2Voxel\/mainPolygon2Voxel_01.png\" hspace=\"5\" vspace=\"5\" \/><\/p>\n<pre style=\"background: #F9F7F3; padding: 10px; border: 1px solid #c8c8c8;\">xyz = gallery(<span style=\"color: #a020f0;\">'uniformdata'<\/span>, [30 3], 1)*100;\r\nAS = alphaShape(xyz, 40);\r\nfigure\r\nplot(AS)<\/pre>\n<p><img decoding=\"async\" src=\"http:\/\/blogs.mathworks.com\/images\/pick\/Sean\/mainPolygon2Voxel\/mainPolygon2Voxel_02.png\" hspace=\"5\" vspace=\"5\" \/><\/p>\n<p>Now let&#8217;s convert these to voxels. I will use just the freeBoundary, i.e. the surface, of the delaunay triamgulation.<\/p>\n<pre style=\"background: #F9F7F3; padding: 10px; border: 1px solid #c8c8c8;\">[f, t] = freeBoundary(DT);\r\nConvexImage = polygon2voxel(struct(<span style=\"color: #a020f0;\">'faces'<\/span>, f, <span style=\"color: #a020f0;\">'vertices'<\/span>, t), [100 100 100], <span style=\"color: #a020f0;\">'auto'<\/span>);<\/pre>\n<p>For an alphaShape, <tt>boundaryFacets<\/tt> will return the surface.<\/p>\n<pre style=\"background: #F9F7F3; padding: 10px; border: 1px solid #c8c8c8;\">[f, t] = boundaryFacets(AS);\r\nConcaveImage = polygon2voxel(struct(<span style=\"color: #a020f0;\">'faces'<\/span>, f, <span style=\"color: #a020f0;\">'vertices'<\/span>, t), [100 100 100], <span style=\"color: #a020f0;\">'auto'<\/span>);<\/pre>\n<p>We can watch both volumes as a movie.<\/p>\n<pre style=\"background: #F9F7F3; padding: 10px; border: 1px solid #c8c8c8;\">implay([ConvexImage ConcaveImage])<\/pre>\n<p><embed src=\"https:\/\/blogs.mathworks.com\/pick\/files\/Unfilled.mp4\" height=\"400\" width=\"400\"><\/p>\n<p>The holes in the volumes have not been filled; <tt>imfill<\/tt> can do that for us.<\/p>\n<pre style=\"background: #F9F7F3; padding: 10px; border: 1px solid #c8c8c8;\">ConvexImage = imfill(ConvexImage, <span style=\"color: #a020f0;\">'holes'<\/span>);\r\nConcaveImage = imfill(ConcaveImage, <span style=\"color: #a020f0;\">'holes'<\/span>);\r\n\r\nimplay([ConvexImage ConcaveImage])<\/pre>\n<p><embed src=\"https:\/\/blogs.mathworks.com\/pick\/files\/Filled.mp4\" height=\"400\" width=\"400\"><\/p>\n<p>Now logical indexing can be used to find the overlapping volume:<\/p>\n<pre style=\"background: #F9F7F3; padding: 10px; border: 1px solid #c8c8c8;\">VoxelIntersection = ConvexImage &amp; ConcaveImage;\r\nfigure\r\nisosurface(VoxelIntersection)<\/pre>\n<p><img decoding=\"async\" src=\"http:\/\/blogs.mathworks.com\/images\/pick\/Sean\/mainPolygon2Voxel\/mainPolygon2Voxel_05.png\" hspace=\"5\" vspace=\"5\" \/><\/p>\n<p>The volume could be computed as just the sum of the voxels:<\/p>\n<pre style=\"background: #F9F7F3; padding: 10px; border: 1px solid #c8c8c8;\">volumevoxels = sum(VoxelIntersection, <span style=\"color: #a020f0;\">'all'<\/span>);\r\ndisp(<span style=\"color: #a020f0;\">\"The intersecting volume is: \"<\/span> + volumevoxels + <span style=\"color: #a020f0;\">\"units\"<\/span>)<\/pre>\n<pre style=\"font-style: oblique;\">The intersecting volume is: 272741units\r\n<\/pre>\n<h3>Comments<a name=\"9\"><\/a><\/h3>\n<p>Do you have a use case for 3D polyshapes in MATLAB? Give it a try and let us know what you think <a href=\"http:\/\/blogs.mathworks.com\/pick\/?p=10926#respond\">here<\/a> or leave a <a href=\"http:\/\/www.mathworks.com\/matlabcentral\/fileexchange\/24086#comments\">comment<\/a> for Dirk-Jan.<\/p>\n<p><script language=\"JavaScript\">\n<!--\n\n    function grabCode_a02a184672f542b49b72f4bc08693c91() {\n        \/\/ Remember the title so we can use it in the new page\n        title = document.title;\n\n        \/\/ Break up these strings so that their presence\n        \/\/ in the Javascript doesn't mess up the search for\n        \/\/ the MATLAB code.\n        t1='a02a184672f542b49b72f4bc08693c91 ' + '##### ' + 'SOURCE BEGIN' + ' #####';\n        t2='##### ' + 'SOURCE END' + ' #####' + ' a02a184672f542b49b72f4bc08693c91';\n    \n        b=document.getElementsByTagName('body')[0];\n        i1=b.innerHTML.indexOf(t1)+t1.length;\n        i2=b.innerHTML.indexOf(t2);\n \n        code_string = b.innerHTML.substring(i1, i2);\n        code_string = code_string.replace(\/REPLACE_WITH_DASH_DASH\/g,'--');\n\n        \/\/ Use \/x3C\/g instead of the less-than character to avoid errors \n        \/\/ in the XML parser.\n        \/\/ Use '\\x26#60;' instead of '<' so that the XML parser\n        \/\/ doesn't go ahead and substitute the less-than character. \n        code_string = code_string.replace(\/\\x3C\/g, '\\x26#60;');\n\n        author = 'Sean de Wolski';\n        copyright = 'Copyright 2019 The MathWorks, Inc.';\n\n        w = window.open();\n        d = w.document;\n        d.write('<\/p>\n\n\n\n\n<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\n\n\n\n<p>\\n');\n      \n      d.title = title + ' (MATLAB code)';\n      d.close();\n      }   \n      \n-->\n<\/script><\/p>\n<p style=\"text-align: right; font-size: xx-small; font-weight: lighter; font-style: italic; color: gray;\">\n<p>Published with MATLAB\u00ae R2019b<\/p>\n<\/div>\n<p><!--\na02a184672f542b49b72f4bc08693c91 ##### SOURCE BEGIN #####\n%% Polygon2Voxel\n%\n% <http:\/\/www.mathworks.com\/matlabcentral\/profile\/authors\/3208495 Sean>'s pick this week is\n% <http:\/\/www.mathworks.com\/matlabcentral\/fileexchange\/24086 |polygon2voxel|> by\n% <http:\/\/www.mathworks.com\/matlabcentral\/profile\/authors\/1097878 Dirk-Jan Kroon>.\n%\n\n%%\n%\n% I was recently asked: \"How can you find the intersecting volume of two\n% polyhedra?\"\n%\n% R2017b introduced <https:\/\/www.mathworks.com\/help\/releases\/R2019a\/matlab\/ref\/polyshape.html |polyshape|> to MATLAB to make working with two\n% dimensional polygons easy.  However, it does not support 3D.\n%\n% One way to get an approximation to the volume would be to convert the\n% polyhedra to a three dimensional image and find the overlapping voxels.\n% Dirk-Jan's |polygon2voxel| helps with the first step of this.\n%\n% First, let's create two 3D shapes, one convex, and one not.\n\nxyz = gallery('uniformdata', [30 3], 0)*100;\nDT = delaunayTriangulation(xyz);\nfigure\ntetramesh(DT)\n%%\nxyz = gallery('uniformdata', [30 3], 1)*100;\nAS = alphaShape(xyz, 40);\nfigure\nplot(AS)\n\n%%\n% Now let's convert these to voxels.  I will use just the freeBoundary,\n% i.e. the surface, of the delaunay triamgulation.\n[f, t] = freeBoundary(DT);\nConvexImage = polygon2voxel(struct('faces', f, 'vertices', t), [100 100 100], 'auto');\n\n%%\n% For an alphaShape, |boundaryFacets| will return the surface.\n[f, t] = boundaryFacets(AS);\nConcaveImage = polygon2voxel(struct('faces', f, 'vertices', t), [100 100 100], 'auto');\n\n%%\n% We can watch both volumes as a movie.\nimplay([ConvexImage ConcaveImage])\n\n%%\n% The holes in the volumes have not been filled; |imfill| can do that for\n% us.\nConvexImage = imfill(ConvexImage, 'holes');\nConcaveImage = imfill(ConcaveImage, 'holes');\n\nimplay([ConvexImage ConcaveImage])\n\n%%\n% Now logical indexing can be used to find the overlapping volume:\nVoxelIntersection = ConvexImage & ConcaveImage;\nfigure\nisosurface(VoxelIntersection)\n\n%%\n% The volume could be computed as just the sum of the voxels:\n\nvolumevoxels = sum(VoxelIntersection, 'all');\ndisp(\"The intersecting volume is: \" + volumevoxels + \"units\")\n\n%% Comments\n%\n% Give it a try and let us know what you think\n% <http:\/\/blogs.mathworks.com\/pick\/?p=10926#respond here> or leave a\n% <http:\/\/www.mathworks.com\/matlabcentral\/fileexchange\/24086#comments % comment> for Dirk-Jan.\n%\n\n##### SOURCE END ##### a02a184672f542b49b72f4bc08693c91\n--><\/p>\n","protected":false},"excerpt":{"rendered":"<div class=\"overview-image\"><img decoding=\"async\"  class=\"img-responsive\" src=\"http:\/\/blogs.mathworks.com\/images\/pick\/Sean\/mainPolygon2Voxel\/mainPolygon2Voxel_01.png\" onError=\"this.style.display ='none';\" \/><\/div>\n<p>\nSean&#8216;s pick this week is polygon2voxel by Dirk-Jan Kroon.<br \/>\nI was recently asked: &#8220;How can you find the intersecting volume of two polyhedra?&#8221;<br \/>\nR2017b introduced polyshape to MATLAB&#8230; <a class=\"read-more\" href=\"https:\/\/blogs.mathworks.com\/pick\/2019\/08\/23\/polygon2voxel\/\">read more >><\/a><\/p>\n","protected":false},"author":87,"featured_media":0,"comment_status":"open","ping_status":"closed","sticky":false,"template":"","format":"standard","meta":[],"categories":[16,6],"tags":[],"_links":{"self":[{"href":"https:\/\/blogs.mathworks.com\/pick\/wp-json\/wp\/v2\/posts\/10926"}],"collection":[{"href":"https:\/\/blogs.mathworks.com\/pick\/wp-json\/wp\/v2\/posts"}],"about":[{"href":"https:\/\/blogs.mathworks.com\/pick\/wp-json\/wp\/v2\/types\/post"}],"author":[{"embeddable":true,"href":"https:\/\/blogs.mathworks.com\/pick\/wp-json\/wp\/v2\/users\/87"}],"replies":[{"embeddable":true,"href":"https:\/\/blogs.mathworks.com\/pick\/wp-json\/wp\/v2\/comments?post=10926"}],"version-history":[{"count":6,"href":"https:\/\/blogs.mathworks.com\/pick\/wp-json\/wp\/v2\/posts\/10926\/revisions"}],"predecessor-version":[{"id":10942,"href":"https:\/\/blogs.mathworks.com\/pick\/wp-json\/wp\/v2\/posts\/10926\/revisions\/10942"}],"wp:attachment":[{"href":"https:\/\/blogs.mathworks.com\/pick\/wp-json\/wp\/v2\/media?parent=10926"}],"wp:term":[{"taxonomy":"category","embeddable":true,"href":"https:\/\/blogs.mathworks.com\/pick\/wp-json\/wp\/v2\/categories?post=10926"},{"taxonomy":"post_tag","embeddable":true,"href":"https:\/\/blogs.mathworks.com\/pick\/wp-json\/wp\/v2\/tags?post=10926"}],"curies":[{"name":"wp","href":"https:\/\/api.w.org\/{rel}","templated":true}]}}