Digital Image Processing Version 2.0 - What's New?

A short tour of new features and major enhancements

Mariusz Jankowski

University of Southern Maine

Astract

The soon to be released version 2.0 of the Digital Image Processing application package will feature a substantial number of new functions and important enhancements to many existing functions, both in speed and capabilities, while retaining ease of use and maintaining compatibility with the latest versions of Mathematica. This talk will present several of these new capabilities and enhancements.

Overview - key new features

Usability

ImageData object has been generalized to allow representation of image data with arbitrary number of channels and unspecified color space.

Color space transformations

The following frequently requested transformations have been added: CIE XYZ, Lab, and Luv to/from RGB.

Image segmentation

New additions include watershed transform, distance transform, region labeling, Canny and Shen-Castan edge detectors, new threshold operators and derivative filters.

Power spectrum

Several Fourier transform methods for computing the power spectrum of an n-dimensional array were added including the classic periodogram, averaging periodogram (Welch), modified averaging periodogram (Bartlett) and the smoothing periodogram (Blackman and Tukey).

Image restoration

Several of the best-known image restoration algorithms including inverse, minimum mean square error (Wiener), constrained least squares, geometric mean and Lucy-Richardson.

Java connectivity

A translator for Java's BufferedImage class and the Java Advanced Imaging (JAI) library from Sun have been added.

Usability

