Brett's Pick this week is Skeleton3D, by Philip Kollmansberger.A couple years ago, I put together a demo to show how to use MATLAB to calculate the tortuosity, or “twistedness,” of blood vessels. There are different formulas for measuring that property, but perhaps the easiest way of determining tortuosity of a... read more >>

]]>Brett's Pick this week is `Skeleton3D`, by Philip Kollmansberger.

A couple years ago, I put together a demo to show how to use MATLAB to calculate the tortuosity, or “twistedness,” of blood vessels. There are different formulas for measuring that property, but perhaps the easiest way of determining tortuosity of a path is using the arc-chord ratio—i.e., the length of the curve divided by the distance between its endpoints.

In my demo, I assumed the vessels were 2-dimensional, and then used some Image Processing Toolbox functionality to: 1) skeletonize the vessels of interest using `bwmorph`/endpoints, and 2) determine the endpoints of the segment (facilitated by a call to `bwmorph`/endpoints. Then, using one or both endpoints of the segment, respectively, to 3) calculate the geodesic distance transform of the vessel (`bwdistgeodesic`), the maximum value of which provides the arc-length; and 4) calculate the distance the endpoints. (You can find a discussion of this approach in a recorded webinar at mathworks.com.)

This approach is reasonably straightforward when you assume planarity, but in 3D, it gets much more challenging—because bwmorph doesn’t currently support the analysis of 3D images. That leaves us struggling to do the skeletonization- and endpoint-detection steps.

Enter Philip’s `Skeleton3D` function. Using an algorithm (properly cited in the code) based on “3-D medial surface/axis thinning algorithms,” `Skeleton3D` finds axial values that can be subsequently analyzed for the determination of tortuosity.

load testvol inds = find(testvol); [y,x,z] = ind2sub(size(testvol),inds); s = scatter3(x,y,z,100,... 'Marker','o',... 'MarkerEdgeColor','none',... 'MarkerFaceAlpha',0.3,... 'MarkerFaceColor','c'); xlabel('x'); ylabel('y'); zlabel('z'); title('Test Volume')

Now we can extract and display points along the skeleton of this test volume:

skel = Skeleton3D(testvol); inds = find(skel); [yskel,xskel,zskel] = ind2sub(size(skel),inds); hold on p2 = scatter3(xskel,yskel,zskel,20,... 'Marker','o',... 'MarkerEdgeColor','none',... 'MarkerFaceColor','r');

For tortuosity determinations, there's still some more work to be done--we still need to detect endpoints and arclengths, for instance. Those are topics for a different day.

For Philip, an implementation comment: `Skeleton3D` is considerably slower than it needs to be. It relies on looping over indices where vectorized approaches can be much faster. And converting to logical indexing for comparisons can also speed things up. For instance:

In the determination of "candidate foreground voxels," you can replace:

cands = intersect(find(cands(:)==1),find(skel(:)==1));

with

cands = find(cands & skel);

(Admittedly, that buys just a little improvement in performance, but I find it more elegant.)

More pointedly, if you were to profile your code, you'd find a large percentage of time spent calculating the skeleton of your testvol reflects nearly 87,000 calls to `p_oct_label`. That code could be greatly sped up by replacing successive indx(row) = `find(cube(row,:)==1)` calls in the analysis of 'cube' with something along the lines of

[row,col] = ind2sub(size(cube),cube==1); indx(row) = col(row == 1);

and iterating over values of `row`.

Nonetheless, this provides a good start to anyone needing to skeletonize structures in 3D. Thanks, Philip, for this very nice contribution!

As always, I welcome your thoughts and comments.

Get
the MATLAB code

Published with MATLAB® R2016a

Sean's pick this week is GetAuthentication by Stefan. I was recently asked about a simple way to password protect an executable compiled with MATLAB Compiler. The customer's goal was... read more >>

]]>Sean's pick this week is GetAuthentication by Stefan.

I was recently asked about a simple way to password protect an executable compiled with MATLAB Compiler. The customer's goal was that this would be on a shared machine and only authorized users would be able to use the executable.

It occurred to me that I frequently use a MATLAB-based password box every time I upload a file to our FTP or SFTP site to share with customers. Digging into the code there, I discovered this file and was excited to see it publicly available on the File Exchange.

To use it to protect your executable, simply call it in a wrapper function to decide whether to open the main executable or quit. Here's a simple example that gives access to my four dogs to see a picture of themselves.

