Image Processing in Microscopy

Markus van Almsick, Wolfram Research Inc.

Join the Conversation #WolframTechConf

Abstract

Light microscopy is an essential tool in the life and material sciences. To obtain, enhance, and utilize the fascinating images of the world at μm-scale an entire workflow of processing steps is required. In this talk we introduce the tools in the Wolfram Language that contribute to this workflow ranging from brightness equilibration all the way to 3D reconstruction.

Microscopy Work-Flow

Image Import

BIO Formats

Bio-Formats is a community driven software project for reading and writing image data using standardized, open formats.

Bio-Formats 5 converts over 140 file formats to the OME-TIFF data standard. It supports High Content Screening, time lapse imaging, digital pathology, and multidimensional image formats.

BioFormats Installation Steps

Assuming that you have installed Java on your computer.

Initialize Wolframs JLink`.

Needs["JLink`"];

Reinstall Java with a large memory partition.

JLink`ReinstallJava[JVMArguments->"-Xmx4048m"];

Establish a BioFormats directory in $UserAddOnsDirectory:

BioFormatsJarStorage=FileNameJoin[{$UserAddOnsDirectory,"ImageProcessing","Java"}]

If[Not@FileExistsQ@BioFormatsJarStorage, CreateDirectory[BioFormatsJarStorage]]

Download the BioFormats jar version 5.15:

bioFormatsJar=URLSave["http://downloads.openmicroscopy.org/bio-formats/5.1.5/artifacts/bioformats_package.jar",BioFormatsJarStorage<>"bioformats_package.jar"]

Activate required Java classes.

JLink`AddToClassPath[bioFormatsJar];

JLink`LoadJavaClass["java.awt.image.BufferedImage"];

JLink`LoadJavaClass["loci.formats.gui.BufferedImageReader"];

Implementation of BioFormats Import Functions

Utility function to convert pixel types.

BioFormatToMathematicaTypeAndColorSpace[type_Integer]:=
Switch[type,
0, {"Real", Automatic},
10,{"Byte", "Grayscale"},
11,{"Bit16","Grayscale"},
12,{"Byte", "Grayscale"},
_,{"Byte", "RGB"}
];

ImportBioFormatImage function to import image data.

Clear[ImportBioFormatImage];
ImportBioFormatImage[file_String]:=
JavaBlock@Module[
{ir=JLink`JavaNew[JLink`LoadJavaClass["loci.formats.ImageReader"]],
reader,image,raster,series,channels,timeSeries,slices, index,
height, width, type,stack, colorspace = Automatic},
ir@setId[file];
reader=BufferedImageReader`makeBufferedImageReader[ir];
series=reader@getSeriesCount[];
channels=reader@getEffectiveSizeC[];
timeSeries=reader@getSizeT[];
slices = reader@getSizeZ[];
Table[
stack=Table[
ColorCombine[
Table[
index=reader@getIndex[z,c,t];
image=reader@openImage[index];
raster=image@getRaster[];
height=image@getHeight[];
width=image@getWidth[];
{type,colorspace}=BioFormatToMathematicaTypeAndColorSpace[image@getType[]];
Image[ArrayReshape[raster@getPixels[0,0,width,height,JLink`JavaNew["[I",width*height]],{height,width}], type],
{c,0,channels-1}
],
colorspace
],
{z,0,slices-1}
];
If[slices>1,Map[Image,stack,{-1}],Image@@stack],
{t,0,timeSeries-1}
]
]

ImportBioFormatMetadata to import image metadata.

