{"id":2275,"date":"2013-06-20T16:14:21","date_gmt":"2013-06-20T21:14:21","guid":{"rendered":"https:\/\/blogs.mathworks.com\/community\/?p=2275"},"modified":"2021-10-28T11:38:26","modified_gmt":"2021-10-28T15:38:26","slug":"paul-prints-the-l-shaped-membrane","status":"publish","type":"post","link":"https:\/\/blogs.mathworks.com\/community\/2013\/06\/20\/paul-prints-the-l-shaped-membrane\/","title":{"rendered":"Paul Prints the L-Shaped Membrane"},"content":{"rendered":"<p><em>As promised, Paul Kassebaum is back this week with an in-depth discussion of how to get from a mathematical object in MATLAB to a solid object you can hold in your hand. Paul is a maker in the truest sense of the word. If there are a hundred weird and unexpected obstacles between him and the thing he wants to create, he works through them methodically and documents his work carefully. I get incredibly frustrated trying to tune up my lawn mower. Paul, on the other hand, could probably build a working lawn mower from a pair of scissors, a box of binder clips, and a broken copy machine.<\/em><\/p>\n<h2>3D Printing from MATLAB<\/h2>\n<div class=\"content\">\n<p><i>by Paul Kassebaum<\/i><\/p>\n<p>Many 3D-printed objects are designed by combining and sculpting primitive shapes such as spheres and cones using CAD software. This post will show how to use MATLAB to create 3D-printable objects based on equations or data, a task that most CAD programs are not designed to do.<\/p>\n<p>This post will be structured as follows. First, we will explain the file format read by 3D printers, called stereo-lithography files, or STL-files. We will then show how to generate these files starting from matrix-based surface plots created from such sources as photographs and solutions to partial differential equations such as the beloved L-shaped membrane. Finally, we will be able to hold a MATLAB plot in our hands.<\/p>\n<h3>Contents<\/h3>\n<div>\n<ul>\n<li><a href=\"#1f9d71f8-8030-43d7-a639-a0388376118e\">Stereo-lithography (STL) Files<\/a><\/li>\n<li><a href=\"#0d2c8d8b-352f-4b12-a0ba-5b5922a8795e\">Printing Surface Plots<\/a><\/li>\n<li><a href=\"#99ba97f6-522b-4298-98a4-6b8233c16a14\">MATLAB Materialized<\/a><\/li>\n<\/ul>\n<\/div>\n<h4>Stereo-lithography (STL) Files<a name=\"1f9d71f8-8030-43d7-a639-a0388376118e\"><\/a><\/h4>\n<p>STL-files describe a closed surface in terms of triangular faces. These files consist of a list of triangles, with each triangle described by the cartesian coordinates of its three vertices and a normal vector oriented outward from the closed surface. Here is an excerpt from an ASCII STL-file:<\/p>\n<pre>facet normal nx ny nz\r\nouter loop\r\nvertex v1x v1y v1z\r\nvertex v2x v2y v2z\r\nvertex v3x v3y v3z\r\nendloop\r\nendfacet<\/pre>\n<p>In practice, it might look like this:<\/p>\n<pre>facet normal 6.6998060E-01 -6.6246430E-01 3.3506277E-01\r\nouter loop\r\nvertex 7.9166667E-01 2.5000000E-01 9.4328269E-01\r\nvertex 7.5000000E-01 2.5000000E-01 1.0265980E+00\r\nvertex 7.9166667E-01 2.0833333E-01 8.6090207E-01\r\nendloop\r\nendfacet<\/pre>\n<p>Before being useful for a 3D printer, a surface described by an STL-file must be sliced into layers defining the path traced out by the printer. Many 3D printers come with their own slicing program to perform this translation, so we will not address this stage. It is worth noting that depending upon the type of printer being used, the orientation of the surface may factor into the quality of the print.<\/p>\n<h4>Printing Surface Plots<a name=\"0d2c8d8b-352f-4b12-a0ba-5b5922a8795e\"><\/a><\/h4>\n<p>The most natural way to create surface plots in MATLAB is to use a matrix of height values. Consider the surface plot of the MathWorks logo, based on the L-shaped membrane.<\/p>\n<pre class=\"codeinput\">n = 12; <span class=\"comment\">% number of partitions in each dimension.<\/span>\r\n[X,Y] = meshgrid(linspace(0,1,2*n+1));\r\nL = (40\/51\/0.9)*membrane(1,n);\r\nsurf(X,Y,L);\r\ncolormap <span class=\"string\">pink<\/span>;\r\nset(gca,<span class=\"string\">'dataAspectRatio'<\/span>, [1 1 1]);\r\n<\/pre>\n<p><img decoding=\"async\" src=\"https:\/\/blogs.mathworks.com\/community\/files\/matlab3DPrinting_01.png\" alt=\"\" hspace=\"5\" vspace=\"5\" \/><\/p>\n<p>The first step in translating a matrix-based plot into an STL-file is to break up each square element in the mesh into two triangular elements. We can make use of the function DELAUNAY to create a Delaunay triangulation of the rectilinear mesh.<\/p>\n<pre class=\"codeinput\">faces   = delaunay(X,Y);\r\npatches = trisurf(faces,X,Y,L);\r\nset(gca,<span class=\"string\">'dataAspectRatio'<\/span>, [1 1 1]);\r\n<\/pre>\n<p><img decoding=\"async\" src=\"https:\/\/blogs.mathworks.com\/community\/files\/matlab3DPrinting_02.png\" alt=\"\" hspace=\"5\" vspace=\"5\" \/><\/p>\n<p>Notice that this mesh has no thickness, so it cannot yet be used to make a solid 3D print. We can make a shell out of this surface as follows. First, project all the vertices of the triangles downward along normal vectors to create a second surface with no thickness beneath the original. Then connect the two surfaces along their boundaries to define a third surface.<\/p>\n<p>Let us begin by creating a 3-by-3-by-n tensor called &#8216;facets&#8217; that will allow us to easily index through the 3 cartesian coordinates of the 3 vertices of the n triangles that make up our surface.<\/p>\n<pre class=\"codeinput\">vertices = get(patches,<span class=\"string\">'vertices'<\/span>);\r\nfacets = vertices';\r\nfacets = reshape(facets(:,faces'), 3, 3, []);\r\n<\/pre>\n<p>We will say that the normal at a vertex is the average of the normals to the facets that share this vertex. So, we will need to know the normals for each facet. We can calculate the normals of each facet by taking the cross-product of two edges of the facet, taking note of their orientation to ensure that the normals of the whole surface are consistent.<\/p>\n<pre class=\"codeinput\"><span class=\"comment\">% SQUEEZE compacts empty dimensions.<\/span>\r\nedge1 = squeeze(facets(:,2,:) - facets(:,1,:));\r\nedge2 = squeeze(facets(:,3,:) - facets(:,1,:));\r\nnormals = edge1([2 3 1],:) .* edge2([3 1 2],:)<span class=\"keyword\">...<\/span>\r\n        - edge2([2 3 1],:) .* edge1([3 1 2],:);\r\nnormals = bsxfun(@times,<span class=\"keyword\">...<\/span>\r\n  normals, 1 .\/ sqrt(sum(normals .* normals, 1)));\r\n<\/pre>\n<p>We can compute the normal at each vertex by averaging the normals of its neighboring facets.<\/p>\n<pre class=\"codeinput\">meanNormal = zeros(3,length(vertices)); <span class=\"comment\">% pre-allocate memory.<\/span>\r\n<span class=\"keyword\">for<\/span> k = 1:length(vertices)\r\n  <span class=\"comment\">% Find all faces shared by a vertex<\/span>\r\n  [sharedFaces,~] = find(faces == k);\r\n  <span class=\"comment\">% Compute the mean normal of all faces shared by a vertex<\/span>\r\n  meanNormal(:,k) = mean(normals(:,sharedFaces),2);\r\n<span class=\"keyword\">end<\/span>\r\nmeanNormal = bsxfun(@times, meanNormal,<span class=\"keyword\">...<\/span>\r\n  1 .\/ sqrt(sum(meanNormal .* meanNormal, 1)));\r\n<\/pre>\n<p>Now we merely need to copy and shift all of the vertices downward along their normals. we will call these new vertices &#8216;underVertices&#8217;.<\/p>\n<pre class=\"codeinput\">shellThickness = 0.05;\r\nunderVertices  = vertices - shellThickness*meanNormal';\r\n<\/pre>\n<p>Let us see what we have got so far.<\/p>\n<pre class=\"codeinput\">underFaces = delaunay(underVertices(:,1),underVertices(:,2));\r\ntrisurf(underFaces,<span class=\"keyword\">...<\/span>\r\n   underVertices(:,1),<span class=\"keyword\">...<\/span>\r\n   underVertices(:,2),<span class=\"keyword\">...<\/span>\r\n   underVertices(:,3));\r\nset(gca,<span class=\"string\">'dataAspectRatio'<\/span>, [1 1 1],<span class=\"keyword\">...<\/span>\r\n  <span class=\"string\">'xLim'<\/span>,[0 1],<span class=\"string\">'yLim'<\/span>,[0 1]);\r\n<\/pre>\n<p><img decoding=\"async\" src=\"https:\/\/blogs.mathworks.com\/community\/files\/matlab3DPrinting_03.png\" alt=\"\" hspace=\"5\" vspace=\"5\" \/><\/p>\n<p>Notice that there are some extra facets that we did not intend. These arise from the 2D Delaunay triangulation, which connects vertices contained within the convex hull of the whole set of vertices by default. This did not pose a problem in our original surface because the boundaries consisted of a square, which is its own convex hull. The boundary of our new surface is the original square boundary deformed by the translation of each vertex through the normal vectors we calculated. Thus, our new surface has a boundary that is not its own convex hull. MATLAB&#8217;s DELAUNAYTRI class allows us to specify the boundary in 2D, which solves our problem.<\/p>\n<p>To resolve this problem, we will first find the indices of the boundaries of the original surface, which will also index the boundaries of the lower surface. These boundary indices will serve to constrain the Delaunay triangulation and will be used to connect the two surfaces with a new surface. First, let us find the boundary indices.<\/p>\n<pre class=\"codeinput\">boundaryIndices = <span class=\"keyword\">...<\/span>\r\n [find(vertices(:,2) == min(vertices(:,2))); <span class=\"comment\">% min y<\/span>\r\n  find(vertices(:,1) == max(vertices(:,1))); <span class=\"comment\">% max x<\/span>\r\n  find(vertices(:,2) == max(vertices(:,2))); <span class=\"comment\">% max y<\/span>\r\n  find(vertices(:,1) == min(vertices(:,1)))];<span class=\"comment\">% min x<\/span>\r\n<\/pre>\n<p>Next we will rearrange the indices so they parameterize the boundary. That is, we traverse the boundary in a counterclockwise sense as the index increases.<\/p>\n<pre class=\"codeinput\">boundaryIndices = [<span class=\"keyword\">...<\/span>\r\n  boundaryIndices(1:floor(end\/4-1)); <span class=\"comment\">% semi open interval [1, end\/4).<\/span>\r\n  boundaryIndices(floor(end\/4+1):floor(end\/2));<span class=\"comment\">%[end\/4, end\/2)<\/span>\r\n  boundaryIndices(floor(end*3\/4-1):-1:floor(end\/2+1));<span class=\"comment\">%[end\/2,end*3\/4)<\/span>\r\n  boundaryIndices(end-1:-1:floor(end*3\/4+1))]; <span class=\"comment\">%[end*3\/4, end)<\/span>\r\n<\/pre>\n<p>The DELAUNAYTRI constructor expects each constrained edge to be defined in terms of its terminal vertices. We can create a sequence of edges by staggering the boundary vertices.<\/p>\n<pre class=\"codeinput\">constrainedEdges = [boundaryIndices(1:end-1), boundaryIndices(2:end)];\r\nunderFaces = DelaunayTri(<span class=\"keyword\">...<\/span>\r\n  [underVertices(:,1),underVertices(:,2)],constrainedEdges);\r\n<\/pre>\n<p>The DELAUNAYTRI constructor has created a triangulation that consists of two parts: one region within our constrained edges, and another outside of them. Since we only care about the inside region, we will pick it out.<\/p>\n<pre class=\"codeinput\">inside = underFaces.inOutStatus; <span class=\"comment\">% 1 = in, 0 = out.<\/span>\r\nunderFaces = underFaces.Triangulation(inside,:);\r\n<\/pre>\n<p>The normals of the lower surface have the same orientation as the original surface because of the way we constructed it. However, we will be making a closed surface by connecting these two, at which point the normals of the lower surface will need to be reversed to point outward. So, let us flip the lower surface&#8217;s normals.<\/p>\n<pre class=\"codeinput\">underFaces = fliplr(underFaces);\r\n<\/pre>\n<p>Let us see what we have got so far.<\/p>\n<pre class=\"codeinput\">trisurf(underFaces,<span class=\"keyword\">...<\/span>\r\n   underVertices(:,1),<span class=\"keyword\">...<\/span>\r\n   underVertices(:,2),<span class=\"keyword\">...<\/span>\r\n   underVertices(:,3));\r\nset(gca,<span class=\"string\">'dataAspectRatio'<\/span>, [1 1 1],<span class=\"keyword\">...<\/span>\r\n  <span class=\"string\">'xLim'<\/span>,[0 1],<span class=\"string\">'yLim'<\/span>,[0 1]);\r\n<\/pre>\n<p><img decoding=\"async\" src=\"https:\/\/blogs.mathworks.com\/community\/files\/matlab3DPrinting_04.png\" alt=\"\" hspace=\"5\" vspace=\"5\" \/><\/p>\n<p>This looks much better. Now we can move on to connect these two surfaces with a third surface that we will call &#8216;wall&#8217;. The wall will be made up of the boundary vertices of the top and bottom surfaces. We will triangulate these vertices by defining each facet of the wall to have one vertex on one surface and two vertices on the other surface.<\/p>\n<pre class=\"codeinput\">wallVertices = [vertices(boundaryIndices,:);\r\n           underVertices(boundaryIndices,:)];\r\n<span class=\"comment\">% Number of wall vertices on each surface (nwv).<\/span>\r\nnwv          = length(wallVertices)\/2;\r\n<span class=\"comment\">% Allocate memory for wallFaces.<\/span>\r\nwallFaces    = zeros(2*(nwv-1),3);\r\n<span class=\"comment\">% Define the faces.<\/span>\r\n<span class=\"keyword\">for<\/span> k = 1:nwv-1\r\n  wallFaces(k      ,:) = [k+1  ,k      ,k+nwv];\r\n  wallFaces(k+nwv-1,:) = [k+nwv,k+1+nwv,k+1];\r\n<span class=\"keyword\">end<\/span>\r\n<\/pre>\n<p>Let us see what we have got so far.<\/p>\n<pre class=\"codeinput\">trisurf(wallFaces,<span class=\"keyword\">...<\/span>\r\n   wallVertices(:,1),<span class=\"keyword\">...<\/span>\r\n   wallVertices(:,2),<span class=\"keyword\">...<\/span>\r\n   wallVertices(:,3));\r\nset(gca,<span class=\"string\">'dataAspectRatio'<\/span>, [1 1 1],<span class=\"keyword\">...<\/span>\r\n  <span class=\"string\">'xLim'<\/span>,[0 1],<span class=\"string\">'yLim'<\/span>,[0 1]);\r\n<\/pre>\n<p><img decoding=\"async\" src=\"https:\/\/blogs.mathworks.com\/community\/files\/matlab3DPrinting_05.png\" alt=\"\" hspace=\"5\" vspace=\"5\" \/><\/p>\n<p>Now let us assemble our three surfaces into one closed surface that we will call &#8216;shell&#8217;.<\/p>\n<pre class=\"codeinput\"><span class=\"comment\">% Shift indices to concatenate with the original surface.<\/span>\r\nunderFaces = underFaces +   length(vertices);\r\nwallFaces  = wallFaces  + 2*length(vertices);\r\n<span class=\"comment\">% Concatenate the results.<\/span>\r\nshellVertices = [vertices; underVertices; wallVertices];\r\nshellFaces    = [faces;    underFaces;    wallFaces];\r\n<\/pre>\n<p>Most 3D printers require that all of the vertex coordinates be non-negative, so we will shift our shell up to satisfy this convention.<\/p>\n<pre class=\"codeinput\">minZ = min(shellVertices(:,3));\r\nshellVertices = shellVertices<span class=\"keyword\">...<\/span>\r\n  - repmat([0 0 minZ],length(shellVertices),1);\r\n<\/pre>\n<p>Let us look at the final result.<\/p>\n<pre class=\"codeinput\">trisurf(shellFaces,<span class=\"keyword\">...<\/span>\r\n   shellVertices(:,1),<span class=\"keyword\">...<\/span>\r\n   shellVertices(:,2),<span class=\"keyword\">...<\/span>\r\n   shellVertices(:,3));\r\nset(gca,<span class=\"string\">'dataAspectRatio'<\/span>, [1 1 1],<span class=\"keyword\">...<\/span>\r\n  <span class=\"string\">'xLim'<\/span>,[0 1],<span class=\"string\">'yLim'<\/span>,[0 1]);\r\n<\/pre>\n<p><img decoding=\"async\" src=\"https:\/\/blogs.mathworks.com\/community\/files\/matlab3DPrinting_06.png\" alt=\"\" hspace=\"5\" vspace=\"5\" \/><\/p>\n<p>The final stage is to convert this surface into an STL-file. We refer to the excellent entry in the File Exchange called <a href=\"https:\/\/www.mathworks.com\/matlabcentral\/fileexchange\/20922-stlwrite-write-ascii-or-binary-stl-files\">stlwrite<\/a> for converting our surface to an STL-file. Stlwrite is suitable for closed surfaces defined by data, but not for surfaces with boundaries like the one we started with here, because stlwrite will generate a mesh with no thickness.<\/p>\n<pre class=\"codeinput\"><span class=\"comment\">% Name your STL-file<\/span>\r\nfilename = <span class=\"string\">'MathWorksLogo.stl'<\/span>;\r\n\r\n<span class=\"comment\">% Create the facets.<\/span>\r\nshellFacets = shellVertices';\r\nshellFacets = reshape(shellFacets(:,shellFaces'), 3, 3, []);\r\n\r\n<span class=\"comment\">% Compute their normals.<\/span>\r\nedge1 = squeeze(shellFacets(:,2,:) - shellFacets(:,1,:));\r\nedge2 = squeeze(shellFacets(:,3,:) - shellFacets(:,1,:));\r\nshellNormals = edge1([2 3 1],:) .* edge2([3 1 2],:)<span class=\"keyword\">...<\/span>\r\n        - edge2([2 3 1],:) .* edge1([3 1 2],:);\r\nshellNormals = bsxfun(@times,<span class=\"keyword\">...<\/span>\r\n  shellNormals, 1 .\/ sqrt(sum(shellNormals .* shellNormals, 1)));\r\n\r\n<span class=\"comment\">% Associate each facet with its normal.<\/span>\r\nshellFacets = cat(2, reshape(shellNormals, 3, 1, []), shellFacets);\r\n\r\n<span class=\"comment\">% Open the file for writing<\/span>\r\nfid = fopen(filename,<span class=\"string\">'wb+'<\/span>);\r\n\r\n<span class=\"comment\">% Write the file contents.<\/span>\r\n<span class=\"comment\">% Write HEADER.<\/span>\r\nfprintf(fid,<span class=\"string\">'solid %s\\r\\n'<\/span>,filename);\r\n<span class=\"comment\">% Write DATA.<\/span>\r\nfprintf(fid,[<span class=\"keyword\">...<\/span>\r\n<span class=\"string\">'facet normal %.7E %.7E %.7E\\r\\n'<\/span> <span class=\"keyword\">...<\/span>\r\n<span class=\"string\">'outer loop\\r\\n'<\/span> <span class=\"keyword\">...<\/span>\r\n<span class=\"string\">'vertex %.7E %.7E %.7E\\r\\n'<\/span> <span class=\"keyword\">...<\/span>\r\n<span class=\"string\">'vertex %.7E %.7E %.7E\\r\\n'<\/span> <span class=\"keyword\">...<\/span>\r\n<span class=\"string\">'vertex %.7E %.7E %.7E\\r\\n'<\/span> <span class=\"keyword\">...<\/span>\r\n<span class=\"string\">'endloop\\r\\n'<\/span> <span class=\"keyword\">...<\/span>\r\n<span class=\"string\">'endfacet\\r\\n'<\/span>], shellFacets);\r\n<span class=\"comment\">% Write FOOTER.<\/span>\r\nfprintf(fid,<span class=\"string\">'endsolid %s\\r\\n'<\/span>,filename);\r\n\r\n<span class=\"comment\">% Close the file.<\/span>\r\nfclose(fid);\r\n<\/pre>\n<h4>MATLAB Materialized<a name=\"99ba97f6-522b-4298-98a4-6b8233c16a14\"><\/a><\/h4>\n<p>With your STL-file in hand, you can now make your print! If you do not have your own 3D printer, you can look for your nearest maker-space to see if they have one, or send your STL-file off to a 3D printing service provider such as <a href=\"http:\/\/www.shapeways.com\">Shapeways<\/a>, <a href=\"http:\/\/www.ponoko.com\">Ponoko<\/a>, <a href=\"http:\/\/www.sculpteo.com\">Sculpteo<\/a>, <a title=\"http:\/\/www.i.materialise.com (broken)\">I.Materialize<\/a>, ZoomRP, or RedEye. I am a member of the maker-space <a href=\"http:\/\/artisansasylum.com\/\">Artisan&#8217;s Asylum<\/a>, which has a nice 3D printer, among many other tools. Here is a photo of my completed L-shaped membrane.<\/p>\n<p><img decoding=\"async\" src=\"https:\/\/blogs.mathworks.com\/community\/files\/shape2.png\" alt=\"\" hspace=\"5\" vspace=\"5\" \/><\/p>\n<p>Good luck and have fun turning your MATLAB plots into tangible objects!<\/p>\n<p style=\"text-align: right; font-size: xx-small; font-weight: lighter; font-style: italic; color: gray;\">\n<p class=\"footer\">Published with MATLAB\u00ae R2013a<\/p>\n<\/div>\n","protected":false},"excerpt":{"rendered":"<div class=\"overview-image\"><img decoding=\"async\"  class=\"img-responsive\" src=\"https:\/\/blogs.mathworks.com\/community\/files\/matlab3DPrinting_01.png\" onError=\"this.style.display ='none';\" \/><\/div>\n<p>As promised, Paul Kassebaum is back this week with an in-depth discussion of how to get from a mathematical object in MATLAB to a solid object you can hold in your hand. Paul is a maker in the truest&#8230; <a class=\"read-more\" href=\"https:\/\/blogs.mathworks.com\/community\/2013\/06\/20\/paul-prints-the-l-shaped-membrane\/\">read more >><\/a><\/p>\n","protected":false},"author":69,"featured_media":0,"comment_status":"open","ping_status":"closed","sticky":false,"template":"","format":"standard","meta":[],"categories":[256],"tags":[],"_links":{"self":[{"href":"https:\/\/blogs.mathworks.com\/community\/wp-json\/wp\/v2\/posts\/2275"}],"collection":[{"href":"https:\/\/blogs.mathworks.com\/community\/wp-json\/wp\/v2\/posts"}],"about":[{"href":"https:\/\/blogs.mathworks.com\/community\/wp-json\/wp\/v2\/types\/post"}],"author":[{"embeddable":true,"href":"https:\/\/blogs.mathworks.com\/community\/wp-json\/wp\/v2\/users\/69"}],"replies":[{"embeddable":true,"href":"https:\/\/blogs.mathworks.com\/community\/wp-json\/wp\/v2\/comments?post=2275"}],"version-history":[{"count":20,"href":"https:\/\/blogs.mathworks.com\/community\/wp-json\/wp\/v2\/posts\/2275\/revisions"}],"predecessor-version":[{"id":8109,"href":"https:\/\/blogs.mathworks.com\/community\/wp-json\/wp\/v2\/posts\/2275\/revisions\/8109"}],"wp:attachment":[{"href":"https:\/\/blogs.mathworks.com\/community\/wp-json\/wp\/v2\/media?parent=2275"}],"wp:term":[{"taxonomy":"category","embeddable":true,"href":"https:\/\/blogs.mathworks.com\/community\/wp-json\/wp\/v2\/categories?post=2275"},{"taxonomy":"post_tag","embeddable":true,"href":"https:\/\/blogs.mathworks.com\/community\/wp-json\/wp\/v2\/tags?post=2275"}],"curies":[{"name":"wp","href":"https:\/\/api.w.org\/{rel}","templated":true}]}}