Image Processing in Mathematica 4

This how-to demonstrates how Mathematica 4 can be used to do image analysis tasks. In addition to many common operations, the ListInterpolation function can be used to do new and exiting manipulations and analyses.

Note: To save bandwidth, the original imported picture is not shown in every subsection. Click here to see it

Importing, displaying, and exporting graphics files

The Import command provides Mathematica users with an efficient and intuitive way to read files in most standard graphics formats, including bitmap, JPEG, GIF, TIFF, and numerous platform specific formats. In many cases, Mathematica even automatically detects the correct graphics format; otherwise, users can specify it themselves. Note that one has to use double backslashes in the file path. To interactively select the file, use the GetFilePath command in the Input menu. This command opens a standard file dialog and returns a correctly formatted path to the file at the current text cursor position.

dog  =  
        Import["C:\\desk\\40 Demos\\schlapp2.jpg"];

The graphic can now be displayed using Show,

Show[dog]

and it can be exported in a different format. For a complete list of import and export formats, consult the documentation for  Import and Export.

Export["C:\\windows\\desktop\\dog.tif", dog, "TIFF"]

tifdog = Import["C:\\windows\\desktop\\dog.tif"];

Show[tifdog]

Color conversions

ColorOutput converts Mathematica graphics between different color models. The three standard models are RGBColor, CMYKColor and GrayLevel. Note that ColorOutput changes the image representation, not just the rendering.

Show[dog,  ColorOutput -> GrayLevel];

Accessing image data

Mathematica handles many different kinds of objects: mathematical formulas, lists, and graphics, to name but a few. Although they often look very different, Mathematica represents all in one uniform way: they are all expressions. The imported graphics files are, for example, stored as a Graphics object.

The object f in an expression f[x, y, … ] is known as the head of the expression. You can extract it using Head. Particularly when you write programs in Mathematica, you will often want to test the head of an expression to find out what kind of thing the expression is.

Both the imported JPEG and TIFF files are, of course, Graphics objects.

Head[dog]

Head[tifdog]

The easiest way to access the pixel data is to change the head of the expression to a Mathematica list and extract the array of pixel values. The following command replaces the Graphics head with a List head. Explaining the mechanism is beyond the scope of this document. If you are interested, consult the documentation for Replace.

dog  = dog /. Graphics -> List;

The list contains an array of pixel values and some supplemental information on image size and color model. Short gives an overview.

Short[dog, 12]

The pixel values are always the first element of the data structure. They can now be treated like any other list in Mathematica. A RGB image, for example, is represented as an array of {r, g, b} triples.

f(x, y)=(({{r(0, 0), g(0, 0), b(0, 0), ..., r(0, N - 1), 
                  g(0, N - 1), b(0, N - 1)}, {r(1, 0), g(1, 0), b(1, 0), ⋯, r(1, N - 1), g(1, N - 1), 
                  b(1, N - 1)}, {⋮, ⋱, ⋮}, {r(M - 1, 0), g(M - 1, 0), b(M - 1, 0), ⋯, r(M - 1, N - 1), g(M - 1, N - 1), 
                  b(M - 1, N - 1)}})).

The next line extracts the array of pixel values. To learn more about matrix manipulation, consult Section 1.8.4. of The Mathematica Book or the Matrix Manipulation how-to.

pixelvalues = dog[[1, 1]];

As you can see, the color values are now stored in a 302×291×3array.

Dimensions[pixelvalues]

For example, here is a ListDensityPlot of the red channel of the image. All is a new function in Mathematica 4 that simply instructs Mathematica to take all elements in the respective dimension of the array.

ListDensityPlot[pixelvalues[[All, All, 1]], 
        Mesh -> False];

Simple math

All Mathematica functions that can be applied to lists can also work on the image data. For example, the following two lines are two different ways of inverting a grayscale image.

ListDensityPlot[-pixelvalues[[All,  All,  1]], 
        Mesh -> False];

ListDensityPlot[256 - pixelvalues[[All, All, 1]], 
        Mesh -> False];

The same procedure works for color images, too. Using Apply is a shortcut to convert all three color values in one step. Also, RGBColor expects its values to be between 0 and 1, so pixelvalues gets divided by 256.