ImportBioFormatMetadata[file_String]:=
JavaBlock@Module[
{ir=JLink`JavaNew[JLink`LoadJavaClass["loci.formats.ImageReader"]],reader,metastring},
ir@setId[file];
reader=BufferedImageReader`makeBufferedImageReader[ir];
metastring=reader@getCoreMetadataList[]@toString[];
Map[
First@StringReplace[#, StartOfString~~token__~~" = "~~value__~~EndOfString:> (StringTrim[token]-> If[StringMatchQ[value,NumberString],ToExpression[value],value])]&,
Map[Rest,StringSplit[StringSplit[StringTake[metastring,{2,-2}],", "],"\n"]],
{2}
]
]

An Import Example

Embryo cells provided by Kevin O’Connell from NIDDK at NIH. The image format is Openlab LIFF.

filePath=FileNameJoin[{NotebookDirectory[],"Data","embryo2","embryo2"}]

Group the list of grayscale images into groups of 3 and convert them to RGB color.

embryo=BlockMap[ColorCombine,First@ImportBioFormatImage[filePath],3];

Manipulate[ImageAdjust[embryo[[k]]],{k,1,Length[imgs],1}]

Obtain the meta-information of the image sequence.

First@ImportBioFormatMetadata[filePath]

Camera Connection

Utilize SKDs and frameworks to connect directly to a camera.

Devices Framework

FindDevices[]

devices=Devices["Camera"]

microscope = First[devices];

DeviceOpen[microscope]

img=DeviceRead[microscope,"Image"]

DeviceClose[microscope]

OpenCV via LibraryLink

LibraryLink C-code using OpenCV to implement the functions capture_microscopy_image and close_microscope.

#include "WolframLibrary.h"
#include "WolframImageLibrary.h"
#include "opencv2/opencv.hpp"

using namespace cv;

int iface = 0;
VideoCapture cap;


// Library Link //
 
 EXTERN_C DLLEXPORT mint WolframLibrary_getVersion ( ) {
       return WolframLibraryVersion;
   }

EXTERN_C DLLEXPORT int WolframLibrary_initialize ( WolframLibraryData libData) {
      cap = VideoCapture (iface); // open the default camera
         return 0;
  }

EXTERN_C DLLEXPORT void WolframLibrary_uninitialize ( WolframLibraryData libData) {
      try {
             cap.~VideoCapture ();
         }
        catch (Exception) {}
        return;
  }


//  Library Functions  //
 
 template < typename T >
  static void ifree_capture (void *out0, Mat & I) {
        T *out = reinterpret_cast < T * > (out0);
         
         // accept only char type matrices
          CV_Assert (I.depth () == CV_ 8 U);
        
        int channels = I.channels ();
        
        int nRows = I.rows;
        int nCols = I.cols * channels;
        
        if (I.isContinuous ())
         {
              nCols *= nRows;
              nRows = 1;
          }
         
         int i, j;
        uchar* p;
        for (i = 0; i < nRows; ++i)
         {
              p = I.ptr < uchar > (i);
              for (j = 0; j < nCols; ++j)
               {
                    out[i*nRows + j] = (T) p[j];
                }
          }
    }

EXTERN_C DLLEXPORT int close_microscope (WolframLibraryData libData, mint Argc,
                               MArgument *Args, MArgument res) {
      try {
             cap.~VideoCapture ();
         }
        catch (Exception) {}
        return 0;
  }


EXTERN_C DLLEXPORT int capture_microscopy _image (WolframLibraryData libData, mint Argc,
                                      MArgument *Args, MArgument res) {
      
      mint I0;
          I0 = MArgument_getInteger (Args[0]);
          
          if (iface != I0) {
                 iface = I0;
                 cap = VideoCapture (iface);
             }
            if (! cap.isOpened ()) {
                 iface = I0;
                 cap = VideoCapture (iface);
             }
            if (! cap.isOpened ())
                return - 1;
          
          Mat frame;
          cap >> frame; // get a new frame from camera
           
           mint length;
         MImage image_out = 0;
         void *data_out;
         int err = LIBRARY_FUNCTION _ERROR;
         imagedata_t type;
         WolframImageLibrary_Functions imgFuns = libData -> imageLibraryFunctions;
         
         imgFuns -> MImage_new2D (frame.cols, frame.rows, frame.channels (), MImage_Type _Bit8, MImage_CS _RGB, true, & image_out);
         
         data_out = imgFuns -> MImage_getRawData (image_out);
         
         if (data_out == NULL)
              goto cleanup;
         
         ifree_capture < raw_t _ubit8 > (data_out, frame);
         
         MArgument_setMImage (res, image_out);
         // return err;
        
        return 0;
    cleanup :
          imgFuns -> MImage_free (image_out);
        return err;
        
        // the camera will be deinitialized automatically in VideoCapture destructor
         return 0;
  }

First, you have to load this package.

Needs["CCompilerDriver`"];

This creates a full library path name.

lib = FileNameJoin[{NotebookDirectory[],"libMicroscopyIO.dylib"}]

The library was written as a Wolfram LinkLibrary. You can load functions from it using LibraryFunctionLoad:

capture=LibraryFunctionLoad[lib, "capture_microscopy_image",{Integer},LibraryDataType[Image|Image3D]];

Command to capture an image from the microscope:

CaptureMicroscopyImage[n_Integer:-1] := ColorCombine[Reverse@ColorSeparate[capture[n]],"RGB"];

Command to close the connection to the microscope:

CloseMicroscopyImage = LibraryFunctionLoad[lib, "close_microscope",{},{Integer}];

Brightness Equilibration

A typical image with uneven lighting:

img =

Determine the background region that is supposed to be evenly illuminated.

mask=ColorNegate@AlphaChannel@RemoveBackground[img]

Determine the flat-field correction by fitting a brightness model to the background.

flatfield=FindFlatFieldCorrection[img, Method-> {"Polynomial",4}, Masking-> mask]

Determine a RGB flat-field:

ColorCombine@Map[img|->FindFlatFieldCorrection[img, Method-> {"Polynomial",4}, Masking-> mask], ColorSeparate[img]]

Apply one flat-field correction to an image.

correctionBW=ImageLocalBrightnessAdjust[img, Method-> {"Polynomial",4}, Masking-> mask];

Apply separate flat-field corrections to each channel of an image.

correctionRGB=ColorCombine@Map[img|->ImageLocalBrightnessAdjust[img, Method-> {"Polynomial",4}, Masking-> mask],ColorSeparate[img]];

View the results:

TabView[{"Image"->img,"Adjusted"->correctionBW, "RGB Adjusted"->correctionRGB}]

Interactive Viewing

Images tend to be larger than the notebook. DynamicImage provides a panning and zooming pane to view large images.

