Showing posts with label MATLAB. Show all posts
Showing posts with label MATLAB. Show all posts

Monday, September 5, 2016

The Pentium Bug

I remember the Pentium Bug in 1994 (I had just started my PhD) but didn't know enough about scientific computing to understand the significance.

Here's a fascinating blog post from Cleve Moler describing what happened.

Saturday, September 3, 2016

POVRay sarcomeres - final touches

This post tidies up some loose ends about making movies of muscle sarcomeres in POV-Ray. I blogged earlier about

Putting all these techniques together allowed me to render an image of chain of sarcomeres.

The most difficult part here was defining the geometry and making sure that I was drawing each object in the correct x,y,z location. I automated this process using MATLAB. The whole scene is made up of spheres, cylinders, and boxes. A few objects are merged together as blobs.

Once I had that working for a single image, I just moved the camera around the scene until I had a sequence that I thought looked interesting. Then I rendered the images (it took a few hours on my PC) and stitched the sequence together to make a movie.

Here's the final result.



I am happy to share the code but it may take me some time to get it posted online. If you want to get it sooner, just send me a note via a comment below.

Thursday, August 25, 2016

How to learn to code?

Coding (perhaps like most things) is all about perspective.

My academic appointment is the Department of Physiology in the Medical School at the University of Kentucky. Only a few people in the department know how to program and, as a result, most of my colleagues think I'm 'pretty good with computers'.

My collaborators in Engineering and Chemistry across the street are always polite but think I code like a Geico Caveman.

Why do so few physiologists code? I used to think it's because biomedical scientists didn't think it was useful but I've changed my mind. I now think that our students don't code because they've never been taught. They've typically taken lots of courses in biology, chemistry, and genetics, but nobody has shown them how to write code, or how to break a task down into little steps so that they can create an algorithm.

What's the best way of helping these students get started?

I've blogged about how I  learned to program before but my advice won't work for everybody and I'm always looking for new resources that I can tell students about.

Today, I came across the MATLAB OnRamp. This resource has actually been around for several years but I hadn't seen it before and it looks interesting. If you have any experience with it, or can suggest better ways of getting started in scientific computing, please share your thoughts in the comments below. 


Monday, August 22, 2016

Filtering binary images by object size in MATLAB

As usual, I learned something new when I read Steve Eddin's latest blog post.

This time around, I discovered that MATLAB has some new features that make it easy to filter binary objects based on their size. There were, of course, old ways to do this but the new bwareafilt command (at least, new in MATLAB 2014B) simplifies some workflows.


Sunday, August 21, 2016

POV-Ray animation

This post continues my introduction to generating images and movies using POV-Ray.

In the last post, I showed how to render a simple scene. Now I am going to show how to animate it.

POV-Ray includes a lot of options for animation and this page provides some great tutorials and examples. However, I have always found it easier to generate animations by automatically generating different *.pov files and then rendering them as part of my workflow.

The last part is easy because you can call POV-Ray directly from the command line. (See here for many options.)

For example, once you are in the directory that contains the main POV-Ray executable, you can produce a PNG file that is 340 by 280 pixels from the scene described in c:\temp\out.pov using this command.

pvengine +W340 +H280 c:\temp\out.pov /exit

Here's some MATLAB code that will spin the camera around the previous image.

function povray_animation

r = -25;
no_of_points = 50;
x = r*cosd(linspace(0,360,no_of_points))

