# MATLAB Spoken Here

## MATLAB EXPO UK 2013 – Something for Everyone

Guest post from Dr. Tanya Morton.

Last month, MathWorks held MATLAB EXPO UK at the Silverstone Circuit, the home of British motor racing. It was a great opportunity for the MATLAB and Simulink community to come together and the event provided something for everyone across multiple industries, applications and level of experience. The background roar of cars racing round the track certainly added to the atmosphere of the event.

A photo montage from MATLAB EXPO

The breadth of use of MATLAB in industry and academia never fails to impress me. This event featured top customer speakers from the aerospace, automotive, power, finance, and education industries:

• Rolls-Royce discussed simulation and modelling for aerospace engines, including how MATLAB helped to model the effects of volcanic ash, which caused so much disruption in 2010
• Ford Motor Company highlighted the use of MATLAB to reduce carbon dioxide emissions in automotive vehicles
• HSBC presented a policy simulator for modelling global financial risk, built using object oriented MATLAB
• University of Cambridge demonstrated their MATLAB package for enhancing mathematical education

Rolls-Royce takes the stage

The range of applications and technology showcased at the event was simply fantastic. Some of my personal favourites were:

• The full-size model of the BLOODHOUND SuperSonic Car in the exhibition hall with a real-time demo of the control system
• The NAO humanoid robot which can mimic human motion. Delegates had the opportunity to literally dance with a robot!
• Jasmina Lazic, an application engineer at MathWorks, gave a talk on Machine Learning for data analytics with a classification example. The session introduced delegates to when, why and how to select from a variety of machine learning methods
• Ellie Dobson, application engineer at MathWorks, talked the audience through techniques to speed up MATLAB algorithms through code analysis, multicore CPUs, GPUs, cloud technologies and MATLAB Coder
• Marta Wilczkowiak, application engineer at MathWorks, discussed the myriad of ways of sharing your MATLAB algorithms including via Apps and Web Services

A NAO Humanoid Robot takes a break

This year we added some introductory sessions for delegates wanting a primer on MathWorks core technology. These sessions proved very popular and we plan to expand them for next year. There was plenty for the advanced user too with a great ‘What’s new in MATLAB?’ talk from Chief Scientist Joe Hicklin, and masterclass tracks for both MATLAB and Simulink.

An animated Chief Scientist makes his point

There was also plenty of geeky fun on offer at MATLAB EXPO through games and competitions. In particular,

• MATLAB EXPO Conversations offered a prize to the best networkers by analysis of attendee interactions using active RFID tags with Arduino chips, Raspberry Pi base stations, and a web front-end powered by MATLAB Production Server – a riot of cool technology!
• The magnetic MATLAB puzzle boards challenged attendees to solve a Cody checkerboard challenge problem
• The Model-Based Design track took a fun Simulink quiz with 13 multiple-choice questions (and yes, the answer to all of them was b to celebrate the R2013b release!)
• Lucky attendees were able to try to get the fastest lap-time on the F1 Simulator and competition was fierce!

Concentrating hard in the F1 Simulator

The presentations from MATLAB EXPO are now available for download from the website:
http://www.matlabexpo.com/uk/

I hope to see some of you at next year’s community event!

## MATLAB Mobile Autocompletes You

With today’s release of MATLAB Mobile, typing commands just got easier. With the new autocomplete feature, MATLAB Mobile offers suggestions for functions, variables, methods, properties, files and folders.

To see this in practice, type two or more characters in the command window:

Suggestions for MATLAB Functions (iPhone)

Tap on an option to complete it to the window.

You can leverage autocomplete for a variety of tasks, like navigating folder hierarchy:

Navigating Folder Hierarchy (Android)

Or displaying the methods that an object exposes:

Suggestions for an Object’s Methods (iPhone)

Autocomplete is available on the latest releases of MATLAB Mobile on the App Store and Google Play.

Has autocomplete improved your experience? Let us know by leaving a comment below.

## MathWorks Support Solutions in MATLAB Answers

MathWorks-authored technical support solutions are now part of Answers in MATLAB Central.

In the past MATLAB users had to go to two different places to find questions and answers.  MathWorks-authored technical support solutions were located in the support area on the MathWorks website, while community questions and answers were located in Answers on MATLAB Central.  Now both are located in MATLAB Answers. This creates a one-stop shop for MATLAB users to find answers to their questions, making the process faster and easier.