img = ;

DynamicImage[img]

Interactive Viewing

Construct interactive tools to measure or manipulate image content:

Press ⌥ to more the circle. Move slider to vary the radius.

measurementCircle[Dynamic[center_],Dynamic[r_]]:={Yellow,Thick,Dynamic@Circle[center,r],Dynamic@Text[Style[StringJoin[ToString@Round@r,"px"],FontSize-> 18],Scaled@{0.1,0.9},Background-> RGBColor[0,0,0,0.5]]};

DynamicModule[
{center=ImageDimensions[img]/2,r},
Manipulate[
EventHandler[
DynamicImage[
img,
Epilog-> measurementCircle[Dynamic[center],Dynamic[r]]
],
{"MouseDown":>If[CurrentValue["OptionKey"],center=MousePosition["Graphics"]]},
PassEventsDown-> Dynamic[Not[CurrentValue["OptionKey"]]]
],
{{r,32},8,800},
FrameMargins->0
]
]

Interactive Viewing

Measure the length of the DNA by tracing it. Place a new point via ⌥-click. Delete line-points via ⇧-⌥-click.

dna=;

measurementLine[Dynamic[pts_],Dynamic[length_]]:={Yellow,Thick,Dynamic@Line[pts],Dynamic@Text[Style[StringJoin[ToString@Round@length,"px"],FontSize-> 18],Scaled@{0.1,0.9},Background-> RGBColor[0,0,0,0.5]]}

DynamicModule[
{pts={}, length = 0},
Panel@EventHandler[
DynamicImage[
dna,
Epilog-> measurementLine[Dynamic[pts],Dynamic[length]]
],
{"MouseDown":>If[CurrentValue["OptionKey"],If[CurrentValue["ShiftKey"],pts=Most[pts],AppendTo[pts,MousePosition["Graphics"]]];length=Total@Apply[EuclideanDistance,Partition[pts,2,1],{1}]]},
PassEventsDown-> Dynamic[Not[CurrentValue["OptionKey"]]]
]
]

Image Resolution

A microscope has a very shallow depth of field. Thicker samples are never entirely in focus.

The point spread function of a microscope usually exhibits a lateral resolution that is about 4-times higher than the resolution in the vertical direction.

Moving the focal plane through the microscopic sample reveals its depth.

SetDirectory[FileNameJoin[{NotebookDirectory[],"Data","FocusStack"}]];

stack = Map[file |-> ImageResize[Import[file],640],FileNames["*.bmp"]];

ListAnimate[stack, ImageSize-> 640]

Focus Stacking

One can apply focus stacking to overcome the shallow depth of field.
Simply assemble only those parts of a stack into one image that are in focus.
The latter can be detected by a Laplacian filter.

focusResponses=
Map[
ImageNorm@LaplacianGaussianFilter[#,{3,1}]&,
stack
];

depthVol=MaxChannelSuppression[focusResponses];

inFocus=ImageAdd@@MapThread[ImageMultiply,{stack,depthVol}]

3D Reconstruction

The focus stack also provides the information for a depth map.

We utilize only the reliable depth information in a focus stack.

confindenceMap=ImageAdd@@MapThread[ImageMultiply[Image`ImageDivide[#1,ImageAdd[ColorConvert[#2,"Grayscale"],0.2]],#3]&,{focusResponses,stack,depthVol}]

depthMap1=ImageAdjust[ImageAdd@@MapIndexed[ImageMultiply[#1,First[#2]]&,depthVol]]

depthMap2=ImageMultiply[depthMap1,Binarize[confindenceMap,0.02]]

Reducing noise by a median filter and fixing wholes in the depth map with a filling transformation.

depthMap3=FillingTransform@MedianFilter[depthMap2,3]

Displaying a 3D reconstruction.

ListPlot3D[
Reverse@ImageData[Blur@ImageResize[depthMap3,Scaled[1/8]]],
InterpolationOrder->1, 
PlotStyle->Texture[SetAlphaChannel[inFocus,Binarize[depthMap3,0.1]]] , 
Mesh-> False,
Lighting->{{"Ambient",White}},
Axes-> None,
BoxRatios->{1,ImageAspectRatio[depthMap3],1/3},
ImageSize-> 640
]

A Forensic Classification Example

Machine learning has a huge potential in the automatic interpretation of microscopic image data.

Here we train a classifier to distinguish the original handwriting with a pen from laser copies.

originals={,,,,,,,,,};

copied={,,,,,,,,,};

Train the classifier:

cf=Classify[<|"original"->originals,"copied"->copied|>]

Apply the classifier:

cf[]

cf[]

The Spikey at μm-Scale

Initializations

Get[FileNameJoin[{NotebookDirectory[],"ImageProcessing","BrightnessAdjust.wl"}]]

ImageNorm=img|->Image`ImagePower[ImageAdd@@ColorSeparate[Image`ImagePower[img,2]],1/2];

MaxChannelSuppression[list:{__Image}]:=
ImageMultiply@@@Partition[
Map[(Binarize[#, 0]&)@*ImageSubtract, Permutations[list, {2}]], 
Length[list] - 1
]