function authenticateApp() % This function authenticates a username and password before opening an % application. % Username and password for each of my dogs unames = {'Zoey', 'Mya', 'Lexi', 'Lilly'}; passwords = {'MilkB0ne', 'Dent@Stix', 'Pupp3roni', 's0ck$'}; % Try three times for ii = 1:3 % Get authentication [user, pass] = GetAuthentication(); % Does a username match? useridx = ismember(unames,user); if nnz(useridx) % If yes, does the password match for that user? if isequal(useridx, ismember(passwords,pass)) % Open application here: imshow('fourdogs.jpg') break end end end end

Calling it:

authenticateApp

Now if you want to compile this, make the authentication wrapper the main function and your app will have this protection around it.

My only suggestion would be for `GetAuthentication` to pass back a third output argument about whether the user hit the cancel button. Right now, if cancel is hit, it just
returns empty for password and username. This is ambiguous to hitting okay with empty password and username which is a challenge
to handle if you want to use it in a loop like I did.

Give it a try and let us know what you think here or leave a comment for Stefan.

Get
the MATLAB code

Published with MATLAB® R2016b

Jiro's pick this week is Volumetric Data Explorer by my fellow MathWorker, Adam Filion.Volume visualization refers to a graphical representation of data sets defined on three dimensional grids. There are two types of volume data sets that MATLAB typically deals with - scalar and vector data. Scalar data contain scalar... read more >>

]]>Jiro's pick this week is Volumetric Data Explorer by my fellow MathWorker, Adam Filion.

Volume visualization refers to a graphical representation of data sets defined on three dimensional grids. There are two types of volume data sets that MATLAB typically deals with - scalar and vector data. Scalar data contain scalar information at every three dimensional grid points, much like temperature.

$T = f(x,y,z)$

Vector data contain three values at each grid point, much like wind flow.

$[u,v,w] = f(x,y,z)$

There are many functions for visualizing volume data. Adam's Volumetric Data Explorer is an app that lets you interactively visualize scalar data, using some of these builtin functions under the hood. Here are some of the cool things about this app.

- The app lets you visualize orthogonal slices and isosurfaces of scalar data. Slices let you cut into your scalar data, and isosurfaces show you the surfaces with equal values, sort of like contour lines but for volume data.

- If you provide an additional dimension (time), the app will animate the data along the time.

- You can create recordings (AVI and MP4) of the animations. The app comes with a wealth of options for specifying animation frames, frame rates, recording regions, etc.

- The app comes with a detailed documentation, which you can access from the help menu.

I remember seeing one of the first versions of this app several years ago, and Adam has been adding feature after feature since then. Thanks for a great app, Adam!

**Comments**

Give this a try and let us know what you think here or leave a comment for Adam.

Get
the MATLAB code

Published with MATLAB® R2016a

Will's pick this week is Julia Sets (Live Editor Example) by Joshua Meyer. The MATLAB Live Editor is a wonderful, new medium for expressing your ideas in MATLAB. In this contribution, Joshua educates us on Julia sets. The Live Editor enables him to do more than merely supply code with commentary.... read more >>

]]>The MATLAB Live Editor is a wonderful, new medium for expressing your ideas in MATLAB. In this contribution, Joshua educates us on Julia sets. The Live Editor enables him to do more than merely supply code with commentary. We also receive a rich document in which Julia sets are clearly explained in a narrative. Equations are displayed as you'd find them in a textbook before getting written as MATLAB code to produce figures. Enjoy some of the exciting results!

Let us know what you think here or leave a comment for Joshua.

]]>

Sean's pick this week is Forbidden Donut by Jerome Briot. In honor of National Donut Day, which to my surprise is a real thing and not just made up for... read more >>

]]>Sean's pick this week is Forbidden Donut by Jerome Briot.

In honor of National Donut Day, which to my surprise is a real thing and not just made up for restaurant marketing, and to celebrate Homer Simpson's first live interview last month: Enjoy the forbidden donut!

Below is a quick GIF of me trying to catch the donut. For the full experience with audio and a short intro from Homer, download the file.

fdonut

Thanks Jerome for a humourous and enjoyable impossible game that is also a good simple example of a user interface that deals with mouse interaction.

Give it a try and let us know what you think here or leave a comment for Jerome.

Get
the MATLAB code

Published with MATLAB® R2016a

Jiro's pick this week is numcmp by Carlos Adrian Vargas Aguilera.

Some (many?) of you may have heard about issues with floating point math, and may have even read this reference. Even in our blogs, many of us have written about floating point arithmetic.

Here's a short example that illustrates this.

