Jiro's pick this week is smoothn by Damien Garcia.As the title of Damien's entry states, smoothn is a fast and easy smoothing function for n-dimensional data. It uses a generalized cross-validation method to estimate the smoothing parameter, which affects the quality of the output. In layman's terms, it automatically does... read more >>

]]>Jiro's pick this week is `smoothn` by Damien Garcia.

As the title of Damien's entry states, `smoothn` is a fast and easy smoothing function for n-dimensional data. It uses a generalized cross-validation method to estimate the smoothing parameter, which affects the quality of the output. In layman's terms, it automatically does things for you so that you just pass in the noisy data and out comes the smoothed data. Damian includes links to two of his papers for technical reference for those interested in the detailed theory.

From what I've seen, his code is very well written (Code Analyzer is green!), with plenty of error-checking and good use of vectorized code. It has extensive help, and he includes many different types of examples. Here is one that I liked.

**Noisy image with missing data**

Filling in missing data in an image?! What?!

n = 256; % Original image y0 = peaks(n); % Noisy image with missing data y = y0 + randn(size(y0))*2; I = randperm(n^2); y(I(1:n^2*0.5)) = NaN; % lose 1/2 of data y(40:90,140:190) = NaN; % create a hole % Smoothed image z = smoothn(y); % <-- That's all!! % View subplot(2,2,1:2) imshow(y,[]) title('Noisy corrupt data') subplot(2,2,3) imshow(z,[]) title('Recovered data ...') subplot(2,2,4) imshow(y0,[]) title('... compared with the actual data') colormap(parula)

**Comments**

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

Get
the MATLAB code

Published with MATLAB® R2016a

Contents Hardware Support for Image Aquisition from Kinect v2 Using the Kinect for Windows v2 in MATLAB Initializing the device Get a 3-D point cloud from the device Viewing the point cloud stream from the Kinect v2 Detect planes in the 3-D point cloud Release the device Avi's pick of the... read more >>

]]>colorDevice = imaq.VideoDevice('kinect',1) depthDevice = imaq.VideoDevice('kinect',2)

colorDevice = imaq.VideoDevice with properties: Device: 'Kinect V2 Color Sensor (kinect-1)' VideoFormat: 'BGR_1920x1080' ROI: [1 1 1920 1080] ReturnedColorSpace: 'rgb' ReturnedDataType: 'uint8' DeviceProperties: [1x1 imaq.internal.DeviceProperties] depthDevice = imaq.VideoDevice with properties: Device: 'Kinect V2 Depth Sensor (kinect-2)' VideoFormat: 'Depth_512x424' ROI: [1 1 512 424] ReturnedColorSpace: 'grayscale' ReturnedDataType: 'uint16' DeviceProperties: [1x1 imaq.internal.DeviceProperties]

step(colorDevice); % Initialize color/RGB sensor step(depthDevice); % Initialize depth sensor colorImage = step(colorDevice); % Load one color/RGB frame depthImage = step(depthDevice); % Load one depth frame

ptCloud = pcfromkinect(depthDevice,depthImage,colorImage);

player = pcplayer(ptCloud.XLimits,ptCloud.YLimits,ptCloud.ZLimits,... 'VerticalAxis','y','VerticalAxisDir','down'); xlabel(player.Axes,'X (m)'); ylabel(player.Axes,'Y (m)'); zlabel(player.Axes,'Z (m)'); % Then stream a live 3-D point cloud from the Kinect v2. for i = 1:500 colorImage = step(colorDevice); depthImage = step(depthDevice); ptCloud = pcfromkinect(depthDevice,depthImage,colorImage); view(player,ptCloud); end

colorImage = step(colorDevice); depthImage = step(depthDevice); ptCloud = pcfromkinect(depthDevice,depthImage,colorImage); % Grab a new point cloud maxDistance = 0.02; referenceVector = [0,0,1]; maxAngularDistance = 5; [model,inlierIndices,outlierIndices] = pcfitplane(ptCloud,maxDistance,referenceVector,maxAngularDistance); figure; pcshow(ptCloud); hold on; plot(model)

release(colorDevice); release(depthDevice);]]>

Sean's pick this week is Treemap by Joe Hicklin. Contents Single Treemap ... read more >>

]]>Sean's pick this week is Treemap by Joe Hicklin.

I've recently been working with one of our marketing teams to visualize where MATLAB is used across an organization or university
based on department. One nice way to do this is with a treemap because it allows for hierarchical visualization. Let's take
a look at a single layer treemap for a synthetic university. The data is stored in table called *activations* with variables *College* and *Department*.

collcounts = countcats(activations.College); % How many in each college? r = treemap(collcounts); % Build tree rectangles plotRectangles(r,categories(activations.College)) % Plot tree map

That looks pretty good, but it's just one layer. We also have information on the departments within a college. To include
this, I will use the `findgroups/splitapply` workflow to split the departments based on their colleges.

% Which college collegeidx = findgroups(activations.College); % Count the number of activations in each department and bin by college deptcounts = splitapply(@(x){nonzeros(countcats(x))},activations.Department,collegeidx); % Keep the names for labeling deptnames = splitapply(@(x){categories(removecats(x))},activations.Department,collegeidx);

Now we'll copy and modify Joe's example code to plug in my new data.