for i=1:no_of_points
    
    pov_file = fopen(fullfile(cd,sprintf('out%.0f.pov',i)),'w');
    fprintf(pov_file,'#include "colors.inc"\n');
    fprintf(pov_file,'background {color White}\n');
    fprintf(pov_file,'camera {location <%f,25,25>\n',x(i));
    fprintf(pov_file,'\tsky <0,0,1>\n');
    fprintf(pov_file,'\tlook_at <10,0,0>}\n');
    fprintf(pov_file,'light_source {<%f,25,25> color White}\n',x(i));
    fprintf(pov_file, ...
        'cylinder {<0,0,0>,<10,0,0>,1 texture {pigment {color Red}}}\n');
    fprintf(pov_file, ...
        'cylinder {<0,0,0>,<0,10,0>,1 texture {pigment {color Green}}}\n');
    fprintf(pov_file, ...
        'cylinder {<0,0,0>,<0,0,10>,1 texture {pigment {color Blue}}}\n');
    fprintf(pov_file, ...
        'sphere {<10,10,10>,3 texture {pigment {color OrangeRed}}}\n');
    fclose(pov_file);
    
    cd_string = 'cd c:\program files\pov-ray\v3.7\bin';
    command_string = sprintf('%s\npvengine +W340 +H280 %s +A /exit', ...
        cd_string, ...
        fullfile(cd,sprintf('out%.0f.pov',i)));
    batch_file = fopen('pov.bat','w');
    fprintf(batch_file,'%s\n',command_string);
    fclose(batch_file);
    system('pov.bat');
end

It's now easy to stitch the *.png files together to make a movie.




Wednesday, August 10, 2016

Word clouds

A quick aside - I will get back to POV-Ray soon.

I wanted to make a Word Cloud for something at work. I've never done this before so I adopted my usual strategy - look in the MATLAB FileExchange to see if somebody has already created a tool.

As is nearly always the case, I found something useful, WordData Visualization. It was even a Pick of the Week back in 2015.

Five minutes later, and I had what I wanted. Then I went on to generate a figure based on the text from my PhD thesis. It doesn't get much more exciting :-)





Sunday, August 7, 2016

POV-Ray - old but gold

I wanted to render some 3D images recently. I downloaded SketchUp and Blender and did some quick googling.

Both packages look really interesting but they look as if they are primarily GUI-based. I realize that's what most people probably want but I needed to render a scene that contains thousands of objects. I didn't want to have to add these into my scene one at a time.

I'm sure somebody will explain how to create complex scenes in SketchUp and Blender (please write a comment if you do) but I decided to go back to POV-Ray. It's old (the last release was in 2013) but it's very functional. I used MATLAB to create my scene in the POV-Ray format (*.pov) and then rendered the image using the POV-Ray engine.

I'll explain how I automated the process in the next few posts and created 'fly-through movies' with a single click but, for now, here's a quick example of the sort of image I was able to generate.

I think it's pretty cool.


Monday, May 23, 2016

Numerical integration

I love reading everything Cleve Moler writes about scientific computing. He makes things sound so simple and I always learn something new. In fact, I nearly always learn lots of things.

Here's his latest blog post on numerical integration.

http://blogs.mathworks.com/cleve/2016/05/23/modernization-of-numerical-integration-from-quad-to-integral/

Monday, May 9, 2016

New resources for learning about differential equations

Students often ask me how they can start to learn about scientific computing.

I normally tell them to do what I did - start with a problem that they are interested in (maybe one that they are working on?) and begin to tinker. Here's how that worked out ...



But I realize that there are probably better approaches. For example, many people now like to learn by watching videos and tutorials.

This blog post caught my eye this morning and the videos, lectures, and online material look super interesting. I hope to give these courses a try myself.

http://blogs.mathworks.com/cleve/2016/05/09/strang-and-moler-video-course-on-differential-equations/

Monday, May 2, 2016

Neural networks for beginners

I've been thinking a little about machine learning recently so this blog post caught my attention.

http://blogs.mathworks.com/loren/2015/08/04/artificial-neural-networks-for-beginners/?s_eid=PSM_da

It's a beautiful post that starts with a concrete example focused on using artificial networks to recognize digits and ends with a program that solves Soduko puzzles.

Thursday, April 28, 2016

New options for thresholding images in MATLAB

I've written several posts recently about thresholding images. I was thus pretty excited to see Steve Eddins announce a new series of posts on his blog.

http://blogs.mathworks.com/steve/2016/04/28/image-binarization-new-functional-designs/

My impression from a very quick scan of the documentation is that MATLAB 2016a includes a new function imbinarize which can be called with a locally adaptive threshold option.

I don't think this will lead to a breakthrough in image processing - adaptive thresholding is not a new technique - but it will simplify coding. The new function will allow users to replace multiple lines of 'old' MATLAB code with a single command.