At first glance, these appear to be questions and answers with the same contributor.  In reality these Support Answers are questions that were previously asked to MathWorks Technical Support by customers.  The “MathWorks Support Team” contributor name will appear as both the questioner and answerer, even though the question was originally asked by someone else.

A lot of Support Answers were recently added to MATLAB Answers, and we want to make sure it is still easy for Community members to find what they are looking for.  There is a toggle filter located in the left navigation panel titled “Refine by Source.” Clicking on one of the two options will make search results display either Community or MathWorks Support contributions.

Another way to differentiate between Community and MathWorks Support content is this little blue MathWorks icon located in the top right corner of all MathWorks Support contributions:

Below we show an example of a simple product download question that MATLAB users have asked.  The three most relevant solutions found in MATLAB Answers are shown for “How to download older version of MATLAB”

As you can see, the first and third questions with accepted answers is from the MathWorks Support Team, and the second question is from someone within the MATLAB Central Community.  If you drill into the top question the next screen will show the question previously asked by a MATLAB user, and the answer MathWorks Support provided the user.  Any additional Community comments and answers will also be visible on this page.

We think this will make MATLAB Answers even more useful as a resource for MathWorks-related inquiries.

## September 24th, 2013

If you’ve been following Cody for any time at all, you know that we have some superstar competitors who have stayed in the top ranks since the very beginning. Some time ago I started tracking the top 10 players using Trendy. Here’s what that chart looks like.

For the longest time, and dating back to the day Cody launched, Bryant Tran (@bmtran) was the undisputed king of the hill. We even featured him in an interview on this blog. Eventually Richard Zapor was to surpass him. Richard is a tireless competitor and has made something of a specialty of creating high-quality challenges for the community (173 and counting!).

If you look at the blue line that represents Richard, you can see that he has steadily moved upward. But now look at Alfonso Nieto-Castanon. His green line progresses in enormous leaps punctuated by periods of relative inactivity. Alfonso finally passed Bryant Tran earlier this year, and today I was surprised to see that he was also converging rapidly with Richard at the number one spot.

When I saw this graph, I dashed straight over to the Cody leaderboard and what should I see but an exact tie for first place between Richard and Alfonso!

I was so impressed I took this snapshot. It’s probably changed by now and who knows who will be in the lead at the end of the day? But still, this snapshot depicts a great moment in Cody and sporting history.

Well done, gentlemen! It’s an honor to compete in the same game, though not in the same league.

## September 18th, 2013

Sometimes it takes more than words to explain a question or demonstrate an answer. And starting today it’s much easier to enhance your questions and answers by inserting images and attaching files.

MATLAB Markup allows you to display images and link to files. But those images and files have to be hosted somewhere on the web. In the past, that meant finding a place to host your file, figuring out the URL to the file, and adding the URL to your markup. And then hoping that the hosting service does not  disable or remove the file. But now you can host your files directly on MATLAB Answers!

Suppose I have an image that I want to use in my question. It’s saved on my computer. I’d like to add it to my question (or answer–the process is the same). All I have to do is click on the image icon in the markup toolbar, choose the image from my computer, and insert it. We will upload the image in the background and add the correct markup:

Attaching files works the same way. However, when you attach a file we don’t add any markup. Instead, a link to the file is added below the text of your question or answer (think of the way that email attachments work).

## The old way still works

If you’d like to use an image hosted somewhere else, that’s cool. We don’t mind. In fact, now you have two ways to include it:

1. Directly type the markup. To insert an image all you have to do is include the image URL in double angle brackets. For example, <<http://mathworks.com/matlabcentral/images/surf.gif>>.
(Note: Make sure that there is a blank line before and after your image markup.)
2. Use the new insert image button, but select “From the Web.” Type the image URL and we’ll insert the markup.

We hope that these new tools make it easier for you to ask and answer questions!

## Two Hundred Creators and Counting

About six weeks ago John Kelly wrote a post about some remarkable achievements in Cody. In that post he offered a prize to whomever would win the 200th creator badge.

