Steve on Image Processing

Contents

Dealing with "Really Big" Images: Block Processing

Hi! This is Brendan Hannigan, back again to continue our discussion on working with very large images in MATLAB. In the previous blog I discussed a couple of different ways to view and explore large images using the Image Processing Toolbox. Today I'll take the next logical step down the "large image workflow" path and we'll explore how to process images that are too large to load into memory.

"Sounds good. let's do this!"

Since the entire image cannot be loaded into memory at one time, we opted for an incremental, file-to-file solution. Basically we want to read a part of your "input" image into memory, process it in some way, and then write the results back to a new file, the "output" image. We continue to do this until the entire image has been processed, avoiding Out of Memory errors.

"This seems familiar to me..."

The flow of data described here is very similar to an existing IPT function, blkproc, which allows for block processing of images, but like most other IPT functions, it only supports in-memory processing. Originally, we considered expanding the scope of blkproc to support file-to-file workflows as well, but there were some syntactic and behavioral issues that I wasn't comfortable with and we would've likely introduced some backwards incompatibilities.

We instead opted for a completely new function in release R2009b. The result? blockproc! Very creative name right? blockproc has all of the capabilities of blkproc and a LOT more!

"Wait, didn't I read on CSSM that blockproc is way slower than blkproc ?"

Ahem... well... Ok fine, yes, it USED to be slower. The initial release of blockproc went out with what I will refer to simply as "performance growth opportunities" (sorry about that).

However, with the release of R2010b, all that has changed! blockproc performance was dramatically improved and is now comparable with blkproc for similar tasks.

Let's have a quick look at how it works!

% read input image
A = imread('peppers.png');
imshow(A);
% define block size and function to run on each block of data
block_size = [64 64];
my_function = @(block_struct) block_struct.data(:,:,[3 2 1]);
% call blockproc!
B = blockproc(A,block_size,my_function);
figure;
imshow(B);

"What just happened?"

What we've done here is swapped the red and blue color channels of our RGB peppers image with some simple indexing. blockproc did this one 64x64 block at a time and assembled the results into a single output image. Things that were "mostly red" became "mostly blue" and vice versa while the purple background remained mostly unchanged.

Basically all blockproc needs is an input image, a block size, and a function handle to run (similar to blkproc);

You may notice that I used an "anonymous function" to create the function handle my_function, but that's not necessary. You can also provide a function handle to a function that's either defined on your path or sitting in your current directory. The only requirement is that the function must accept a "block struct" as it's sole input argument (that we will pass to it, from inside of blockproc), and must return the processed data.

"Uhm... 'block struct'? That seems overly complicated"

The block struct is a MATLAB struct that contains the block's image data (in the .data field) as well as several other pieces of useful information. Most important among these is the .location field, which contains the location (in your input image) where the block came from. This .location information opens up a vast new world of potential uses for blockproc that we won't get into here.

Most of the time you use blockproc you'll probably just use the .data field, but for any operation that changes depending on which block of the image you are processing, the .location field is key, and you'll be thanking me then! You can check out the blockproc doc to learn more about what other information we package into the block struct.

Now for a more "real" example, applying a low-pass filter to an image. First, the conventional way:

% read the original photo into memory
origp = imread('cameraman.tif');
imshow(origp);
% create a Gaussian low-pass filter
h = fspecial('gaussian',5,2);
% compute the derived photo
derp1 = imfilter(origp,h);
imshow(derp1);

Now, with blockproc. We will re-use the original photo, origp, as well as the Gaussian filter.

% create a function handle that we will apply to each block
myFun = @(block_struct) imfilter(block_struct.data,h);
% setup block size
block_size = [64 64];
% compute the new derived photo
derp2 = blockproc(origp,block_size,myFun);
imshow(derp2);

"What's with those lines all over your result?"

What we see here are artifacts from block processing. What happened? Well, as our function, imfilter, processes its input it will require some padding values near the edges of the image. By default, imfilter will use "zero padding" to supply these "pixels" that lie beyond the boundary of the actual image data.

But remember, blockproc processes each block totally independently from its neighbors, so as we processed each block, imfilter was diligently padding each one with zeros, causing these black lines to appear when the results were assembled into our final output image.

"Ok, so how do I fix it?"