x = 2 - 1/3 - 1/3 - 1/3

x = 1.0000

But if we compare the result to its theoretical answer, we get that they aren't equal.

tf = isequal(x, 1)

tf = 0

If we display `x` in hexadecimal format to inspect the full precision,

```
format hex
x
```

x = 3ff0000000000001

On the other hand, the value "1" has a hexadecimal representation of

1

ans = 3ff0000000000000

Note the difference in the last bit.

Because of this, usually it is not a good practice to do a straight up comparison of floating point arithmetics. Instead you may do things like this.

```
format short
tf = abs(x - 1) <= eps(1)
```

tf = 1

Carlos's `numcmp` allows you to compare two values within a specific tolerance. You can choose from a set of various comparisons, such as '==', '~=', '<', '>', etc. You can also select the tolerance, specified by a positive integer `TOL` which represents `10^(-TOL)`.

tf = numcmp(x,'==',1,10) % tolerance of 1e-10

tf = 1

Thanks for the entry, Carlos. One comment. I like that it is vectorized, but you use `repmat` to match the size of the inputs. You even have a note saying that you "avoided using `bsxfun`". I actually recommend using `bsxfun`. It is much more efficient in terms of speed and memory, especially for larger data.

**Comments**

Give this a try, and let us know what you think here or leave a comment for Carlos

Get
the MATLAB code

Published with MATLAB® R2016a

Jiro's pick this week is CityName by Richard Moore.When I used to travel to various places in the U.S. for work, I would wonder which cities I was flying over during the flight. On some planes, they have seat back monitors that show the flight path with names of cities... read more >>

]]>Jiro's pick this week is CityName by Richard Moore.

When I used to travel to various places in the U.S. for work, I would wonder which cities I was flying over during the flight. On some planes, they have seat back monitors that show the flight path with names of cities the plane is flying over.

With Richard's CityName, I can do the same. I simply call the function like this

[city, country, dist] = CityName(lat, lon)

and it returns the closest city to the latitude and logitude location based on the GeoNames 5000 cities. It also returns the distance (in km) from the city center.

Lets say that I am flying from Boston, Massachusetts to Los Angeles, California.

boston = [42.36 -71.06]; los = [34.052 -118.244];

We can plot the map of the U.S. using the Mapping Toolbox.

ax = usamap('conus'); states = shaperead('usastatehi', 'UseGeoCoords', true,... 'Selector', {@(name) ~any(strcmp(name,{'Alaska','Hawaii'})), 'Name'}); geoshow(ax, states, 'DisplayType', 'polygon', 'FaceColor', [0.5 1 0.5]) framem off, gridm off, mlabel off, plabel off

Next, we can calculate the great circle track between the two cities, again using `track2` from the Mapping Toolbox. Yes, this is not the path a plane would take, but for the purpose of this post, this will do.

[lat,lon] = track2(boston(1),boston(2),los(1),los(2));

Now, we can write a loop for each coordinate to find the closest city within 10 km of the path.

% Figure out view limits for animation purpose sLim = [1386896 3288804 4332357 5404342]; eLim = [-2438739 -536831 3335946 4407931]; limits = [linspace(sLim(1),eLim(1))',linspace(sLim(2),eLim(2))',... linspace(sLim(3),eLim(3))',linspace(sLim(4),eLim(4))']; right = true; % Used for text placement for id = 1:100 % Plot a segment for the path if id < 100 plotm(lat(id:id+1),lon(id:id+1),'Color','r','LineWidth',2) end % Set axis limits axis(limits(id,:)) % Find closest city [c,~,dist] = CityName(lat(id),lon(id)); % If the distance is within 10 km, display name of the city if dist < 10 c = strtrim(c); % Place point plotm(lat(id),lon(id),'Color','b','Marker','o','MarkerFaceColor','b') % Display name of city. Alternate placing text on right and left. if right textm(lat(id),lon(id),c,... 'HorizontalAlignment','left','VerticalAlignment','top',... 'FontSize',12,'FontWeight','bold','Color','r') else textm(lat(id),lon(id),c,... 'HorizontalAlignment','right','VerticalAlignment','bottom',... 'FontSize',12,'FontWeight','bold','Color','r') end right = ~right; end drawnow end

**Comments**

Give this a try and let us know what you think here or leave a comment for Richard.

Get
the MATLAB code

Published with MATLAB® R2016a

Sean's pick this week is Flatten (Nested) Cell Arrays by Yung-Yeh. Background I recently have been using webread to scan websites and mine information to do data analysis... read more >>

]]>Sean's pick this week is Flatten (Nested) Cell Arrays by Yung-Yeh.