And in my opinion, that's why MATLAB is so useful for scientific computing. It's not that the software gives you completely new options; you can code anything in C if you want to. It's just that MATLAB makes it easy to try things without investing huge amounts of time. That's why it accelerates science.

Saturday, April 16, 2016

Image segmentation - part 3

This is the third post in our discussion of image segmentation. The first post talked about splitting an image by RGB color using ImageJ / Fiji. Then I talked about using k-means segmentation and the L*a*b color space. You might remember that I was pretty excited about the technique. Is it a general solution to all segmentation problems?

Let's try it on this image.


When I look at this picture, I see 3 colors: white, black/blue and lots of red. Let's convert from RGB to L*a*b and then split into 3 regions using k-means segmentation. The MATLAB code to do this is

function k_means_clustering

image_file_string = '3mokotri.tif';

% Load image
im_rgb = imread(image_file_string);

% Transform to lab color space
cform = makecform('srgb2lab');
im_lab = applycform(im_rgb,cform);

% Create a n x 2 matrix of [a b] values
a = im_lab(:,:,2);
a = double(a(:));
b = im_lab(:,:,3);
b = double(b(:));
ab = [a b];

% Construct the 2D histogram with appropriate limits
c{1}=[min(ab(:)) : max(ab(:))];
c{2}=c{1};
n = hist3(ab,c);
n1 = n';
n1(size(n,1) + 1, size(n,2) + 1) = 0;

% Display
figure(1);
clf;
imagesc(log10(n1));
hold on;
title('2D histogram of pixel values, colored based on log10 of pixel count')
xlabel('a value');
ylabel('b value');

[id,cl]=kmeans(ab,3);

figure(2);
clf;
hold on;
cm = paruly(3);
for i=1:3
    vi = find(id==i);
    plot(ab(vi,1),ab(vi,2),'+','Color',cm(i,:));
    plot(cl(i,1),cl(i,2),'md','MarkerFaceColor','m');
end
xlabel('a value');
ylabel('b value');
set(gca,'YDir','reverse');

[r,c] = size(im_rgb(:,:,1));
im_cluster = reshape(id,r,c);

figure(3);
imagesc(im_cluster);

The end result is
Hmm, that was not what I was hoping for.

Compare with the original again.

It looks like the algorithm did a good job of finding the white lines but it didn't find the black dots I was looking for. Instead, it separated the red areas in the original picture into dark red and light red (maybe pink?) areas.

I have some ideas about how to improve the algorithm for this particular image (which I intend to address in future posts) but I think I've just demonstrated a general result of image processing. It's hard to find a single approach that works for every type of task. It's also hard to develop techniques that don't have a few fiddle factors (parameters that you adjust until 'something works').

Image processing is never easy but as you gain experience and learn more about the techniques that are available, you begin to make faster progress.

Next up, some thoughts on t-tests.

Sunday, April 10, 2016

Image segmentation - part 2

This continues the previous post where we talked about segmenting a test image into 3 parts. As a reminder, here's the test image again.



k-means clustering takes your data points and tries to group them into k clusters. It's a general technique that can be used with many sorts of data but it can be useful for segmenting images.

If you want to segment by color, it's often easier to switch to the L*a*b color-space. When you do this, the second and third planes contain all of the color data. The first L plane is to do with brightness.

Here's some MATLAB code that
  1. loads our test image
  2. transforms it into the lab color space
  3. calculates the 2D histogram of pixel intensities
  4. displays the heat map of the histogram with a log color scale

function k_means_clustering

image_file_string = 'hubert_8_3_small.png';

% Load image
im_rgb = imread(image_file_string);

% Transform to lab color space
cform = makecform('srgb2lab');
im_lab = applycform(im_rgb,cform);

% Create a n x 2 matrix of [a b] values
a = im_lab(:,:,2);
a = double(a(:));
b = im_lab(:,:,3);
b = double(b(:));
ab = [a b];

% Construct the 2D histogram with appropriate limits
c{1}=[min(ab(:)) : max(ab(:))];
c{2}=c{1};
n = hist3(ab,c);
n1 = n';
n1(size(n,1) + 1, size(n,2) + 1) = 0;