I am happy to report that the prize has been claimed by Kevin Hellemans. You may remember that last Friday was the 13th day of September. Appropriately, our friend Kevin created problem 1877, entitled Friday the 13th. So we might say that Friday the 13th was Kevin’s lucky day.

An interesting footnote to this story is how we detected that Kevin was the 200th winner of the Creator Badge. Did we visit the Creator Badge page every morning to see if the number was greater than 200? No! That would be tedious, and besides there’s a much nicer way to do it. For this, we invoke another MATLAB Central application: Trendy.

What I did was create a trend to keep track of the Number of People Awarded the “Creator” Badge. Then I set a trigger so that I would get an email when the number collected was greater than 199.

This morning when I got to work, there was an email message from Trendy telling me that the magic number had been reached. And that is the story of how Kevin Hellemans won his prize.

So, thanks Kevin! John is sending a prize your way.

## Formula Student Germany – A Design Competition

This week we bring back guest blogger Joachim Schlosser, who leads a team of Education Technical Specialists from his office in Germany. That gave him a special vantage point for a thrilling student competition.

## Formula Student Germany

by Joachim Schlosser

The air is shimmering at 37°C. From a cloudless sky, the sun beats down on the asphalt track. We are at an automobile race track in Hockenheim, in southwest Germany. Numerous groups of students wearing matching t-shirts crowd around the cars the pits. It is a young crowd, with students specializing in mechanical engineering, electrical engineering, mechatronics, business and finance, and computer science. Among them is Julian Vossen, a student of mechanical engineering at RWTH Aachen University. He is examining the data from the recent run of his car’s cooling system. His team, Ecurie Aix Aachen, built a car, and they came here to win. The event is Formula Student Germany.

Formula Student Germany brings together 3000 students and 115 teams from 30 nations all over the world. Modeled on a similar event in North America, it has been running since 2006. The event features two different kinds of cars: conventional combustion engines (75 teams) and purely electrical vehicles (40 teams). Over the last year these students have sacrificed evenings, weekends, and sometimes the attendance of university seminars and classes to finish their project.

What is so special about Formula Student Germany? Both Tim Hannig, the chairman of the board, and Julian highlight the special atmosphere. There is plenty of competition, but the teams are incredibly helpful to each other when something breaks and another team needs something. And all agree that it is a really good feeling to build a car from scratch. They develop a concept for the chassis and the drive train. They have to think about how many and which electric drives to use. They develop gearboxes. They evaluate different concepts for cooling. Some teams even build the car body out of composite carbon fiber materials, which requires a huge effort.

So they design and construct a real car. And therefore, they want to use the software tools that car manufacturers and their suppliers like for real car design, too. More and more teams were approaching MathWorks, and so the company now offers complimentary use of MATLAB & Simulink and more than 35 add-ons to all teams participating in Formula Student Germany.

## Enter the cooling man

Back to Julian Vossen. In the team from RWTH Aachen, he is responsible for building the cooling system of the engine. The engine generates a lot of heat that has somehow to be transferred away. Otherwise it could lead to an engine meltdown. So the cooling system is important for the overall performance. The lower the temperature of the engine, the better its efficiency. The problem ultimately is to control the temperature over the course of a 22-kilometer endurance test.

Photo courtesy Formula Student Germany, Kimmo Hirvonen.

Julian modeled the whole system in Simulink and Simscape, creating a block diagram that represents the physical parts of the system like pumps, motors, heat sinks. He then fed the data about the race course that he had into the system and so could simulate the temperature load of the whole race without having built any real part of the car yet. Julian saw the different curves of temperature and used these results to dimension the radiator and other parts of the cooling system.

About MATLAB & Simulink, Julian says: “It was great to have a tool in which you can simulate transients. We can see for every second in the race temperature and other parameters.” In Simscape, the first principles add-on for physical modeling, he chooses the cooler size and pump size that exactly fit the needs of his system. The result? Their car survived the most challenging race: the endurance test, where you have to keep the car running for 22 km. They even made a pretty good lap time and ranking, which takes into account the energy consumption.

So cooling the engines is important for generating lots of torque. But how do you actually bring that torque on the pavement?

## Measuring and simulating vehicle dynamics

Since 2010 the team of the University in Bayreuth competes in the Formula Student Germany Electric series. They have a pretty good electric drive pack, so the big challenge is to build a chassis that brings the enormous torque to the wheels and on the street. This highly depends on how you lay out the system and how exactly the forces act on the different poles of the wheel suspension.

