Building a Microscopy Application in Mathematica

September 9, 2011

Jon McLoone, International Business & Strategic Development

As a change from my usual recreational content, today I thought I would describe a real Mathematica application that I wrote. The project came from my most important Mathematica user—not because she spends a lot of money with Wolfram Research, but because I am married to her!

Her company, Particle Therapeutics, works on needle-free injection devices that fire powdered drug particles into the skin on a supersonic gas shock wave. She was trying to analyze the penetration characteristics on a test medium by photographing thin slices of a target under a microscope and measuring the locations of the particles.

The problem was that her expensive image processing software was doing a poor job of identifying overlapping particles and gave her no manual override for its mistakes.

Faced with the alternative of holding rulers up to her screen and recording each value by hand, I promised that I could do better in Mathematica, with the added advantage that now her image processing tool would be integrated into her analysis code to go from image file to report document in a single workflow.

Here is a detail from a typical image:

image=ImageTake[Import["022108shot06.jpg"],All,{1,800}]

The task breaks into three parts:

1. Image processing
2. An optimal GUI for adjusting parameters
3. Report generation

Step 1 seems like it should be the hardest, but it can be achieved in very little code. First I need to clean up the image by deleting dirt and texture. DeleteSmallComponents removes small white dots within the particles; inverting and doing it again with a parameter lets me delete large black components up to the parameter size.

image2=DeleteSmallComponents[ColorNegate[DeleteSmallComponents[Binarize[image]]], 200]

Once cleaned, I need to erode the components back to a single pixel. But I need multiple particles that slightly overlap to break apart and erode back to separate points. This “ultimate erosion” can be done by finding the local maxima of the DistanceTransform of the image.

ImageAdjust[DistanceTransform[]]

MaxDetect[DistanceTransform[]]

I can then use ComponentMeasurements to get information on the identified points—in this case the coordinates of the center.

ComponentMeasurements[MaxDetect[DistanceTransform[image2]],"Centroid"]

As a final cleanup, I use DeleteDuplicates to remove spurious points that come from badly shaped particles. I give a user parameter for how close points should be before they are considered the same. Here is all of that packaged into a function:

identifyBlobs[image_, minBlob_, minSeperation_] := DeleteDuplicates[Last /@ ComponentMeasurements[MaxDetect[DistanceTransform[DeleteSmallComponents[ColorNegate[DeleteSmallComponents[Binarize[image]]], minBlob]]], "Centroid"], EuclideanDistance[#1, #2] < minSeperation &]

It turns out that, already, this does better than her image processing software, which was only capable of standard Erosion that treats all touching particles as one.

The hardest part is the GUI, because I want to give lots of control. Here is a video of the final application in action:



Let’s walk through each of those elements. 

First I want this to be in a palette of its own, not in a notebook. But I want multiple instances to be able to be open at once. So my code starts by localizing variables into the window instance (not globally to the Mathematica session).

CreatePalette[DynamicModule[{img=Image[{{1}}],file="",data},...]]

The main guts of the user interface (UI) are then just a Manipulate that first calls the identifyBlobs command, and then displays the results.

Manipulate[
data=identifyBlobs[img,minSize,minSeperation];
{w,h}=ImageDimensions[img];
Show[img,Graphics[{Thickness[.001],Cyan,Line[{{0,h l},{w,h r}}]}]],...]

identifyBlobs will be called automatically when any of the parameters change, and the image will update automatically without any extra work.

The clever bit is the use of locators, which means that the “+” signs for the identified positions are not just for display, but can be changed by dragging them. With the LocatorAutoCreate option, I can Alt-click to add or remove points and automatically update the particle location data variable. That’s a lot of UI for one line of code.

{data, Locator, LocatorAutoCreate -> True, Appearance -> Style["+", 20, White]},...

The sliders for the left and right side of the surface line are just Manipulate controls, except that I choose vertical sliders rather than horizontal and place them in specific order. (It is possible to auto-detect the surface with ImageLines, but in some images the surface is so faint that this manual control is more practical.)

{{l, 0.9, ""}, 0, 1, ControlPlacement -> Left, ControlType -> VerticalSlider},
{{r, 0.9, ""}, 0, 1, ControlPlacement -> Right, ControlType -> VerticalSlider},...

Then we have a button to load the data. It’s mostly error trapping code, but note the 
Method -> “Queued” that prevents the Button from timing out if I spend time dithering about which file to open.

Button["Import image", If[MatchQ[#, _Image], img = #] &[Quiet[Import[file = SystemDialogInput["FileOpen"]]]], Method -> "Queued"]

I hide the advanced controls using OpenerView, and I specify the layout of the controls inside that.

OpenerView[{"auto detection settings", Column[{
    Control[{{minSize, 50, "minumum particle size"}, 0, 200}],
    Control[{{minSeperation, 7, "minumum seperation"}, 0, 20}]
    }]}]

Finally, to make the scrollable area, all I have to do is set the ContentSize to manual, with an initial size:

ContentSize -> {700, 400},...

Now on to reporting: I can't share the details of the analysis she was doing, but here is a rough outline for how reporting works. Mathematica provides a complete document description language, so all I have to do is describe the layout and style of the document in terms of the individual visualizations and analysis that I have done.

CreateDocument[{
  TextCell["Particle impact analysis", "Title"],
TextCell[Row[{"Sample: ", path, "\n", "Report generated ", DateString[], "\n", "Particles found: ", Length[data]}], "Text"],
  ExpressionCell[DistributionChart[penetrationData[[All, 1]], PlotLabel -> "Spread profile", BarOrigin -> Left], "Output"],...}]

This creates a new document that starts with three cells: first a static title, then some text summarizing results, and then a chart.

And that’s it. 

In a sensible width window, UI, image processing, analysis, and report generation, including a few extra charts, amount to about 30 lines of code and can be deployed to Mathematica or Mathematica Player Pro. Here is the entire application put together:

CreatePalette[DynamicModule[{img=Image[{{1}}],file="",data},

Manipulate[
data=identifyBlobs[img,minSize,minSeperation];
{w,h}=ImageDimensions[img];
Show[img,Graphics[{Thickness[.001],Cyan,Line[{{0,h l},{w,h r}}]}]],

(*Interface*)
{data,Locator,LocatorAutoCreate->True,Appearance->Style["+",20,White]},

Button["Import image",If[MatchQ[#,_Image],img=#]&[Quiet[Import[file=SystemDialogInput["FileOpen"]]]],Method->"Queued"],

Button["generate report",GenerateReport[file,data,h l,h r,w]],

OpenerView[{"auto detection settings",Column[{
Control[{{minSize,50,"minumum particle size"},0,200}],
Control[{{minSeperation,7,"minumum seperation"},0,20}]
}]}],

{{l,0.9,""},0,1,ControlPlacement->Left,ControlType->VerticalSlider},
{{r,0.9,""},0,1,ControlPlacement->Right,ControlType->VerticalSlider},

ContentSize->{700,400},

Initialization:>(
(*Image recognition function*)
identifyBlobs[image_, minBlob_, minSeperation_] := DeleteDuplicates[Last /@ ComponentMeasurements[MaxDetect[DistanceTransform[DeleteSmallComponents[ColorNegate[DeleteSmallComponents[Binarize[image]]], minBlob]]], "Centroid"], EuclideanDistance[#1, #2] < minSeperation &];

depth[{x_,y_},w_,l_,r_]:={x,l+(x/w)(r-l)-y};

(*Report generation function*)
GenerateReport[path_,data_,l1_,l2_,w_]:=Block[{penetrationData},
penetrationData=depth[#,w,l1,l2]&/@data;

If[Length[data]>0,CreateDocument[{
TextCell["Particle impact analysis","Title"],

TextCell[Row[{"Sample: ",path,"\n","Report generated ",DateString[],"\n","Particles found: ",Length[data]}],"Text"],

ExpressionCell[DistributionChart[penetrationData[[All,1]],PlotLabel->"Spread profile",BarOrigin->Left],"Output"],

ExpressionCell[DistributionChart[penetrationData[[All,2]],PlotLabel->"Depth profile",BarOrigin->Top],"Output"],

ExpressionCell[ListLinePlot[Map[Mean[#[[1]]]&,DeleteCases[BinLists[penetrationData,w/25,10^10],{},∞]],PlotLabel->"Mean depth"],"Output"],

ExpressionCell[DensityHistogram[penetrationData,8,PlotLabel->"Density"],"Output"],

TextCell["Particle coordinates","Section"],

ExpressionCell[Sort[Round[penetrationData]],"Input"]}]]];
)]],WindowTitle->"Particle analyzer"];