Luckily, we thought of that. blockproc has several optional parameters that we can specify to control all aspects of padding and block borders. One in particular is 'BorderSize', which lets us specify some "overlap" between blocks. Since our filter is size 5x5, we need a 2 pixel overlap between all of our blocks. Here's how we do that:

border_size = [2 2];
derp3 = blockproc(origp,block_size,myFun,'BorderSize',border_size);
imshow(derp3);

"Ok great, you fixed it. What does this have to do with large images?"

Right. Here's the "large data" hook: blockproc can take string filenames as its "input image", and can subsequently write the "output" to a file by specifying a new 'Destination' parameter. When you specify a new output 'Destination', blockproc will not return the result image to the MATLAB workspace.

% specify a string filename as our input image
origp_filename = 'cameraman.tif';
% specify a string filename as the 'Destination' of our output data
derp_filename = 'output.tif';
% don't request an output argument this time!
blockproc(origp_filename,block_size,myFun,...
    'BorderSize',border_size,'Destination',derp_filename);
imshow(derp_filename);

Now that the input and output are specified as files instead of workspace variables, you can run this command regardless of image size. Images are read, processed, and written incrementally, one block at a time.

Need to apply a filter on a 3 gigabyte image? No problem! Trying to segment vegetation from a few terabytes of satellite imagery? No problem!

"Hmm, ok yea that's cool. What's the catch?"

There's no catch! Ok, there's a small catch. blockproc only supports reading and writing to TIFF and JPEG2000 format files "natively".

"You're killing me with these file format restrictions! I don't use TIFF!"

Hey don't worry! We have you covered. I said that blockproc only supports TIFF and JPEG2000 "natively". What I meant by that is blockproc has "built-in" support for those file formats, but the function can "adapt" (hint,hint) to many other formats... which I will talk about next time.

Stay tuned!

Thanks, Brendan. -SE

Get the MATLAB code

Published with MATLAB® 7.12

Contents

I want to thank Steve for letting me guest blog about some recent work we did to help our large data users work with their imagery in MATLAB using the Image Processing Toolbox.

Several releases ago, we had been hearing from customers that they had "really big" images that they couldn't work with because they were too large to load into memory. How big are these images? "VERY big", they assured us. They couldn't even LOOK at them, much less operate on them in MATLAB! They could do nothing! They were very frustrated. Well, over the never several releases we worked hard to provide better tools to help these customers out.

Over the next few weeks I'll discuss a series of (relatively) new features that have taken these users from hitting Out of Memory errors at every turn, to being able to visually explore their large imagery in a smooth, interactive tool, to being able to process their images incrementally from disc, all while working with arbitrarily large images of arbitrary file format! A much improved workflow!

For this example we'll be working with a (relatively) large TIFF file that ships with the Mapping Toolbox, boston.tif. This file is only about 37 MB in memory, so not problematic in itself, but the tools I discuss here will scale to any size image.

tiff_filename = 'boston.tif';
imfinfo(tiff_filename)

ans = 

                     Filename: 'C:\MATLABs\R2011a\toolbox\map\mapdemos\boston.tif'
                  FileModDate: '04-Jun-2007 17:12:10'
                     FileSize: 38729900
                       Format: 'tif'
                FormatVersion: []
                        Width: 4481
                       Height: 2881
                     BitDepth: 24
                    ColorType: 'truecolor'
              FormatSignature: [77 77 0 42]
                    ByteOrder: 'big-endian'
               NewSubFileType: 0
                BitsPerSample: [8 8 8]
                  Compression: 'Uncompressed'
    PhotometricInterpretation: 'RGB'
                 StripOffsets: [12x1 double]
              SamplesPerPixel: 3
                 RowsPerStrip: 256
              StripByteCounts: [12x1 double]
                  XResolution: 300
                  YResolution: 300
               ResolutionUnit: 'Inch'
                     Colormap: []
          PlanarConfiguration: 'Chunky'
                    TileWidth: []
                   TileLength: []
                  TileOffsets: []
               TileByteCounts: []
                  Orientation: 1
                    FillOrder: 1
             GrayResponseUnit: 0.0100
               MaxSampleValue: [255 255 255]
               MinSampleValue: [1 1 1]
                 Thresholding: 1
                       Offset: 38729291
             ImageDescription: '"GeoEye"'
                     DateTime: '2007:02:23 21:46:13'
                    Copyright: '"(c) GeoEye"'
           ModelPixelScaleTag: [3x1 double]
             ModelTiepointTag: [6x1 double]
           GeoKeyDirectoryTag: [24x1 double]
            GeoAsciiParamsTag: 'State Plane Zone 2001 NAD = 83|'