Photo courtesy Formula Student Germany, S. Schindler.

So the task is to get some measurement data from the mechanics that allows calculating the individual force vectors. Elias Schmidek is student of mechanical engineering and team leader of vehicle dynamics at Elefant Racing. He applied strain gauge strips at all the 12 poles that have contact to the wheels—3 per wheel— to measure the tiny stretching that appears during the curves on the track. Each round at the test track produces tons of log files, including a trigger signal to get the time and the voltage. Furthermore, he logs the steering angle and yaw rate. All that data feeds into MATLAB, where Elias’ script calculates the forces’ absolute values for each of the connection points of the poles—the joints. Since he knows the characteristics of the materials used for the poles, he can generate the proper equations in MATLAB.

Now that he as the absolute values for all the twelve forces, it is time to calculate the force vectors. For this, Elias and his teammate Bernhard Hauffen created a model of the wheel suspensions in SimMechanics, the multi-body simulator built on the Simscape technology in Simulink. If you think of the suspension as several poles, connected by rotating joints and a prismatic joint to allow explicitly for the vertical spring damper movement, you get an accurate model of the real thing if you just design every length and position of a joint accurately. In SimMechanics, Elias feeds the absolute values for each of the joints into the simulation. The simulation now exactly puts out the three dimensional vector of all joint movements in time, and of course also the three dimensional vectors for the movement of the four wheels.

For every test on the track, the calculation of the force vectors takes just a couple of minutes, and ultimately gives insight into the forces on the wheels. With this information Elias, Bernhard and their team can construct and calibrate their chassis. And since they have all lengths of the rods as parameters in the MATLAB workspace, they can easily simulate different possibilities in their Simulink and SimMechanics model for further optimization, and eventually apply numerical optimization methods that are available in MATLAB via the Optimization Toolbox.

Why spend the effort to build such a simulation? As Elias says, having proper data on the car is the key. A drivetrain, as important as it is to be as powerful as possible, is only worth as much as you can get the power on the street. If the wheels do not have proper contact to the ground, power is lost by wheel spin. Consider the acceleration test: from a standsill the cars have to cover 75 meters as fast as possible. The best cars do this in less than four seconds, and the car from Bayreuth is among these, scoring 7th of the 40 teams.

The sheer breadth of applications of MATLAB & Simulink and the add-ons we were able to see during the week and discuss with lots of engaged members of many different teams is just amazing. What these students create in just a couple of months is the real thing.

We at MathWorks are proud to help the teams being more successful by donating a fairly large software package around MATLAB & Simulink to all teams.

The Formula Student Germany 2013 ended Sunday night with the awards ceremony for both the combustion competition as well as the electric competition. In Formula Student Combustion, Global Formula Racing (@TeamGFR) won ahead of Rennteam Stuttgart (@Rennteam0711) and Rennstall Esslingen (@Rennstall). In Formula Student Electric, DUT Racing Team (@DUT_Racing) scored 1st ahead of AMZ Racing from ETH Zürich (@amzracing) and ka.race.ing from KIT (@KA_RaceIng). Congratulations from MathWorks to those and to all other teams that made it to Hockenheim!

## Cody Recognition

Cody has recently reached several significant milestones. I’d like to take a moment to recognize these milestones and also offer a challenge.

First let’s acknowledge two people: the winner of the 5,000th Solver badge, and the author of 1,000th problem.

Badges in Cody are the Olympic medals of MATLAB achievement. Well, maybe not that grand, but as players contribute to Cody, badges and points are earned in recognition of various accomplishments. A Cody player’s reputation increases with each badge. David (that’s the only name he’s made public) is the 5000th person to solve a problem on Cody. This achievement also makes him the 5000th person to earn the Solver Badge. David solved Problem 1. Times 2 – START HERE, which is the first problem we recommend to all new Cody Players.

There are currently more than 1100 problems on Cody. Who created the 1000th? Counting problems is a little tricky, because some have been deleted along the way. You might think the problem with ID=1000 (Richard Zapor’s Zernike Coefficients) would be the right choice. But that was actually only the 698th problem in order of creation. The prize for Number One Thousand goes to Claudio Gelmi for his problem Number of lattice points within a circle.

