Wolfram Blog
Piotr Wendykier

Gigapixel Images in Mathematica

April 29, 2013 — Piotr Wendykier, Mathematica Algorithm R&D

Professional cameras offer a resolution of 50 megapixels and more. In addition, projects like GigaPan allow one to create gigapixel panoramas with billions of pixels. How can we process these images on a desktop computer with 8 GB of RAM?

One of Mathematica 9′s new and exciting features is out-of-core image processing. What does the out-of-core term really mean? It is a way to process very large images that are too big to fit into main memory. Let’s say we have a machine with 8 GB of RAM, and let’s assume that Mathematica can use up to 7.2 GB of that memory (the remaining 0.8 GB will be used by the operating system). Freshly started, Mathematica 9 on Windows 8 takes up about 200 MB of memory, so the kernel can use about 7 GB of RAM. What is the maximal size of the image that we can load into the kernel (we don’t want to visualize it at this point)? If we assume that the image is in the RGB color space and a single byte encoding, then the following formula gives a maximal width (and height) of an image that can be loaded at once into the memory:

maxMemoryInBytes0 = 7 1024^3 ; maxWidth0 = Floor[Sqrt[maxMemoryInBytes0/3.0]]

50053

50 k by 50 k is a very large image, but remember that we haven’t left any memory for processing and visualization. If we want to perform some operations on that image and visualize the results, then the realistic memory limit would be something around 0.5 GB. Why only 0.5 GB? We should leave at least 2 GB for the front end, and we have to keep in mind that any nontrivial image-processing algorithm requires additional storage for temporary objects. So, when we restrict a single image storage to 0.5 GB, then the maximal dimensions would be:

maxMemoryInBytes1 = 0.5 1024^3 ; maxWidth1 = Floor[Sqrt[maxMemoryInBytes1/3.0]]

13377

The above value is a good approximation provided that we only use image processing functions that preserve byte encoding and do not require internal conversion to Real32 or Real type. One example of such a function is ColorNegate. What if we want to call ImageConvolve? Then, not only the internal computations have to be performed using floating points, but also the resulting image is stored in the Real type (8 bytes per value). Our memory limit decreases to something like 0.1 GB, and the maximal width and resolution (in megapixels) become:

maxMemoryInBytes2 = 0.1 1024^3 ; maxWidth2 = Floor[Sqrt[maxMemoryInBytes2/3.0]] maxResolution =   Quantity[Round[(maxWidth2 * maxWidth2)/10^6], "MegaPixels"]

5982

36 megapixels

Thus, to do image processing on images larger than about 36 megapixels, we need out-of-core processing. Before we show how to do out-of-core processing on an image, let’s write some code that will allow us to visualize large images. In this post we will use an image from the NASA Visible Earth collection. It is a 222 megapixel JPEG file:

inputFile = FileNameJoin[{Directory[], "world_topo.jpg"}]; URLSave["http://eoimages.gsfc.nasa.gov/images/imagerecords/74000/\ 74518/world.topo.200412.3x21600x10800.jpg", inputFile]; dims =   Import[inputFile, "ImageSize"]

{21600, 10800}

Let’s start with inspecting the data and visualizing a thumbnailed version of our large image. The function LargeImageInspect allows us to specify the range of rows and columns of the image that we want to see. Notice that we use the “FirstRow” and “LastRow” options of the Import function that were introduced in Mathematica 9:

LargeImageInspect[file_, {r1_, r2_}, cspec_: All] :=    ImageTake[Import[file, "Image", "FirstRow" → r1, "LastRow" → r2],     All, cspec];

LargeImageInspect[inputFile, {3000, 3200}, {4000, 4600}]

Thumbnail version of the large image

Now, let’s create a function that shows a thumbnail of the large image. The function LargeImageThumbnail takes two arguments: a file path and optionally the width of a thumbnail image. It uses Import with the “FirstRow” and “LastRow” options to read a strip of a large image, re-sizes it according to the value of the second argument, and places the result in a list. At the end, ImageAssemble is called on that list to create a final thumbnail.

Creating a function that shows a thumbnail of the large image

LargeImageThumbnail[inputFile]

Thumbnail of the large image

To explore large images, we need one more utility—an out-of-core image viewer that will allow us to visualize actual pixel values stored in a file. Here I implement the function LargeImageViewer that shows the thumbnail of a large image in the top pane, and a part of the actual image data in the bottom pane. To see the implementation and all the code, please download this post as a Computable Document Format (CDF) file.

LargeImageViewer[inputFile]

Large image explorer

Now that we know how to explore large images, let’s take a closer look at the out-of-core image processing functions. Mathematica 9 offers three functions (ImageFileApply, ImageFileFilter, and ImageFileScan) that allow processing of large images in chunks and provide file-to-file workflow. All three functions work with TIFF, JPEG, and PNG files.

ImageFileApply is an out-of-core equivalent of the ImageApply function introduced in Mathematica 7. It applies a function specified as a first argument to the list of channel values for each pixel of the image.