m = 39; % Department colors for ii = 1:numel(deptcounts) colors = (3*repmat(rand(1,3),m,1) + rand(m,3))/4; rNew = treemap(deptcounts{ii}',r(3,ii),r(4,ii)); rNew(1,:) = rNew(1,:) + r(1,ii); rNew(2,:) = rNew(2,:) + r(2,ii); plotRectangles(rNew,deptnames{ii},colors) end % Outline college boxes outline(r)

On my large monitor blown up that looks good. But it's kind of busy compressed. I'll enhance it to make it so that you can expand one school or "Drill down". First, repeat all of the above steps but leave the labels off the departments.

figure plotRectangles(r,categories(activations.College)); for ii = 1:numel(deptcounts) colors = (3*repmat(rand(1,3),m,1) + rand(m,3))/4; rNew = treemap(deptcounts{ii}',r(3,ii),r(4,ii)); rNew(1,:) = rNew(1,:) + r(1,ii); rNew(2,:) = rNew(2,:) + r(2,ii); plotRectangles(rNew,[],colors) end outline(r)

Second, I'll disable the ability for the patches to interact by turning off their *'HitTest'* property.

patches = findall(gca,'Type','patch'); set(patches,'HitTest','off')

Finally, add a button down callback to the axes which calls my custom explode function.

function explodeBlock(evt,r,nsubs,sublabels) % Blow up one block in a tree map % Where was hit xyhit = evt.IntersectionPoint(1:2); % Which rectangle? inx = xyhit(1) > r(1,:) & (xyhit(1) < r(1,:)+r(3,:)); iny = xyhit(2) > r(2,:) & (xyhit(2) < r(2,:)+r(4,:)); idx = inx & iny; % Blow up that rectangle in new figure rnew = treemap(nsubs{idx}); figure plotRectangles(rnew,sublabels{idx}) end

set(gca,'Visible','on') set(gca,'ButtonDownFcn',@(~,evt)explodeBlock(evt,r,deptcounts,deptnames))

Using the treemap utility is pretty straight-forward and Joe has provided nice examples of all of the features. I do have a couple of enhancements though:

`plotRectangles`and`outline`could return handles to the patches they create. This would allow me to skip the sloppy call to`findobj`that I have above.- The text, as you probably noticed, tends to overlap a bit. It would be nice for the text to be jittered or rotated so as to avoid the collision.

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

Get
the MATLAB code

Published with MATLAB® R2016b

Sean's pick this week is Alternative Box Plot by Christopher Hummersone. Background I was recently doing a seminar for some MATLAB users where one of the plots I... read more >>

]]>Sean's pick this week is Alternative Box Plot by Christopher Hummersone.

I was recently doing a seminar for some MATLAB users where one of the plots I showed was a `boxplot`. The box plot was used to quickly get a feel for the distributions of some energy data based on time of day. Here it is:

figure boxplot(Power, Hour, 'Notch', 'on') xlabel('Hour of Day') ylabel('Power (MW)')

There is a lot of information in this plot. The whiskers are the limits excluding outliers which are the red +s. The boxes represent the 25th and 75th percentiles, the red line the median, and the notches tell you if two medians are significantly different if they don't overlap.

An audience member then asked me if there was a way to control the box limits to be an arbitrary percentile, they wanted 10% and 90%. I searched the parameter-value pairs but did not find anything.

Alternative box plot provides many of the same options as MATLAB's boxplot but also a few more and the percentile is one of them.

Christopher has his tools in a package, which is nice to avoid naming conflicts with other functions. I'll `import` the package to use the short name:

import iosr.statistics.* y = tab2box(Hour,Power); % reshapes to be boxPlot compliant bp = boxPlot(0:23,y,'Notch',true,'Percentile',[10 90]);

Very nice!

You may have also noticed that I grabbed an output argument. This allows me to customize the box plot with easy tab completion after it's already made.

bp.lineColor = 'b'; bp.medianColor = 'r';

In fact, I could even have all four percentiles by adding additional ones.

bp = boxPlot(0:23,y,'Notch',true); bp.addPrctiles = [10 90]; bp.addPrctilesColors = {'r';'b'}

bp = boxPlot with properties: addPrctiles: [10 90] addPrctilesColors: {2×1 cell} addPrctilesLabels: {2×1 cell} addPrctilesMarkers: {2×1 cell} addPrctilesTxtSize: [] boxAlpha: 1 boxColor: {'none'} boxWidth: 'auto' groupLabelFontSize: 9 groupLabelHeight: 'auto' groupLabels: {1×0 cell} groupWidth: 0.7500 limit: '1.5IQR' lineColor: {'k'} lineStyle: {'-'} lineWidth: 1 meanColor: {[0 0.4470 0.7410]} meanMarker: {'+'} meanSize: 6 medianColor: {[0 0.4470 0.7410]} method: 'R-5' notch: 1 notchDepth: 0.4000 notchLineColor: {'k'} notchLineStyle: {':'} notchLine: 0 outlierSize: 36 percentile: [25 75] sampleFontSize: 9 sampleSize: 0 scaleWidth: 0 scatterAlpha: 1 scatterColor: {[0.5000 0.5000 0.5000]} scatterLayer: 'top' scatterMarker: {'x'} scatterSize: 36 showLegend: 0 showMean: 0 showOutliers: 1 showScatter: 0 style: 'normal' symbolColor: {[0 0.4470 0.7410]} symbolMarker: {'o'} theme: 'default' xSeparator: 0 xSpacing: 'x' handles: [1×1 struct] x: [1×24 double] y: [334×24 double] weights: [334×24 double] statistics: [1×1 struct]

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

Get
the MATLAB code

Published with MATLAB® R2016b

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`/skel, 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