Cody is a game, and like most games, there is a score. The scores change all the time, but let’s recognize the amazing achievements of our current top three. We can call these our current virtual Gold, Silver, and Bronze medal winners. But watch closely… in another two days the top three might shift around again. Alfonso is charging hard!

Finally, here’s the challenge. The next big milestone we want to recognize is the 200th person to receive the Creator Badge. Just as a store might honor its one millionth customer, MathWorks will honor the 200th player to receive the Creator Badge by giving out a small prize. What will the prize be? How small will it be? Maybe it will be a t-shirt, maybe a water bottle, or maybe even a small electric car*. I will have to check to see what is in the MathWorks prize vault.

*Editor’s Note: Prize will not be a car

The three virtual medal winners (Richard Zapor, @bmtran, and Alfonso Nieto-Castanon) as well as David and Claudio will be sent a MathWorks prize as a thank-you for their participation on Cody.

As I mentioned I will also be sending a prize to the 200th community member to receive their Creator Badge. If you haven’t tried Cody yet, now is the time. The next prize may be yours!

## Anamorphic 3D Printing

You’ve seen funny mirrors at carnivals, the curvy ones that make your head look tiny and your legs impossibly long. Now imagine a person who looks strange in real life, say with a giant head and extremely short legs. Put this unusual character in front of the same curvy mirror and now you see someone who looks normal. That’s the idea behind an anamorphic representation. In his follow-up to the last post, Paul Kassebaum shows us how to compute the inverse funhouse mirror effect for our old friend, the L-shaped membrane.

## Reflections in a Cylinder

by Paul Kassebaum

Using MATLAB to make models for 3D printing gives you the ability to perform all sorts of mathematical transformations to define or deform your creation.

Anamorphic drawing by Istvan Orosz, entitled Anamorph with column 2.

In my previous post, we explored how to turn equations and data into 3D prints, using the L-shaped membrane as an example. In this post, we’ll have a little physics fun with our model in MATLAB before sending it off to be printed, deforming it to create an anamorphic sculpture. In particular, we’ll create a 3D print of the L-shaped membrane that is meant to be viewed in the reflection of a cylinder.

### Contents

#### Underlying physics

To make an anamorphic sculpture look intelligible when reflected off a curved mirror, the only physics you need is the law of reflection. The law of reflection notes that if a ray of light hits a mirror at an angle measured off a line perpendicular to the mirror, then the reflected ray of light will come off the mirror at that same angle.

Simple case of the law of reflection where the mirror is flat.

Reflections off flat mirrors create every-day optical illusions. The mirror appears to be a window to an ethereal world on its other side. However, every ray of light that appears to be coming from the other side of the mirror is in fact coming from the same side we are standing on, only bent at angles that we can calculate using the law of reflection.

To make our anamorphic sculpture work with a cylindrical mirror, we have to imagine that the L-shaped membrane is sitting on the other side of the window the mirror appears to make. This is the same as imagining that the sculpture is sitting in the center of the cylinder. Let me explain further by referencing the drawing below. For every point P on this apparent sculpture inside the cylinder, we’ll draw a straight line to the eye of the viewer at a point V. This line from P to V is the path that our eyes think a ray of light took. However, the light actually came from some point P’ on our side of the mirror, hit the mirror at some angle, reflected off the mirror at that same angle, and then hit our eyes.

Drawing of how light reflects off a cylinder from P’ to the viewer’s eye at V, creating the apparent image at P. The cylinder is shown from above.

Our task is to figure out where P’ is located given the location of P, V, and the radius of the cylinder. Since our sculpture of the L-shaped membrane is made up of triangular faces, we need to solve this problem for every triangle vertex.

Once the intersection points are found, the problem is essentially two dimensional. Each reflected pair of points P and P’ will have the same altitude. The rest of the problem focuses on lines and points that are projected onto a plane to create drawings like the one shown above.

#### The inverse reflection problem

This problem can be solved in just a few steps.

1. First, we’ll find the point of intersection C between the line PV and the cylinder.
2. Then we’ll find the angle between the line perpendicular to the cylinder and the line PV. This is the angle of the reflected ray of light shown in the drawing above.
3. Finally, we know that the light ray CP should have the same length as CP’, so we can find P’ by rotating the point P about C by an angle that satisfies the law of reflection.

