Seminar.nb

•Digital Image Processing Application Package

The Digital Image Processing application package was designed for engineers, scientists, biomedical professionals, and students who need versatile, comprehensive image processing software. The package contains an extensive collection of image transformation and manipulation functions as well as an interactive manual and tutorial to get new users started quickly. The package is completely integrated into Mathematica, therefore all image processing functions can be fully customized and new image processing algorithms easily developed and added.

This loads the package.

In[1]:=

This sets a useful global variable.

In[2]:=

•Functionality

The package adds over 160 new image processing and related functions to Mathematica's basic functionality.

- Color conversion between supported grayscale, RGB, HSV, and CMYK color models
- Extensive collection of point and area operators
- Image measurement functions
- Intensity profiles over arbitrary polygonal paths
- Histograms and co-occurrence matrices
- Fast integer-factor interpolation and decimation functions
- Resizing, padding, and merging
- Spatial transformations including rotation and arbitrary-order warping
- Area-of-interest processing over arbitrary polygonal regions
- Dozens of predefined filters and morphological operators
- Comprehensive number of filter design algorithms
- Fourier, cosine, Hadamard, and wavelet image transforms

•Structure of ImageData

Here we load a demonstration color image. Note that the default output form of an ImageData expression is the string "".

In[3]:=

Out[3]=

ImageData like other graphics primitives may be displayed using Show command.

In[4]:=

[Graphics:../HTMLFiles/index_7.gif]

The structure of ImageData is very simple. It consists of the raw image data, followed by two rules defining a pixel interleaving method and a color function. Consider the input form of rgb.

In[5]:=

Out[5]//InputForm=

ImageData[{{{255, 0, 0}, {0, 255, 0}, {0, 0, 255}, {255, 255, 255}, {0, 0, 0}},
  {{255, 0, 0}, {0, 255, 0}, {0, 0, 255}, {255, 255, 255}, {0, 0, 0}}, {{196, 225, 242}, {96, 204, 248}, {106, 122, 211},
   {219, 112, 131}, {144, 42, 55}}, {{202, 172, 172}, {111, 229, 152}, {227, 140, 232}, {175, 149, 115}, {249, 229, 3}},
  {{2, 168, 101}, {176, 116, 98}, {206, 37, 113}, {185, 107, 237}, {76, 76, 45}}}, PixelInterleave -> True,
ColorFunction -> RGBColor]

A small number of functions are available to extract the raw image data, switch between interleaving schemes, and split and merge ImageData expressions. Here are a few examples.

Here is the raw data.

In[6]:=

Out[6]//MatrixForm=

( ( 255 ) ( 0 ) ( 0 ) ( 255 ) ( 0 ) ) 0 255 0 255 0 0 0 255 255 0 ( 255 ) ( 0 ) ( 0 ) ( 255 ) ( 0 ) 0 255 0 255 0 0 0 255 255 0 ( 196 ) ( 96 ) ( 106 ) ( 219 ) ( 144 ) 225 204 122 112 42 242 248 211 131 55 ( 202 ) ( 111 ) ( 227 ) ( 175 ) ( 249 ) 172 229 140 149 229 172 152 232 115 3 ( 2 ) ( 176 ) ( 206 ) ( 185 ) ( 76 ) 168 116 37 107 76 101 98 113 237 45

This converts the ImageData expression to a color-interleaved (or planar) format.

In[7]:=

Out[7]//InputForm=