I recently have been using `webread` to scan websites and mine information to do data analysis with.

This requires a lot of regular expressions. Regular expressions are one of those things that are incredibly frustrating but fun at the same time. I'm usually looking
for some text pattern inside of html tags which means I'm going to be grabbing *'tokens'* or the unknown part that matches the expression.

Let's do a simple example where we grab the list of MathWorks' products from the website www.mathworks.com/products.

First, let's identify the pattern we'll be looking for. I like to do this in a web browser thanks to the syntax highlighting and other editor features. Here's a first view:

There are two patterns we need to parameterize. The first, in yellow is the product reference. The second, in green, is the product name, the token we want to capture.

Now let's code this up. First, we'll grab the html.

`html = webread('http://www.mathworks.com/products/');`

Next, we'll build the regular expression. It's always nice to keep the doc page open for this.

- Match the string literal "/product/"
- Match any "[]" characters "\w" or hyphens "\w"
- As many times as possible "*"
- Match the next backslash and closing double quote and greater than sign "/""
- Start the token with parenthesis (
- Match any words, hyphens, or spaces "\s" as many times as possible.
- Close the token )

`expr = '/products/[\w\-]*/">([\w\-\s]*)';`

Run the regular expression capturing tokens.

`tokens = regexp(html, expr, 'tokens');`

This is where Yung-Yeh's file comes in. The output from `regexp` with tokens is a nested cell that can have many nesting levels depending on number of tokens and token nesting level. `cellflat` allows me to flatten it as many levels as necessary into a cell string.

tokens = cellflat(tokens);

Now we can look at the unique products with white space at the end removed.

disp(unique(deblank(tokens))')

'...' 'Aerospace Blockset' 'Aerospace Toolbox' 'Antenna Toolbox' 'Audio System Toolbox' 'Bioinformatics Toolbox' 'Communications System Toolbox' 'Computer Vision System Toolbox' 'Control System Toolbox' 'Curve Fitting Toolbox' 'DO Qualification Kit' 'DSP System Toolbox' 'Data Acquisition Toolbox' 'Database Toolbox' 'Datafeed Toolbox' 'Econometrics Toolbox' 'Embedded Coder' 'Filter Design HDL Coder' 'Financial Instruments Toolbox' 'Financial Toolbox' 'Fixed-Point Designer' 'Fuzzy Logic Toolbox' 'Global Optimization Toolbox' 'HDL Coder' 'HDL Verifier' 'IEC Certification Kit' 'Image Acquisition Toolbox' 'Image Processing Toolbox' 'Instrument Control Toolbox' 'LTE System Toolbox' 'MATLAB' 'MATLAB Coder' 'MATLAB Compiler' 'MATLAB Compiler SDK' 'MATLAB Distributed Computing Server' 'MATLAB Mobile' 'MATLAB Production Server' 'MATLAB Report Generator' 'MATLAB for Home Use' 'Mapping Toolbox' 'Model Predictive Control Toolbox' 'Model-Based Calibration Toolbox' 'Neural Network Toolbox' 'OPC Toolbox' 'Optimization Toolbox' 'Parallel Computing Toolbox' 'Partial Differential Equation Toolbox' 'Phased Array System Toolbox' 'Polyspace Bug Finder' 'Polyspace Code Prover' 'RF Toolbox' 'Robotics System Toolbox' 'Robust Control Toolbox' 'Signal Processing Toolbox' 'SimBiology' 'SimEvents' 'SimRF' 'Simscape' 'Simscape Driveline' 'Simscape Electronics' 'Simscape Fluids' 'Simscape Multibody' 'Simscape Power Systems' 'Simulink' 'Simulink 3D Animation' 'Simulink Code Inspector' 'Simulink Coder' 'Simulink Control Design' 'Simulink Design Optimization' 'Simulink Design Verifier' 'Simulink Desktop Real-Time' 'Simulink PLC Coder' 'Simulink Real-Time' 'Simulink Report Generator' 'Simulink Test' 'Simulink Verification and Validation' 'Spreadsheet Link' 'Stateflow' 'Statistics and Machine Learning Toolbox' 'Symbolic Math Toolbox' 'System Identification Toolbox' 'Trading Toolbox' 'Vehicle Network Toolbox' 'Vision HDL Toolbox' 'WLAN System Toolbox' 'Wavelet Toolbox'

Give it a try and let us know what you think here or leave a comment for Yung-Yeh.

Get
the MATLAB code

Published with MATLAB® R2016a

Greg's pick this week is Calling Shared Libraries from Simulink by Mikhail.

Perhaps you have a DLL that you want to access from a Simulink model. Maybe a colleague has developed an algorithm in another language that can be compiled as a shared library, and you want to take advantage of it.

How can you access it from Simulink?

If you're familiar with Simulink, you're probably groaning as the words "S-Function" cross your mind, but there are other options, which can be easier.

Mikhail's entry that describes a variety of methods you can apply to access the functions and variables that are exposed in the shared library from Simulink.

This explanation of shared libraries covers the basics. Shared libraries make for a relatively portable means of sharing functionality between executables. I say relatively because shared libraries are specific to an operating system. On Windows they are DLL-files and Linux they are SO-files.

Your colleague has written some set of functions in C-code used for some other application. If your colleague compiles these functions to be exposed via a DLL, each of you can leverage the functions from your own applications, without needing to share the source code. (Although you generally need to know the interface definitions of the exposed functions)

Mikhail has a very nice document that describes several methods for accessing a shared library from Simulink.

- C-code S-function using Legacy Code Tool
- MATLAB Function block using coder.ceval
- MATLAB System block using coder.ceval
- Stateflow

As part of the examples, it also demonstrates how to load a shared library into MATLAB using the LOADLIBRARY function.

And (my favorite) demonstrates how to load a library that was generated from Simulink using Embedded Coder.

Yes, search "loadlibrary" on the MATLAB File Exchange to give you an idea of how people might use this type of capability.

To start, this entry is well formed and well documented.

Also for me, coming across this entry is very timely. I have been working on some projects that involve deploying Simulink models as MEX-files. Notice that you can generate a DLL from Simulink, and then load that DLL into MATLAB. I have used that very method, but instead of using LOADLIBRARY, I call the DLL from a MEX-file that has been compiled for MATLAB.

- Do you need to access shared libraries from your Simulink model?
- Do you want access to simulations exported from other tools?
- Should C-code generation support calling a function in a DLL?
- Is it useful to call Simulink models as a MEX-function?

Let us know here.

Get
the MATLAB code

Published with MATLAB® 9.0

Richard is a Consultant at MathWorks focused on the Embedded Coder product for code generation, primarily in the Aerospace industry. Richard's pick this week is... read more >>

]]>*Richard is a Consultant at MathWorks focused on the Embedded Coder product for code generation, primarily in the Aerospace
industry.*