#### The intersection points

Imagine a photon traveling along the line PV from P to V. This photon’s journey can be tracked by a vector J defined by position vectors P and V, and a factor t:

J = P + t*(V-P) where t is in the range [0,1],

so as t increases from 0 to 1, the photon moves from P to V.

The photon will hit the cylinder’s surface at some particular t. At this point, the position vector of the photon measured from the cylinder’s center will have the same length as the radius of the cylinder, or

$|J|^2 = R^2$

Demonstration of the equality above.

This equality sets up a quadratic equation for the value of t at which J intersects the cylinder. This equation can be solved for t using the quadratic formula

A = p(:,1).^2 + p(:,2).^2 - 2*p(:,1)*v(1)...
+ v(1)^2    + v(2)^2    - 2*p(:,2)*v(2);
B = 2*(p(:,1)*v(1) + p(:,2)*v(2) - p(:,1).^2 - p(:,2).^2);
C = p(:,1).^2 + p(:,2).^2 - R^2;
t = (-B + sqrt(B.^2 - 4*A.*C))./(2*A);


The code above solves the problem for every point on the apparent model. The variable p is an n-by-3 matrix representing n points of the apparent model’s triangular vertices, each with 3 cartesian coordinates, and v is a 3 component vector representing the viewing point. In this way, t is now an n component vector with each component of t associated with each point in the matrix p. Notice that although the quadratic formula gives us two roots, we’ve chosen the one that guarantees that t is always positive.

The intersection points c are then

pv = bsxfun(@minus,v,p); % p to v
c  = p + bsxfun(@times,t,pv);


#### The angles of reflection

Now that we know the point of intersection C for each point on the apparent model P, we can find the angles of reflection. Since we would like to know both the angle and the orientation (clockwise or counterclockwise) to rotate the point P about C, we should use the ATAN2 function to calculate the angle between the position vector for C and the vector PC from P to C.

pc = c - p;
theta = atan2(c(:,2),c(:,1)) - atan2(pc(:,2),pc(:,1));


#### Locate the real model’s points

With the reflection angles in hand, we can now specify exactly where the points P’ must be. We’ll rotate the vector PC about C by twice the reflection angle: once to bring it along the vector to C, and once more in accordance with the law of reflection.

For each apparent point P, we will have a matrix equation

$p’ = R(\theta) * (c – p) + c,$

where R is a rotation matrix.

To locate P’, draw position vectors to P and C, create the vector PC, rotate this vector by twice the angle calculated above, then translate the vector along C.

In code, we’ll first set up all of the rotation matrices