ImageData[{{{255, 0, 0, 255, 0}, {255, 0, 0, 255, 0}, {196, 96, 106, 219, 144}, {202, 111, 227, 175, 249},
   {2, 176, 206, 185, 76}}, {{0, 255, 0, 255, 0}, {0, 255, 0, 255, 0}, {225, 204, 122, 112, 42},
   {172, 229, 140, 149, 229}, {168, 116, 37, 107, 76}}, {{0, 0, 255, 255, 0}, {0, 0, 255, 255, 0},
   {242, 248, 211, 131, 55}, {172, 152, 232, 115, 3}, {101, 98, 113, 237, 45}}}, PixelInterleave -> False,
ColorFunction -> {RGBColor[#1, 0, 0] & , RGBColor[0, #1, 0] & , RGBColor[0, 0, #1] & }]

Here is a visualization of the example image in a planar format.

In[8]:=

[Graphics:../HTMLFiles/index_13.gif]

Out[8]=

Here we split image rgb into three single-color channels. Each channel is defined by an ImageData expression and retains the source color information.

In[9]:=

Out[9]//InputForm=

{ImageData[{{255, 0, 0, 255, 0}, {255, 0, 0, 255, 0}, {196, 96, 106, 219, 144}, {202, 111, 227, 175, 249},
   {2, 176, 206, 185, 76}}, PixelInterleave -> None, ColorFunction -> (RGBColor[#1, 0, 0] & )],
ImageData[{{0, 255, 0, 255, 0}, {0, 255, 0, 255, 0}, {225, 204, 122, 112, 42}, {172, 229, 140, 149, 229},
   {168, 116, 37, 107, 76}}, PixelInterleave -> None, ColorFunction -> (RGBColor[0, #1, 0] & )],
ImageData[{{0, 0, 255, 255, 0}, {0, 0, 255, 255, 0}, {242, 248, 211, 131, 55}, {172, 152, 232, 115, 3},
   {101, 98, 113, 237, 45}}, PixelInterleave -> None, ColorFunction -> (RGBColor[0, 0, #1] & )]}

•Example 1: Segmentation

Amplitude thresholding is one of many segmentation techniques. Two-level or binary thresholding changes a pixel value to one of two desired values (typically 0 and 1) depending on the setting of a threshold. In its simplest form, thresholding is a point-based operation that assigns the values of 0 or 1 to each pixel of an image based on a comparison with some global threshold value T.

Thresholding is an attractive early processing step because it leads to significant reduction in data storage and results in binary images that are simpler to analyze. Binary images are important in many computer vision applications and in the publishing and printing industry. Binary images permit the use of powerful operators for shape and structure-based analysis of image content. There are many techniques for converting multi-tone (n-bpp) images into two-tone (1-bpp). They all require the selection of a threshold value T. The selection of the best threshold value T is an interesting image processing problem. Here is an image segmentation example where we extract the yellow beans by applying the threshold operation to the individual color channels and use a masking operation to eliminate all non-yellow pixels from the example image. The respective thresholds were selected from an examination of the channel histograms.

In[10]:=

Out[10]=

This thresholds the beans image and displays the binary-valued RGB channels.

In[11]:=

[Graphics:../HTMLFiles/index_21.gif]

Here we form a mask that identifies all the yellow pixels in the beans image according to the following formula.

where are the individual color channels of the binarized image tmp.

This computes the mask: (1) the image is converted into planar format, (2) the green channel is complemented, and (3) the three binary-valued matrices are multiplied.

In[12]:=

Here we multiply the original image and the mask to obtain the final result.

In[13]:=

[Graphics:../HTMLFiles/index_26.gif]

Here I add some morphological post-processing to clean the resulting image.

In[14]:=

[Graphics:../HTMLFiles/index_28.gif]

In[15]:=

•Example 2: Region-of-Interest Processing

It is sometimes of interest to process a single sub-region of an image, leaving other regions unchanged. This is commonly referred to as region-of-interest (ROI) processing. Here is an example of a threshold operation performed on a ROI interactively selected by the user.

This loads the example image. The picture shows a fragment of the aortic wall as seen with a polarizing light microscope.

In[16]:=

[Graphics:../HTMLFiles/index_31.gif]

Here we define a region-of-interest. Image regions may be conveniently selected using the mouse and keyboard (see Get Graphics Coordinates... under the Input menu). Use the Copy and Paste commands to paste the recorded list of positions to any cell or expression in the notebook.

In[17]:=

Out[17]=

This shows the result of applying a threshold to the region defined by roi.

In[18]:=

[Graphics:../HTMLFiles/index_35.gif]

This returns the total number of pixels in the region-of-interest.

In[19]:=

Out[19]=

This counts the number of bright pixels in the region-of-interest. White pixels identify lipids in the aortic wall.

In[20]:=

Out[20]=

Here I find edges in the region-of-interest.

In[21]:=

[Graphics:../HTMLFiles/index_41.gif]

•Example 3: Burt-Adelson Pyramid

As a final demonstration, I compute the well-known Burt and Adelson image pyramid. Note that by using Mathematica's powerful functional programming constructs, the Gaussian and Laplacian pyramids are computed with "one-line" programs.

This loads an example image.

In[22]:=

This defines a commonly used generating kernel.

In[23]:=

Out[23]=

This constructs the Gaussian pyramid. The code shown repeatedly convolves an image with the generating kernel and decimates the result.

In[24]:=

Out[24]=

The images have dimensions that decrease by a factor of 4 at each level of the pyramid.

In[25]:=

Out[25]=

This shows the images in a graphics array.

In[26]:=

[Graphics:../HTMLFiles/index_50.gif]

The differences of images at successive levels of the Gaussian pyramid defines the Laplacian pyramid. To calculate the difference, the images at a lower level must be interpolated by a factor of 2 prior to subtraction.

This computes the pyramid.

In[27]:=

Out[27]=

This shows the six levels of the Laplacian pyramid.

In[28]:=

[Graphics:../HTMLFiles/index_54.gif]

The original image may be reconstructed from the Laplacian pyramid by reversing the previous steps. This interpolates and adds the images at successive levels of the pyramid beginning with the lowest level.

In[29]:=

Out[30]=

In[31]:=

[Graphics:../HTMLFiles/index_58.gif]

In[32]:=

•Example 4: Comet Assay

The single-cell gel electrophoresis (or comet) assay is a rapid and sensitive fluorescence microscopic method for the detection of primary DNA damage at the individual cell level. During the last 10 years, it has gained widespread acceptance for genotoxicity testing. Compared to other genotoxicity tests, the comet assay is a cheap and simple method well suited to analysis and measurement by means of computer-assisted techniques [from Introduction in "A public domain image analysis program for the single-cell gel-electrophoresis (comet) assay," by Christoph Helma and Maria Uhl, in Mutation Research, 466 (2000), 9-15].

See http://cometassay.com for additional information and references.

Here I load an example image (data courtesy Blake Whitaker, University of Southern Maine).

In[2]:=

$dataDir = ToFileName[{"d:\\mjkcc\\Data"}] ; img = ImageRead[dataDir <> "comet03.tif"] ;$

This defines the regions-of-interest as user selected rectangular areas in the graphical coordinate system of the Mathematica frontend. The coordinate pairs represent the lower-left and upper-right vertices of the bounding rectangles for each of six regions in the image.

In[3]:=

Out[4]=

Here is a small version of the original image data showing the regions-of -interest.

In[5]:=

[Graphics:../HTMLFiles/index_64.gif]

Out[5]=

I now proceed with extracting the comets from the original image. (Note for new users: the transposition and reversal of the coordinate pairs is needed in particular to satisfy the syntax of the ImageTake function, the /@ symbol represents the Map function, and the symbol & denotes a pure function.)

In[6]:=

Out[6]=

This shows the individual comets.

In[7]:=

Out[7]=

The original images are in 3-channel RGB format. Clearly the green channel signal is strongly dominant. This can be easily verified by displaying the planar format of an example comet.

In[8]:=

Out[8]=

A closer examination reveals the presence of weak signals in the red and blue channels, especially in the tail section of the comet (interestingly, this effect is not discussed in the scientific literature). Since all computer-aided quantification results are traditionally presented in the context of gray-scale images, so we proceed with a color transformation.

In[9]:=

Out[9]=

I now display the histograms of the example comet images.

In[10]:=

In[11]:=

[Graphics:../HTMLFiles/index_78.gif]

Out[11]=

I now proceed with a typical calculation. For example, here it is of interest to obtain the area of the head of the comet ("AH"), the total density of the head ("DNAH") and the mean intensity of the head ("MH"). Usually the same set of measurements is obtained for the tail and the two results are contrasted. In the literature on comet assays the head is identified by the pixels with the highest intensity. Here, as is commonly done, I use a threshold value that is 85% of the maximum intensity in the image.

In[12]:=

Out[12]=

I now binarize the images, adding a morphological filtering operation to smooth the borders of the head regions.

In[13]:=

Out[13]=

In[14]:=

Out[14]=

I now continue with identifying the regions and measuring the required properties, additionally including calculations of the centroids and enclosing circles.

In[15]:=

In[16]:=

Out[16]=

In[17]:=

Out[17]=

In[18]:=

Out[18]=

This displays the results in a tabular fashion.

In[19]:=

Out[19]//TableForm=

image #

area

density

mean intensity

305

285

330

45191

42293

49066

148.1672131147541`

148.39649122807018`

148.6848484848485`

Here I calculate the graphical elements that I will use to annotate the comet images.

In[20]:=

Out[20]=

In[21]:=

Out[22]=

Here is an alternative visualization of a comet image.

In[23]:=

[Graphics:../HTMLFiles/index_101.gif]

Out[23]=

Finally, I calculate and display another property of a comet, the so-called comet profile. It is a plot of the cross-sectional density (density measured perpendicularly to the length-wise dimension of the comet).

In[24]:=

[Graphics:../HTMLFiles/index_104.gif]

Out[25]=

Converted by Mathematica (July 29, 2002)