% Display
figure(1);
clf;
imagesc(log10(n1));
hold on;
title('2D histogram of pixel values, colored based on log10 of pixel count')
xlabel('a value');
ylabel('b value');


The output is

When I look at the test image, I see three colors: a gray/white that forms the background, an orangey-yellow that is muscle tissue, and a darker red that indicates collagen. So let's try to split the pixel values into 3 separate groups. In MATLAB, the command is

[id,cl]=kmeans(ab,3);

id is a vector with values of 1,2, or 3. These link each pixel to a given cluster.
cl is a 3x2 array which shows the centroids of the clusters.

This snippet of code shows which [a,b] values fall into each cluster. The magenta markers show the centroids (cl) of each cluster.

figure(2);
clf;
hold on;
cm = paruly(3);
for i=1:3
    vi = find(id==i);
    plot(ab(vi,1),ab(vi,2),'+','Color',cm(i,:));
    plot(cl(i,1),cl(i,2),'md','MarkerFaceColor','m');
end

set(gca,'YDir','reverse');

giving
It's now easy to set each pixel in the original image based on it's cluster identity. You just have to reshape the id vector.

[r,c] = size(im_rgb(:,:,1));
im_cluster = reshape(id,r,c);

figure(3);
imagesc(im_cluster);

The output is
Compare this result to the original.




The segmentation is not perfect but we didn't set any parameters other than 'we want 3 colors'. It's also super-fast and objective. I think it's pretty cool.

Is k means clustering a perfect solution for image segmentation? We will look at that in the next post.

Thursday, March 31, 2016

Linking MATLAB to REDCap

REDCap is a mature secure web application for building and managing online surveys and databases. It's used extensively in the biomedical community to store and manage clinical data. The Campbell lab uses it to manage clinical data about the patients and organ donors who donate cardiac samples for our research.

Although REDCap has many advantages, I've found exporting data to be tricky. In principle, you can export data in formats for

  • SPSS
  • SAS
  • R
  • Stata
  • CSV / Excel

In practice, everything but the CSV / Excel format requires downloading two or three files. One of these files has the data, the other files have code that tells the statistics packages how to import the information. I haven't experimented with this approach because I'm not keen on having to link my data files to specific import code. That seems complicated.

As noted above, REDCap labels its CSV output format as "CSV / Excel". It's important to note though that the file is not an xlsx spreadsheet. It's a CSV file, which uses commas to separate entries. It might look like this

Date,Name,Age,Comments,BMI
11/22/2015,Mike,22,Nothing,24.8
11/23/2015,Ken,23,Interestingly, Ken had red hair, and also blue eyes,25.2

This looks okay at first, but you quickly get into trouble if any of the fields include commas. For example, the second line would show a BMI of Ken, instead of 25.2.

I ran into this problem with one of our datasets and developed several workarounds. They did the job but I wasn't confident that I would catch every new error going forward. Then I remembered that REDCap had an API and I thought I would give that a look.

This post gave me some useful pointers but I don't have much experience with Python so I had trouble getting things to work. The code in the REDCap API Sandbox also seemed a bit buggy.

The real breakthrough came when I realized that  cURL makes it easy to communicate via http. I used these examples and some snippets from the REDCap API and was quickly able to download a REDCap report in JSON format.

For example, the command generated by this print statement


sprintf('token=%s&content=report&format=json&report_id=2313&rawOrLabel=label&rawOrLabelHeaders=label&exportCheckboxLabel=true&returnFormat=json" https://redcap.uky.edu/redcap/api/ -o %s',my_token,output_file_name);

sucks my REDCap report from the server and saves it to the specified text file.

I then used the MATLAB JSONlab toolbox to turn the data into a MATLAB structure which is what I really wanted for my data processing anyway.

Bottom line, from now on, I think I'm going to give up exporting data from REDCap directly, and suck it into MATLAB using the REDCap API. That way, I get around some of the limitations of the CSV (or REDCap alternative) outputs and I can get the latest version of the database on the fly whenever I need it.