"I can't even LOOK AT my image!"

Our first foray into large data viewing was simply producing an "overview" style image. If your image is too large to load into memory then it is not possible to view it using imshow with any of the standard syntaxes. Before committing the time and resources necessary work with a large image (moving it, copying it, converting to a different format, etc) users just want to do a "sanity" check on the data.

To help them accomplish this task we provided (in release R2008a) a very simple option to the imshow function, the 'Reduce' parameter. When set to true, imshow will subsample a large TIFF file and display a reduced version of that image. The subsampling ratio is based on your reported screen resolution.

This parameter is only valid for TIFF images. We started with the TIFF format in the beginning because MATLAB has excellent support for incremental reading and writing of TIFF images using the Tiff class. We recognize that most users will not be working exclusively with TIFF imagery and will want to view files in other formats as well.

Now... I don't want to scoop my subsequent blog posts too much, but let's just say we hope that we have sufficiently "adapted" (<== hint) to these customers' needs.

Let's view a reduced version of our image.

imshow(tiff_filename,'Reduce',true,'InitialMagnification','fit');

You'll notice that if you zoom in using the figure zoom tool that you will not see higher resolution imagery. Here's a closer look at Fenway Park.

% Simulate interactive zoom
xlim([44 390])
ylim([2638 2858])

It's clear here that you are actually viewing a reduced-resolution version of the original image. This 'Reduce' parameter is intended for users who just want to get a "quick & dirty" overview of a large TIFF image, and is not appropriate for seeing details in large images.

"That's cute, but I need to see my ACTUAL image, not some decimated imposter!"

Next we wanted to provide some "real", industrial strength large image viewing tool. To provide the best interactive viewing experience possible, we must create a multiresolution version of large images so that we can pull in pixels from the appropriate zoom-level as you pan and zoom through the data set.

Most users will be familiar with this behavior; it's very similar to how online mapping tools work. As you navigate around online maps, you will see tiles of data at discrete zoom levels (resolutions) being assembled to make up the current map view.

To provide this feature we delivered a new function (in release R2009a), rsetwrite, which allows you to create a reduced resolution data set (R-Set) from a large TIFF or NITF file. Creating the R-Set can take some time depending on the size of the image, but it is a one-time, up front cost.

To create an R-Set, the syntax is very straightforward. Let's build an R-Set from our large boston.tif image.

rset_filename = 'boston.rset';
rsetwrite(tiff_filename,rset_filename);

Now that you have your R-Set file, you can simply view it using the Image Tool as you would any other supported image file:

   imtool(rset_filename)

Here's a screenshot:

Using the Image Tool's pan & zoom tools you can now zoom into Fenway Park and see that, using the R-Set, we are now displaying the full resolution image data.

Here's another screenshot:

With R-Sets, you can zoom and pan around images of ANY SIZE, and the Image Tool will pull in only the data that it needs. Navigation of large imagery has never looked so good!

"Ok yea that's not bad, but I need to do more than just view the data!"

There are still many issues that I have not covered here such as, "I need to do more than just view my images", and "I'm not using TIFF or NITF files, what now?". These are excellent questions that I'll discuss next time. Stay Tuned!

Get the MATLAB code

Published with MATLAB® 7.12

Much of the information I posted in this blog years ago is still useful today. Image processing theory hasn't been completely overturned since then, and I'm still talking about MATLAB after all. For the benefit of readers who have joined the party more recently, I will occasionally recap the useful bits that I posted about five years ago.

In March 2006:

In April 2006:

Get the MATLAB code

Published with MATLAB® 7.12

Thinking about test images (13-Jan-2006, 16-Jul-2011, 19-Jul-2011, 19-Jul-2011) recently prompted me to take a look at the Image Processing Toolbox function checkerboard.

I = checkerboard;
imshow(I)

The implementation of checkerboard creates a four-by-four set of squares and then uses repmat, something like this:

   black = zeros(m);
   white = ones(m);
   tile = [black white; white black];
   I = repmat(tile,p,q);