outputFile = "world_topo_out.jpg"; ImageFileApply[  1 - # &, inputFile, outputFile]; LargeImageThumbnail[outputFile]

Color negative of the image

The above function call is equivalent to Export[outputFile, ColorNegate[Import[inputFile]]], however, ImageFileApply reads, processes, and writes successive blocks of image data so it requires significantly less memory than the in-core workflow. Here is the proof. The following code shows that ImageFileApply can be successfully executed when memory is constrained to 100 MB. Notice that in the call to ImageFileApply we used a non-default value of the “MaxBlockSize” option that defines the maximal size of image blocks in pixels. The default value of 8000000 is too large to run computations with the 100 MB memory constraint.

MemoryConstrained[  ImageFileApply[1 - # &, inputFile, outputFile,    "MaxBlockSize" → 100*dims[[2]]], 100*2^20]

world_topo_out.jpg

Trying to do the same thing with the in-core image processing functionality fails because it runs out of memory.

MemoryConstrained[Export[outputFile, ColorNegate[Import[inputFile]]],   100*2^20]

$Aborted

ImageFileApply‘s function parameter can be arbitrarily complex as long as it yields a number or a list of numbers. Here, we create a “blue snow” effect:

Thumbnail with a "blue snow" effect

ImageFileApply is a pixel operation and cannot be used for neighborhood operations such as filters. For this purpose, Mathematica 9 has a function called ImageFileFilter. We will use that function on an image from the ExampleData paclet:

inputFile = ExampleData[{"TestImage", "Stall"}, "FilePath"]; dims = ExampleData[{"TestImage", "Stall"}, "ImageSize"] LargeImageThumbnail[inputFile]

{4608, 3456}

Newly imported example image

A simple example of ImageFileFilter is to replace every value by the minimum in its range-r neighborhood (MinFilter). In the example below we have not specified the output file. In that case, both ImageFileApply and ImageFileFilter create an output file in the user directory with a name that is formed from the input file name by appending a time stamp.

outputFile = ImageFileFilter[Min, inputFile, 20] LargeImageThumbnail[outputFile]

Stall 04-18-2013 at 19.34.09.097.tif

Image with MinFilter applied

A popular Gaussian filtering is not much harder than the MinFilter. We create a kernel using the GaussianMatrix function, and then we sum up a product of that kernel with a neighborhood of each pixel, that is, Total[ker #,2]&:

r = 20; ker = GaussianMatrix[r]; LargeImageThumbnail[ImageFileFilter[Total[ker #, 2] &, inputFile, r]]

Image after applying Gaussian filtering

The last function that I want to discuss here is ImageFileScan. It is an out-of-core equivalent of ImageScan (also new in Mathematica 9) that evaluates a given function applied to each pixel of an image. Since ImageFileScan discards the results and does not create any file, it is useful in carrying out operations that have a “side effect.” For example, we can use this function to obtain different image statistics such as minimum and maximum values or image histograms:

Building a histogram for the image

Graphs for red, green, and blue values

All three out-of-core functions have one more interesting option called “ImageList”. It allows us to specify which frames in a multiframe TIFF file should be processed. In this example we color negate only the second and third frame of a multiframe image:

inputFile = FileNameJoin[{Directory[], "v43n1a2.tif"}]; URLSave["http://docmorph.nlm.nih.gov/docview/distrib/v43n1a2.tif",    inputFile]; frames = Import[inputFile, "ImageCount"] GraphicsRow@Import[inputFile]

4

Four example documents

GraphicsRow[  Import[ImageFileApply[1 - # &, inputFile, "ImageList" → {2, 3}]],   ImageSize → Small]

Documents after processing

The three functions that I have described here—ImageFileApply, ImageFileFilter, and ImageFileScan—are just the beginning of our work on out-of-core image processing. Ultimately, in future versions of Mathematica, the file-to-file workflow will be possible for all built-in image processing functions.

Download this post as a Computable Document Format (CDF) file.

Note: “FirstRow” and “LastRow” options of JPEG import have been replaced by the option “TakeRows” in Mathematica 10.

Posted in: Image Processing
Leave a Comment

4 Comments


y.matsuda

Where is the implementation of the function LargeImageViewer?

Posted by y.matsuda    May 23, 2013 at 9:20 pm
Jorma Kilpi

I found this blog very useful. Thanks! Question: How would you implement out-of-core ImagePartition[ ] ??

Posted by Jorma Kilpi    May 24, 2013 at 5:21 am
    Piotr Wendykier

    There is no build-in function for out-of-core ImagePartition yet. You can do something close to that by using Import with “FirstRow” and “LastRow” options (works only with TIFF, JPEG and PNG files). If you want to get tiles instead of stripes, then you can call ImagePartition for each stripe. In the last step you need to export each tile to a file.

    Posted by Piotr Wendykier    May 27, 2013 at 4:48 am
      Jorma Kilpi

      Thanks again! I managed to do something like that. What I needed was actually a virtual ImagePartition such that I can work with ij-part of it, i and j as abstract indexes indicating the part I want to work at the time. I do not want to make separate files for each subimage.

      Posted by Jorma Kilpi    May 27, 2013 at 8:32 am


Leave a comment

Loading...

Or continue as a guest (your comment will be held for moderation):