Show[
        Graphics[
          RasterArray[
            Apply[(RGBColor[#1, #2, #3] &), 
              1 - pixelvalues/256. , {2}]]], 
        AspectRatio -> 1];

Histograms

Although Mathematica's histogram functions in Graphics`Graphics` only work for small datasets, histograms can be easily created by using BinCounts and ListPlot.

BinCounts is part of the Statistics`DataManipulation` package that ships with Mathematica.

Read in the package.

<< Statistics`DataManipulation`

Count the occurrences of pixel values in the red channel of the image.

histo = 
        BinCounts[
          Flatten[pixelvalues[[All, All, 1]]], {0, 255}];

ListPlot[histo, PlotJoined -> True]

Histogram operations

Because the image has been converted into a Mathematica list, any Mathematica function can be applied to parts of the list. This makes histogram operations trivial. For example, add a reddish tint by taking the square root of the red channel. Note that the color values for RasterArray have to be scaled between 0 and 1. The #s in the code are pure functions to learn more about them, see The Mathematica Book.

In the following picture, the red channel is decimated by 70 percent.

Show[
        Graphics[
          RasterArray[
            Apply[(RGBColor[0.3*#1, #2, #3] &), 
              pixelvalues/256. , {2}]]], AspectRatio -> 1];

The next picture shows the effect of taking the square root of the color values in the red channel. Because the values are scaled between 0 and 1, this gives the image a reddish tint.

Show[
        Graphics[
          RasterArray[
            Apply[(RGBColor[Sqrt[#1], #2, #3] &), 
              pixelvalues/256. , {2}]]], AspectRatio -> 1];

Selecting pieces of images

Slices

Once the pixel information is converted into a Mathematica list, it can be operated on like any other Mathematica list of lists. In the following picture, the three color channels are extracted by slicing the 266×276×3 data set along the z axis.

Click here to learn more about selecting pieces of lists.

Show[
        GraphicsArray[{ListDensityPlot[pixelvalues[[All, All, 1]]/256, 
              AspectRatio -> 1, 
              ColorFunction -> (RGBColor[#, 0, 0] &), opts1], 
            ListDensityPlot[pixelvalues[[All, All, 2]]/256, 
              AspectRatio -> 1, 
              ColorFunction -> (RGBColor[0, #, 0] &), opts1], 
            ListDensityPlot[pixelvalues[[All, All, 3]]/256, 
              AspectRatio -> 1, 
              ColorFunction -> (RGBColor[0, 0, #] &), opts1]}], 
        AspectRatio -> Automatic, 
        DisplayFunction -> $DisplayFunction];

Or, take a slice through row 200 of the red channel.

ListPlot[pixelvalues[[200, All, 1]], 
        PlotJoined -> True];

To visualize this slice, set all the pixels along this line to white. Copying the array into a new variable first preserves the original image.

new = pixelvalues;

new[[200, All, 1]] = 256;

ListDensityPlot[Evaluate[new[[All, All, 1]]], 
        Mesh -> False];

The easiest way to plot multiple slices at once is to use the MultipleListPlot standard package.

<< Graphics`MultipleListPlot`

MultipleListPlot[pixelvalues[[200, All, 1]], 
        pixelvalues[[200, All, 2]], pixelvalues[[200, All, 3]], 
        PlotJoined -> True, SymbolShape -> None, 
        PlotStyle -> {red, green, blue}];

Because ListDensityPlot and Raster only work with two dimensional arrays, one has to use RasterArray to put markers on the full color image.

new = pixelvalues;

new[[200, All, All]] = {256, 0, 0};

Show[
        Graphics[
          RasterArray[Apply[(RGBColor[#1, #2, #3] &), new/256. , {2}]]], 
        AspectRatio -> 1];

Area selection

Selecting areas of a graphic is also straightforward. To find the values for Take, click on the graphic and move the mouse while holding down the Control key. The pointer location will be displayed in the lower left corner of the window. To select multiple points, click on all that interest you, then copy, unselect the graphic and paste. Mathematica will return a list of all the locations you selected. For more information on how to select coordinates in Mathematica graphics, see the "How do I extract coordinates" how-to.

This line of code selects a 51-by-31 pixel rectangle of the picture.

righteye = Take[pixelvalues, {150, 180}, {100, 150}];

Dimensions[redeye]

Show[
        Graphics[
          RasterArray[
            Apply[(RGBColor[#1, #2, #3] &), 
              righteye/256. , {2}]]]];

The next line sets the red channel to 0.

greeneye = righteye /. {x_, y_, z_} -> {0, y, z};

Show[
      Graphics[RasterArray[
          Apply[(RGBColor[#1, #2, #3] &), 
            greeneye/256. , {2}]]]]

Putting the results back into the original image is also pretty easy and straightforward. First copy the image into a new variable to preserve the original image

new = pixelvalues;

Then replace the corresponding values with the transformed ones.

new[[Range[150, 180], Range[100, 150], All]] = 
        greeneye;

Show[
        Graphics[
          RasterArray[Apply[(RGBColor[#1, #2, #3] &), new/256. , {2}]]], 
        AspectRatio -> 1];

Thresholding

To select pixels based on some property, simply use If with the appropriate cutoff value. For example, the following line selects all pixels with an R value greater than 0.5×256.

Some more information on selecting parts of lists, see the Matrix Manipulation how-to.

Show[
        Graphics[
          RasterArray[
            Apply[(If[#1 > 0.5, RGBColor[1, 0, 0], 
                    RGBColor[#1, #2, #3]] &), pixelvalues/256. , {2}]]], 
        AspectRatio -> 1];

Convolution and correlation

The new Mathematica 4 functions ListConvolve and ListCorrelate make a number of advanced image processing tasks easy to do. The following graphic shows the effect of applying a Sobel edge filter to the image by convolving the image with two 3×3 kernels.

ListDensityPlot[ListConvolve[(({{-1, -2, -1}, {4, 0, 5}, {1, 2, 1}})),ListConvolve[(({{1, 0, 1}, {0, 0, 2}, {-1, 0, 1}})),pixelvalues[[All, All, 1]]/256. ]],Mesh -> False];

List interpolation and data recovery

An interesting feature of Mathematica is its ability to rapidly calculate list interpolations of any degree in arbitrary dimensions. This feature opens up new ways to manipulate images and recover information lost by the discretizing of data.

The Doghill mountains

Before doing anything useful with list interpolation, let's try something fun first. Perform a first-order list interpolation over the data set and plot the color values in 3D. Using the new real time 3D graphics, one can then rotate the graphic and see the image pop into and out of view.

Timing[interred = 
          ListInterpolation[pixelvalues[[All, All, 1]], {{0, 1}, {0, 1}}, 
            InterpolationOrder -> 1];]

<< RealTime3D`

Plot3D[interred[x, y], {x, 0, 1}, {y, 0, 1}, Mesh -> False, 
        Axes -> False,  Boxed -> False, 
        PlotPoints -> 200];

Data recovery

Let's go back to the small piece of the image selected in "Selecting pieces of images".

righteye = Take[pixelvalues, {150, 180}, {110, 140}];

ListDensityPlot[righteye[[All, All, 1]], 
      Mesh -> False]

As you can see, the image above consists of only 51×31 data points, so the resolution is not terribly high. List interpolation provides a convenient way to recover more features.

The following command creates a list interpolation over the red channel of the image and scales the resulting function between 0 and 1 in the x and y coordinates.

interpol1 = 
        ListInterpolation[
          righteye[[All, All, 1]], {{0, 1}, {0, 1}}];

DensityPlot[interpol1[x, y], {x, 0, 1}, {y, 0, 1}, 
        PlotPoints -> 100,  Mesh -> False];

Mathematical operations

One can do all kinds of fancy image manipulations and transformations using ListInterpolation. For example, differentiation find changes in color values. Two dimensional differentiation also creates a nice embossing.

DensityPlot[Evaluate[D[interred[x, y], x]], {y, 0, 1}, {x, 0, 1}, 
        PlotPoints -> {306, 291}, Mesh -> False, 
        Frame -> False];