(Then there's a postprocessing step to turn the white squares into gray squares on the right half of the image.)

This works just fine, but it occurred to me that it's a little inconvenient if you want an image with a size that's not a multiple of the square size. I had an idea for a different implementation of the checkerboard test pattern, one that is based on the function where n is an integer. This function bounces back and forth between -1 and 1.

n = 0:10;
(-1).^n

ans =

     1    -1     1    -1     1    -1     1    -1     1    -1     1

You can do this in two dimensions, also, with a function like .

[n1, n2] = ndgrid(0:4);
(-1).^(n1 + n2)

ans =

     1    -1     1    -1     1
    -1     1    -1     1    -1
     1    -1     1    -1     1
    -1     1    -1     1    -1
     1    -1     1    -1     1

To hurt your eyes a little, make a bigger version of this and display it as an image.

[n1, n2] = ndgrid(0:199);
I = (-1).^(n1 + n2);
imshow(I, [-1 1])

This is like a checkerboard but with extra tiny squares. We can make the squares bigger by using some scaling and rounding of the n1 and n2 values.

I = (-1).^(round(n1/20) + round(n2/20));
imshow(I, [-1 1])

Hmm, that resulted in half-tiles around the edges. Let's try floor instead of round.

I = (-1).^(floor(n1/20) + floor(n2/20));
imshow(I, [-1 1])

Now back to my observation above about the relationship between the size of the image and the number of squares. The implementation based on makes it very easy to specify independently the image height, image width, and square size. You can even turn the squares into rectangles.

image_height = 206;
image_width = 97;

square_height = 25;
square_width = 11;  % I know, I know; it's not a "square". It's a joke, OK?

[n1, n2] = ndgrid(1:image_height, 1:image_width);
I = (-1).^(floor(n1/square_height) + floor(n2/square_width));
% Make it a binary image of 0s and 1s instead of -1s and 1s.
BW = I == 1;
imshow(BW)

Last week's fun with the Jahne test pattern gave me the notion of making a checkboard image with spatially-varying frequency. The technique is the same - square and scale the exponent terms so that the frequency rises as you move away from the center of the image and reaches a maximum at the desired location. (And I'm going to switch back to using round.)

square_size = 20;
[n1, n2] = ndgrid(-100:100);
I = (-1).^(round(n1.^2/1000) + round(n2.^2/1000));
BW = I == 1;
imshow(BW)

I love MATLAB!

Get the MATLAB code

Published with MATLAB® 7.12

Today I want to take the test pattern I created last time and subject it to a variety of frequency-based filters.

In this post I'll be using a variety of frequency design and frequency-response visualization tools in MATLAB, Signal Processing Toolbox, and Image Processing Toolbox, including fir1, ftrans2, freqspace, fwind1, and freqz2. I won't be giving much explanation of these filter design techniques. Perhaps in the future I'll explore filter design here in more detail. For now, though, I really just want to have some fun with this test pattern.

[x,y] = meshgrid(-200:200);
r = hypot(x,y);
km = 0.8*pi;
rm = 200;
w = rm/10;
term1 = sin( (km * r.^2) / (2 * rm) );
term2 = 0.5*tanh((rm - r)/w) + 0.5;
g = (term1 .* term2 + 1) / 2;

imshow(g, 'XData', [-.8 .8], 'YData', [-.8 .8])
axis on
xlabel('Normalized frequency'), ylabel('Normalized frequency')
title('Test pattern')

Example 1: Design a one-dimensional lowpass filter and apply it in both directions.

b = fir1(30,0.3);
h = b' * b;
g1 = imfilter(g, h);

imshow(g1, 'XData', [-.8 .8], 'YData', [-.8 .8])
axis on
xlabel('Normalized frequency'), ylabel('Normalized frequency')

Example 2: Design a one-dimensional highpass filter and apply it in both directions.

b2 = fir1(30, 0.3, 'high');
h2 = b2' * b2;
g2 = imfilter(g, h2);

imshow(g2, [], 'XData', [-.8 .8], 'YData', [-.8 .8])
axis on
xlabel('Normalized frequency'), ylabel('Normalized frequency')

Example 3: Design a one-dimensional bandpass filter and apply it on both directions.