Richard's pick this week is Device Drivers by Giampiero Campa.

My pick this week is **Developing Simulink Device Driver Blocks: Step-by-Step Guide and Examples** .

One of the many capabilities of Simulink is the ability to generate code for an embedded application. However, when hardware is involved, the complexity goes up exponentially. As with any embedded system, the core algorithm must interface with hardware components to be effective. These can range from simple I/O, to complex sensors such as a GPS. In addition, there are myriad communication channels that can be used such as serial, A2D/D2A, and other communication protocols (i.e. I2C). In all of these cases, there are some low-level hardware routines that provide a bridge between the hardware and the core algorithm - device drivers.

That being said, the issue then becomes how do you call these low-level routines? One way would be to write all of the code by hand. This is not very efficient and not recommended. Another would be to write the base code for the application by-hand, including the interface to the hardware, and calling the algorithmic code produced via Simulink. In reality, this is an approach adobted by many users, particularly where it is a highly complex system with many aspects, with a limited number being generated from Simulink. However, for our case, we want to be able to generate code from our Simulink model and run it directly on our hardware, especially if we are using the Run on Target Hardware feature in Simulink rather than using one of the Coder products. To do this, we need a block in Simulink that will translate into a call to the appropriate device driver.

That is where this post comes into play. This post walks you through the steps needed to create these device driver blocks. In particular, the post provides information on the various approaches that can be used such as S-Function blocks and MATLAB Function blocks, and the pros/cons of each approach. The post provides examples for creating blocks for an Arduino board, but the approach is applicable to any board. This is supported by the number of additional posts that were inspired by this post for a wide range of boards and applications.

Once a device driver block is created, the user can place it in a Simulink library for reuse. In this way, various device driver blocks can be packaged and shared for use in different applications. Many of the Target Support Packages leverage this capability.

If you need to interface with hardware, this post can provide the necessary guidance to allow you to create a set of hardware specific blocks. Give it a try and let us know what you think here or leave a comment for Giampiero.

Get
the MATLAB code

Published with MATLAB® R2016a