<< ImageProcessing`

ImageData object has been generalized to allow representation of image data with arbitrary number of channels and a generic color space (using the new option value None for ColorFunction). The function ToGenericColor may be used to simplify the construction of image objects.

ToGenericColor[ToChannels[rgb = ImageRead["rgb.tif"]]] // 
      InputForm

Since there is no color space associated with the sample data in these channels they by default display as an array of GrayLevel rasters.

Show[Graphics[%]];

Usability (con't)

Here is an example. This reads a DICOM image with eight time-lapsed ultrasound images. I combine these images into a single 8-channel ImageData object so that they can all be conveniently manipulated as a single entity.

\(ImageRead[
      ToFileName[{$HomeDirectory, "My Documents", "Data"}, "US-mono2-8-8x-execho.dcm"]]\)

ToGenericColor[RawImageData /@ %]

This shows the effect of thresholding the channels with a single threshold value.

Show[Graphics[Threshold[%, 50]]];

Color space transformations

As a consequence of the generalization of the ImageData object it is now possible to construct image objects in any desider color space. Three well known and frequently requested transformations have been added: CIE XYZ, Lab, and Luv to/from RGB. The following returns three CIELABColor images uniformly sampled in the color plane defined by the a and b color coordinates with three different lightness values.

cie = 
        Table[ToCIELABColor[
            Table[{k, i, j}, {i, -100, 100, 2}, {j, -100, 100, 
                2}]], {k, 25, 75, 25}];

This maps the cie images into the rgb color space and displays the result.

rgb = ToRGBColor[#, GammaCorrection -> True] & /@ 
          cie;

DisplayTogetherArray[Show[Graphics[#, PlotRange -> {0, 1}], 
              Axes -> True, AxesLabel -> {b, a}, 
              Ticks -> {{{25, -50}, {50, 0}, {75, 
                      50}}, {{25, -50}, {50, 0}, {75, 50}}}] & /@ 
          rgb];

Image segmentation

Thresholding

Faster threshold implementation, new HysteresisThreshold and new methods in OptimumThreshold.

Linear filtering

New smoothing and differentiating filters.

Edge detection

Implementations of edge detection techniques following Canny and Shen - Castan.

Region labeling

The classic binary image connected components labeling algorithm (in Java and Mathematica).

Watershed transform

An implementation of the well-known Vincent Soille algorithm with extension to allow user selected markers for improved segmentation (in Java and Mathematica).

Optimum thresholds

The OptimumThreshold function now implements four different methods of threshold selection.

\(\((wing = ImageRead["wing.tif"]; 
      h = ImageHistogram[
          wing]; 
\(OptimumThreshold[h, 
            Method -> #] &\) /@ {"otsu", "kittler", "kapur", "mse"})\)\)

DisplayTogetherArray[Graphics[
              Threshold[wing, #]] & /@ %];

Here is a second example using a different type of image.

(cam = ImageRead["cameraman.tif"]; 
      h = ImageHistogram[cam];)

\(\(DisplayTogetherArray[\(Graphics[
              Threshold[
                cam, #]] &\) /@ \((\(OptimumThreshold[h, 
                  Method -> #] &\) /@ {"otsu", "kittler", "kapur", "mse"})\)];\)\)

New smoothing and differentiating filters

The new smoothing and differentiating filters lay at the foundations of two new edge detectors added to the package. The derivative of the Gaussian-shaped impulse response closely matches the optimal edge filter originally presented by Canny and the derivative of an exponential was proposed by Shen and Castan as an alternative "best" filter under a modified set of optimality conditions.

SetOptions[LinePlot, PointStyle -> PointSize[0.035], 
        Axes -> False, Ticks -> None, 
        PlotRange -> All];

DisplayTogetherArray[{{LinePlot[GaussianFilter[31]], 
            LinePlot[ExponentialFilter[31]]}, {LinePlot[
              DerivativeOfGaussianFilter[31, 1]], 
            LinePlot[DerivativeOfExponentialFilter[31, 1]]}}];

Canny edge detection

Edge detection according to Canny is a sequence of three operations: linear filtering with two orthogonal directional derivative filters, non-maximum suppression and hysteresis thresholding. All the necessary filters and operators are now available.

lena = ImageRead["lena.tif"];

This shows edge magnitudes using the Canny filters with and without non-maximum suppresion.

DisplayTogetherArray[{Graphics[
            EdgeMagnitude[lena, CannyFilter[17, 17]]], 
          Graphics[
            edge = EdgeMagnitude[lena, CannyFilter[17, 17], 
                NonMaxSuppress -> True]]}];

Here I use standard and hysteresis thresholding to obtain edge maps.

DisplayTogetherArray[{Graphics[Threshold[edge, 1]], 
          Graphics[Threshold[edge, 10]], 
          Graphics[HysteresisThreshold[edge, {1, 10}]]}];

Shen-Castan edge detection

Shen and Castan propose a symmetric exponential filter and its derivative for edge detection and describe a recursive realization of their filters. Here, the first derivative filter is implemented as a finite-sized kernel in analogy with other directional derivative filters in the package.

lena = ImageRead["lena.tif"];

This shows edge magnitudes using the Shen-Castan filters with and without non-maximum suppresion.

DisplayTogetherArray[{Graphics[
            EdgeMagnitude[lena, ShenCastanFilter[17, 17]]], 
          Graphics[
            edge = EdgeMagnitude[lena, ShenCastanFilter[17, 17], 
                NonMaxSuppress -> True]]}];

Here I use standard and hysteresis threholding to obtain edge maps.

DisplayTogetherArray[{Graphics[Threshold[edge, 1]], 
          Graphics[Threshold[edge, 8]], 
          Graphics[HysteresisThreshold[edge, {1, 8}]]}];

Region labeling

Connected components labeling and analysis is a relatively simple grouping algorithm that has been used for many years to isolate, analyse and identify disjoint regions in a binary image. The labeling operation assigns a unique name or number to all 1-pixels that belong to the same connected component of the image. As a result of the labeling, individual components can be extracted from the image programmatically and therefore be available for further processing and analysis. In the example below, I label the connected components and mark in red one randomly selected component.

books = ImageRead["books.tif"];

DisplayTogetherArray[{Graphics[bin = Threshold[books, 175]], 
          Graphics[tmp = RegionLabel[bin]], 
          Graphics[
            RegionFill[ToRGBColor[tmp], {Max[tmp], 0, 0}, 
              Where[tmp, Random[Integer, {1, Max[tmp]}]]]]}];

Watershed transformation

The watershed transform is a popular segmentation technique based on the dual concepts of flooding and grayscale images as topographic reliefs. The flooding process forms so-called catchment basins that grow in size as the water level rises. In the field of topography the lines of separation between nearby basins are called watersheds. The construction of catchment basins and watershed lines in digital images often leads to meaningful segmentation results. Here is an example image from http://cmm.ensmp.fr/~beucher/wtshed.html.

\(\(img = 
        ImageRead[
          ToFileName[{$HomeDirectory, "My Documents", "Data"}, "meyer.tif"]];\)\)

This calculates the watershed transform and extracts the watershed lines.

tmp = WatershedTransform[img];  wsheds = Where[tmp, 0];
pts = 
        MapThread[
          Point[Flatten[{Reverse[#1], #2}]] &, {wsheds, 
            Extract[img[[1]], 
              wsheds]}];

DisplayTogetherArray[{ListPlot3D[
            img[[1]], 
            Mesh -> False], Graphics[tmp], 
          DisplayTogether[{ListPlot3D[
                img[[1]], 
                Mesh -> False], 
              Graphics3D[{Hue[0], PointSize[0.02], pts}]}]}];

Watershed transformation

Here is an example of the application of the watershed transform to the segmentation of a comet assay. The single-cell gel electrophoresis (or comet) assay is a rapid and sensitive fluorescence microscopy method for the detection of primary DNA damage at the individual cell level. For basic measurement purposes it is of interest to isolate the "head" regions of the individual cells. This will be done using the marker-driven watershed transform.

\(\(comet = 
        ImageTake[
          ImageRead[
            ToFileName[{$HomeDirectory, "My Documents", "Data"}, "comet03.tif"]], {40, \(-1\)}];\)\)

Note that a direct approach doesn't work due to significant over-segmentation.

grn = 
        255 - ToChannels[
              comet][[2]];
grad = 
        Round[EdgeMagnitude[grn, 
            ShenCastanFilter[31, 31]]];
DisplayTogetherArray[{Graphics[comet], 
          Graphics[WatershedTransform[grn]], 
          Graphics[WatershedTransform[grad]]}];

Here are the user-selected seed points for the marker-based watershed algorithm.

seeds = 
        Reverse /@ 
          Round[{{53.9766, 184.766}, {38.9141, 47.1947}, {128.285, 
                114.474}, {240.752, 149.62}, {277.906, 118.491}, {224.685, 
                33.1364}}];

result = 
        WatershedTransform[Round[grad], seeds];
border = Point /@ Reverse /@ Where[result, 0];

DisplayTogetherArray[{ShowRegionOfInterest[grad, 
            Point /@ Reverse /@ seeds, 
            PointStyle -> {Hue[0], PointSize[0.02]}], Graphics[result], 
          ShowRegionOfInterest[comet, border, 
            PointStyle -> Hue[0]]}];

Power spectrum

Function PowerSpectrum implements three classic power spectrum estimation techniques: direct DFT calculation (default), Welch's method which computes the periodograms of overlapping segments of data and averages them, and the Blackman-Tukey method in which the DFT of a windowed autocorrelation sequence of the data is computed.

img = 
        ImageRotate[
          ToGrayLevel[
            Table[(UnitStep[Sin[(π/64.)  2  i]] - 
                    0.5) Cos[(π/64.)  15.  j], {i, 128}, {j, 
                128}]], 10 Degree];

SetOptions[ListDensityPlot, Mesh -> False, Axes -> False, 
        Frame -> False];
DisplayTogetherArray[{Graphics[img], 
          ListDensityPlot[PowerSpectrum[img]], 
          ListDensityPlot[
            PowerSpectrum[img, {64, 64}, {48, 48}, Method -> "Welch"]], 
          ListDensityPlot[
            PowerSpectrum[img, {64, 64}, 
              Method -> "BlackmanTukey"]]}];

Image restoration - non-iterative methods

The group of non-iterative methods includes inverse filtering, minimum mean square error (Wiener) filtering, generalized Wiener filtering method and constrained least squares. All methods require knowledge of the distortion and a variety of other noise and image characteristics. Here is an example of non-iterative restoration of an example image distorted by motion blur and additve Gaussian noise.

<< Statistics`NormalDistribution`

img = ImageRead["circuit.tif"];
blurred = DiscreteConvolve[img, {BoxFilter[15]}];
noise  =  
        ToGrayLevel[
          RandomArray[NormalDistribution[0, 5. ], 
            ImageDimensions[blurred]]];
noisyblurred  =  noise  +  blurred;

resInverse  =  
        InverseRestoration[
          noisyblurred, {BoxFilter[15]}];
resWeiner  =  
        WienerRestoration[noisyblurred, {BoxFilter[15]}, 
          (10^(-7))];
resConstrained  =  
        ConstrainedRestoration[noisyblurred, {BoxFilter[15]},  
          (10^(-3))];

DisplayTogetherArray[{Graphics[noisyblurred], Graphics[resInverse], 
          Graphics[resWeiner], Graphics[resConstrained]}];

Image restoration - iterative methods

The iterative image restoration functions include the standard and accelerated versions of the Lucy and Richardson (L-R) algorithm. Briefly, the L-R algorithm is an iterative technique used heavily for the restoration of astronomical images in the presence of Poisson noise. The algorithm requires a good estimate of the degradation process for accurate estimation. First, here is a 1D example.

f  =  {6,  7,  2, 5, 1};
h  =  {1, 1, 1}/3. ;
g  =  DiscreteConvolve[f,  h]

tmp = 
        LucyRichardsonRestoration[g,  h, 
          LucyRichardsonIterations -> 100, 
          LucyRichardsonAccelerated -> False, 
          LucyRichardsonVerbose -> True];
tmpAccelerated = 
        LucyRichardsonRestoration[g,  h, 
          LucyRichardsonIterations -> 50, 
          LucyRichardsonVerbose -> True];

DisplayTogetherArray[{DisplayTogether[ListPlot[#, 
                  PlotJoined  ->  True, 
                  PlotRange  ->  {0,  9}, 
                  AspectRatio -> 1/2] & /@ Transpose[tmp], 
            Graphics[{Hue[0], PointSize[0.02], Point /@ Thread[{100, f}]}]], 
          DisplayTogether[ListPlot[#, PlotJoined  ->  True, 
                  PlotRange  ->  {0,  9}, 
                  AspectRatio -> 1/1] & /@ Transpose[tmpAccelerated], 
            Graphics[{Hue[0], PointSize[0.04], 
                Point /@ Thread[{50, f}]}]]}];

Here is an image restoration example. The distortion is horizontal motion blur.

<< Statistics`NormalDistribution`

img = ImageRead["circuit.tif"];
blurred = DiscreteConvolve[img, {BoxFilter[15]}];
noise  =  
        ToGrayLevel[
          RandomArray[NormalDistribution[0, 3.5], 
            ImageDimensions[blurred]]];
noisyblurred  =  noise  +  blurred;

DisplayTogetherArray[{Graphics[img], Graphics[noisyblurred], 
          Graphics[
            LucyRichardsonRestoration[noisyblurred,  {BoxFilter[15]},  
              LucyRichardsonIterations  ->  8]]}];

Java connectivity

BufferedImage translator

A translator for Java's BufferedImage class and Sun's Java Advanced Imaging (JAI) library have been added.

The BufferedImage class was introduced in JDK 1.2 and encapsulates information about an image including the pixel values, the particular organization of the data and the color format.

Java Advanced Imaging (JAI)

The Java Advanced Imaging API provides a set of object-oriented interfaces that support a simple, high-level programming model which allows developers to manipulate images easily. Examples of image processing techniques available in JAI range from simple operations such as contrast enhancement, cropping, and scaling to more complex operations such as geometric warping and frequency domain processing. There are approximately 80 operators, many with native support for optimized performance.

JAI ImageIO

JAI Image I/O Tools provides reader, writer, and stream plug-ins for the Java Image I/O Framework and Image I/O-based read and write operations for Java Advanced Imaging. Reader-writer plug-ins are supplied for the BMP, JPEG, JPEG 2000, PNG, PNM, Raw, TIFF, and WBMP image formats.

For details see "Implementation of BufferedImage Translator" presented at the 2003 Mathematica Developer Conference.

JAI-enhanced functions

JAI-enhanced functions in the DIP package include BlockMedian, Resize, ImageRotate, ImageWarp, GrayscaleErode, GrayscaleDilate. Here is a timing comparison of versions 1.1 and 2.0.

v1 = {0.609, 2.375, 1.453, 0.828, 
          0.844};
v2 = {0.250, 0.188, 0.547, 0.234, 0.250};

\(TableForm[Transpose[{v1, v2, v1/v2}], 
      TableHeadings -> {{"median", "resize", "rotate", "dilate", "erode"}, {"v1.1", "v2.0", "speedup"}}]\)

Table 1. Timing comparison for several functions in DIP versions 1.1 and 2.0.

Reading JPEG2000 files

Here is an example of using the new ImageIO class recently made available to the Java community to read a file format not yet supported by Mathematica. (Note: the ImageIO package must be obtained from Sun and the included *.jar and *.dll files added to the ImageProcessing Java subdirectory.)

<< JLink`

InstallJava[];

LoadJavaClass["javax.imageio.ImageIO"];

This reads and displays a jpeg2000 file.

\(\(f = 
        JavaNew["java.io.File", 
          ToFileName[{$HomeDirectory, "My Documents", "Data"}, "lena.jp2"]];\)\)
bi = ImageIO`read[f]

Show[
        Graphics[img = BufferedImageToImageData[bi], 
          ImageSize -> 200]];

Summary

In summary, there are many new, exciting and useful additions and enhancements in the Digital Image Processing package. Overall there are approximately 40 new functions with possibly as many of the existing functions modified or enhanced in some manner. The Mathematica code size grew by approximately 36%.

Future work:

- segmentation
- iterative restoration
- blind deconvolution
- color
- graphical helpers using GUIKit