b3 = fir1(30, [0.3 0.4]);
h3 = b3' * b3;
g3 = imfilter(g, h3);

imshow(g3, [], 'XData', [-.8 .8], 'YData', [-.8 .8])
axis on
xlabel('Normalized frequency'), ylabel('Normalized frequency')

Example 4: Apply the one-dimensional bandpass filter from the previous step in only one direction.

g4 = imfilter(g, b3);

imshow(g4, [], 'XData', [-.8 .8], 'YData', [-.8 .8])
axis on
xlabel('Normalized frequency'), ylabel('Normalized frequency')

Example 5: Use ftrans2 to turn the one-dimensional bandpass filter from the previous example into a two-dimensional bandpass filter that is approximately circularly symmetric.

h5 = ftrans2(b3);
g5 = imfilter(g, h5);

imshow(g5, [], 'XData', [-.8 .8], 'YData', [-.8 .8])
axis on
xlabel('Normalized frequency'), ylabel('Normalized frequency')

Example 6: Use fwind1 to design a two-dimensional filter that has a wedge-shaped frequency response.

[f1, f2] = freqspace(101, 'meshgrid');
theta = cart2pol(f1, f2);
Hd = ((theta >= pi/6) & (theta < pi/3)) | ...
    ((theta >= -5*pi/6) & (theta < -2*pi/3));
mesh(f1, f2, double(Hd))
title('Ideal frequency response')

h6 = fwind1(Hd, hamming(31));
freqz2(h6)
title('Designed filter frequency response')

g6 = imfilter(g, h6);

imshow(g6, [], 'XData', [-.8 .8], 'YData', [-.8 .8])
axis on
xlabel('Normalized frequency'), ylabel('Normalized frequency')

Example 7: Compare the use of an averaging filter and a Gaussian filter. This example is inspired by Cris Luengo's comment from last week.

g7 = imfilter(g, fspecial('average', 9));
imshow(g7, [], 'XData', [-.8 .8], 'YData', [-.8 .8])
axis on
xlabel('Normalized frequency'), ylabel('Normalized frequency')
title('13-by-13 averaging filter')

You can see that some very high frequencies remain in the output. The Gaussian filter performs much better:

g9 = imfilter(g, fspecial('gaussian', 9, 2));
imshow(g9, [], 'XData', [-.8 .8], 'YData', [-.8 .8])
axis on
xlabel('Normalized frequency'), ylabel('Normalized frequency')
title('13-by-13 Gaussian filter')

OK, I'm going to give this test pattern a rest for now. Phew!

Get the MATLAB code

Published with MATLAB® 7.12

Earlier today I told you that I was feeling a little dense because I couldn't figure out the right parameters to use in the tanh term of this test pattern:

(This is equation 10.63 in Practical Handbook on Image Processing for Scientific Applications by Bernd Jahne.)

I'm grateful to reader Alex H for quickly enlightening me. He described as an approximation to a step function, where a is the location of the step and w is the width of the transition. That was very helpful.

I've also realized that I misinterpreted the meaning of . In the book this is described as the "maximum radius of the pattern," and I assumed this would be fixed as the distance from the center of the square image to one of its corners. But now I realize that the author intended for this to be an adjustable parameter. That is, one can set so that the maximum instantaneous frequency is reached closer to the center than at the image corners.

For example, I can set the parameters so that the maximum instantaneous frequency of is reached in the center of the image edges, and then the tapering function prevents aliasing artifacts from appearing as you move out to the corners. The book's figure 10.23 is based on maximum instantaneous frequency of (a period of 2.5 samples) reached at the edges, so I'll use that.

[x,y] = meshgrid(-200:200);
km = 0.8*pi;
rm = 200;
w = rm/10;
term1 = sin( (km * r.^2) / (2 * rm) );
term2 = 0.5*tanh((rm - r)/w) + 0.5;
g = term1 .* term2;
imshow(g,[])

Finally, a result that looks like the figure in the book!

Thanks again, Alex.

Get the MATLAB code

Published with MATLAB® 7.12

Last week I talked about the so-called Jähne test pattern and how I found a description of it in Practical Handbook on Image Processing for Scientific Applications, by Bernd Jähne, CRC Press, 1997. This equation was given in the book:

And here's a scan of the test pattern from the book:

I showed how to compute a test pattern based just on sine term. Here's the result:

You can see that the sine-based test pattern extends all the way to the image corners with no decrease in amplitude, whereas the book's test pattern tapers off towards the edges of the image. The book says that "multiplication by the tanh function results in a smooth transition of the ring patterns at the outer edge avoid aliasing problems." That made sense to me, but I couldn't get my generated image to look like what's in the book.

Here's what I tried:

km = pi;
[x,y] = meshgrid(-200:200);
r = hypot(x,y);
rm = max(r(:));
term1 = sin( (km * r.^2) / (2 * rm) );

w = 300;
term2 = 0.5 * tanh((rm - r)/w) + 0.5;
imshow(term2, [])
title('tanh term')

g = term1 .* term2;
imshow(g,[])

That's not much of a taper. I couldn't get my test pattern to look like the figure in the book no matter what value I chose for w. I also noticed that the 2nd term goes down to only 0.5, whereas I would expect a tapering function to range between 0.0 and 1.0. But even when I adjusted the scale and offset of the 2nd term so that it would range between 0.0 and 1.0, I couldn't anything that looked like the book's figure.

So I took a closer look at the shape of .

x = -1:.01:1;
plot(x,tanh(1-abs(x)))

That didn't make sense to me for a taper function. Based on the figure in the book, I expected the taper function to be flatter in the center and to trend more smoothly toward 0 at the edges.

Dear reader, I'm feeling a little dense because I never really figured out this puzzle. If I'm missing something here, please feel free to tell me in the comments. [UPDATE: Reader Alex H. helped me figure out the problem. See this follow-up post.]

I thought about constructing another taper function (maybe the Tukey window?), but I'm not sure how necessary this really is. For example, if you're concerned about aliasing effects because of the high frequency pattern at the image corners, you could just lower a bit:

km = 0.5*pi;
g = (sin( (km * r.^2) / (2 * rm) ) + 1)/2;
imshow(g)

I've got one more post planned about this image (fun with bandpass filtering), and then I plan to show you an interesting way to generate a checkerboard image (chessboard image?) based on the function .

Get the MATLAB code

Published with MATLAB® 7.12

Earlier this week I posted a summary of posts from the beginning of this blog, including my January 13, 2006 post on generating test images. Today I want to revisit one of the examples from that old post.

In my 2006 post I showed this test image, consisting of concentric rings with increasing frequency as you move away from the center of the image:

Since then I have seen this type of image referred to as the Jähne test pattern. This week I decided to go looking for the source of this reference, and I found a very similar image described in Practical Handbook on Image Processing for Scientific and Technical Applications, by Bernd Jähne, CRC Press, 1997. Here's a scan of Figure 10.23a, page 348, from that book:

The test image is constructed using equation (10.63):

where is the vector offset from the image center, is therefore the distance from the image center, and is the maximum value of in the image.

Today I'm going to ignore the second term (tanh) and just focus on the first term, and I'm going to substitute for . The instantaneous frequency of the sinusoid is given by:

The maximum instantaneous frequency, which occurs when , is therefore . For a discrete signal defined on a grid with unit spacing, the maximum representable frequency is , which corresponds to a period of 2 units.

OK, let's try it.

km = pi;
[x,y] = meshgrid(-200:200);
r = hypot(x,y);
rm = max(r(:));

g = sin( (km * r.^2) / (2 * rm) );
imshow(g,[])

I used the autoscaling syntax of imshow because the values of g range from -1 to 1. We can scale and shift g to get it into the range [0,1]. Then we don't need the second input to imshow, and the image will write out correctly if we use imwrite.

g = (g + 1)/2;
imshow(g)

Just for fun, let's push the maximum instantaneous frequency way past reasonable and look at the resulting aliasing artifacts.

km = 4*pi;
g = sin( (km * r.^2) / (2 * rm) );
g = (g + 1) / 2;
imshow(g)

I spent really much too long playing with this image yesterday, so I'm going to subject you to one or two more posts about it. For one thing, having the instantaneous frequency vary directly with distance from the image means we can explore interesting effects using bandpass filters. For another, I want to come back to the tanh term in the equation from Jähne's book. I'm a bit puzzled by it.

Get the MATLAB code

Published with MATLAB® 7.12