theta = 2*theta;
Rotate = [cos(theta) -sin(theta) zeros(numel(theta),1),...
sin(theta)  cos(theta) zeros(numel(theta),3),...
-ones(numel(theta),1)];% negative so
% that pp_z = p_z
Rotate = reshape(Rotate',3,3,[]); % 3-by-3-by-n


Then we’ll reshape our vectors and perform the rotations and translations.

p = reshape(p',3,1,[]); % 3-by-1-by-n
c = reshape(c',3,1,[]); % 3-by-1-by-n
pp = zeros(size(p)); % preallocate p'
for k = 1:size(p,3)
pp(:,:,k) = Rotate(:,:,k)*(c(:,:,k) - p(:,:,k)) + c(:,:,k);
end
% Reshape pp to bring it into the original form of p
pp = reshape(pp,3,[])'; % n-by-3


#### The big picture

Putting all of those steps together, we can define the following function

function pp = reflectThroughCylinder(R, v, p)
end


Let’s put this function to use on the L-shaped membrane sculpture! Using the code explained in the previous post, we have the L-shaped membrane sculpture defined by the vertices of triangles stored in the variable shellVertices, and the way these vertices are connected stored in the variable shellFaces.

R = 2*sqrt(2)/2; % Twice the minimum radius to fit the apparent sculpture.
v = 1.5*R*[cos(0.9406),sin(0.9406),1]; % View from 1.5 radii away
% creates a nice distortion.
p = bsxfun(@minus,shellVertices,[0.5,0.5,0]);% Place origin at x-y
% center of sculpture.
pp = reflectThroughCylinder(R,v,p);

% Visualize the anamorphic sculpture
trisurf(shellFaces,...
pp(:,1),...
pp(:,2),...
pp(:,3));
colormap pink;
axis equal;
axis vis3d;


The apparent and reflected L-shaped membranes.

Here’s the final product!

## Paul Prints the L-Shaped Membrane

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.

## 3D Printing from MATLAB

by Paul Kassebaum

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.

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.

### Contents

#### Stereo-lithography (STL) Files

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:

facet normal nx ny nz
outer loop
vertex v1x v1y v1z
vertex v2x v2y v2z
vertex v3x v3y v3z
endloop
endfacet

In practice, it might look like this:

facet normal 6.6998060E-01 -6.6246430E-01 3.3506277E-01
outer loop
vertex 7.9166667E-01 2.5000000E-01 9.4328269E-01
vertex 7.5000000E-01 2.5000000E-01 1.0265980E+00
vertex 7.9166667E-01 2.0833333E-01 8.6090207E-01
endloop
endfacet

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.

#### Printing Surface Plots

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.

n = 12; % number of partitions in each dimension.
[X,Y] = meshgrid(linspace(0,1,2*n+1));
L = (40/51/0.9)*membrane(1,n);
surf(X,Y,L);
colormap pink;
set(gca,'dataAspectRatio', [1 1 1]);


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.

faces   = delaunay(X,Y);
patches = trisurf(faces,X,Y,L);
set(gca,'dataAspectRatio', [1 1 1]);


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.

Let us begin by creating a 3-by-3-by-n tensor called ‘facets’ 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.

vertices = get(patches,'vertices');
facets = vertices';
facets = reshape(facets(:,faces'), 3, 3, []);


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.

% SQUEEZE compacts empty dimensions.
edge1 = squeeze(facets(:,2,:) - facets(:,1,:));
edge2 = squeeze(facets(:,3,:) - facets(:,1,:));
normals = edge1([2 3 1],:) .* edge2([3 1 2],:)...
- edge2([2 3 1],:) .* edge1([3 1 2],:);
normals = bsxfun(@times,...
normals, 1 ./ sqrt(sum(normals .* normals, 1)));


We can compute the normal at each vertex by averaging the normals of its neighboring facets.

meanNormal = zeros(3,length(vertices)); % pre-allocate memory.
for k = 1:length(vertices)
% Find all faces shared by a vertex
[sharedFaces,~] = find(faces == k);
% Compute the mean normal of all faces shared by a vertex
meanNormal(:,k) = mean(normals(:,sharedFaces),2);
end
meanNormal = bsxfun(@times, meanNormal,...
1 ./ sqrt(sum(meanNormal .* meanNormal, 1)));


Now we merely need to copy and shift all of the vertices downward along their normals. we will call these new vertices ‘underVertices’.

shellThickness = 0.05;
underVertices  = vertices - shellThickness*meanNormal';


Let us see what we have got so far.

underFaces = delaunay(underVertices(:,1),underVertices(:,2));
trisurf(underFaces,...
underVertices(:,1),...
underVertices(:,2),...
underVertices(:,3));
set(gca,'dataAspectRatio', [1 1 1],...
'xLim',[0 1],'yLim',[0 1]);


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’s DELAUNAYTRI class allows us to specify the boundary in 2D, which solves our problem.

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.

boundaryIndices = ...
[find(vertices(:,2) == min(vertices(:,2))); % min y
find(vertices(:,1) == max(vertices(:,1))); % max x
find(vertices(:,2) == max(vertices(:,2))); % max y
find(vertices(:,1) == min(vertices(:,1)))];% min x


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.

boundaryIndices = [...
boundaryIndices(1:floor(end/4-1)); % semi open interval [1, end/4).
boundaryIndices(floor(end/4+1):floor(end/2));%[end/4, end/2)
boundaryIndices(floor(end*3/4-1):-1:floor(end/2+1));%[end/2,end*3/4)
boundaryIndices(end-1:-1:floor(end*3/4+1))]; %[end*3/4, end)


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.

constrainedEdges = [boundaryIndices(1:end-1), boundaryIndices(2:end)];
underFaces = DelaunayTri(...
[underVertices(:,1),underVertices(:,2)],constrainedEdges);


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.

inside = underFaces.inOutStatus; % 1 = in, 0 = out.
underFaces = underFaces.Triangulation(inside,:);


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’s normals.

underFaces = fliplr(underFaces);


Let us see what we have got so far.

trisurf(underFaces,...
underVertices(:,1),...
underVertices(:,2),...
underVertices(:,3));
set(gca,'dataAspectRatio', [1 1 1],...
'xLim',[0 1],'yLim',[0 1]);


This looks much better. Now we can move on to connect these two surfaces with a third surface that we will call ‘wall’. 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.

wallVertices = [vertices(boundaryIndices,:);
underVertices(boundaryIndices,:)];
% Number of wall vertices on each surface (nwv).
nwv          = length(wallVertices)/2;
% Allocate memory for wallFaces.
wallFaces    = zeros(2*(nwv-1),3);
% Define the faces.
for k = 1:nwv-1
wallFaces(k      ,:) = [k+1  ,k      ,k+nwv];
wallFaces(k+nwv-1,:) = [k+nwv,k+1+nwv,k+1];
end


Let us see what we have got so far.

trisurf(wallFaces,...
wallVertices(:,1),...
wallVertices(:,2),...
wallVertices(:,3));
set(gca,'dataAspectRatio', [1 1 1],...
'xLim',[0 1],'yLim',[0 1]);


Now let us assemble our three surfaces into one closed surface that we will call ‘shell’.

% Shift indices to concatenate with the original surface.
underFaces = underFaces +   length(vertices);
wallFaces  = wallFaces  + 2*length(vertices);
% Concatenate the results.
shellVertices = [vertices; underVertices; wallVertices];
shellFaces    = [faces;    underFaces;    wallFaces];


Most 3D printers require that all of the vertex coordinates be non-negative, so we will shift our shell up to satisfy this convention.

minZ = min(shellVertices(:,3));
shellVertices = shellVertices...
- repmat([0 0 minZ],length(shellVertices),1);


Let us look at the final result.

trisurf(shellFaces,...
shellVertices(:,1),...
shellVertices(:,2),...
shellVertices(:,3));
set(gca,'dataAspectRatio', [1 1 1],...
'xLim',[0 1],'yLim',[0 1]);


The final stage is to convert this surface into an STL-file. We refer to the excellent entry in the File Exchange called stlwrite 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.

% Name your STL-file
filename = 'MathWorksLogo.stl';

% Create the facets.
shellFacets = shellVertices';
shellFacets = reshape(shellFacets(:,shellFaces'), 3, 3, []);

% Compute their normals.
edge1 = squeeze(shellFacets(:,2,:) - shellFacets(:,1,:));
edge2 = squeeze(shellFacets(:,3,:) - shellFacets(:,1,:));
shellNormals = edge1([2 3 1],:) .* edge2([3 1 2],:)...
- edge2([2 3 1],:) .* edge1([3 1 2],:);
shellNormals = bsxfun(@times,...
shellNormals, 1 ./ sqrt(sum(shellNormals .* shellNormals, 1)));

% Associate each facet with its normal.
shellFacets = cat(2, reshape(shellNormals, 3, 1, []), shellFacets);

% Open the file for writing
fid = fopen(filename,'wb+');

% Write the file contents.
fprintf(fid,'solid %s\r\n',filename);
% Write DATA.
fprintf(fid,[...
'facet normal %.7E %.7E %.7E\r\n' ...
'outer loop\r\n' ...
'vertex %.7E %.7E %.7E\r\n' ...
'vertex %.7E %.7E %.7E\r\n' ...
'vertex %.7E %.7E %.7E\r\n' ...
'endloop\r\n' ...
'endfacet\r\n'], shellFacets);
% Write FOOTER.
fprintf(fid,'endsolid %s\r\n',filename);

% Close the file.
fclose(fid);


#### MATLAB Materialized

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 Shapeways, Ponoko, Sculpteo, I.Materialize, ZoomRP, or RedEye. I am a member of the maker-space Artisan’s Asylum, which has a nice 3D printer, among many other tools. Here is a photo of my completed L-shaped membrane.

Good luck and have fun turning your MATLAB plots into tangible objects!

News from the intersection of MATLAB, Community, and the web.

These postings are the author's and don't necessarily represent the opinions of MathWorks.