# Steve on Image Processing

## Generating code from Image Processing Toolbox functions

I'd like to welcome guest blogger Vignesh Krishnan for today's post. Vignesh is a developer on the Image Processing Toolbox team. -Steve

Let's try C code generation from Image Processing Toolbox™ functions. I am going to use an example to fill borders of a binary image from Steve Eddins' post.

In this post, the goal is to fill holes in an image and generate code using the Image Processing Toolbox and MATLAB Coder™. First, I create a MATLAB function with code adapted from the blog post mentioned earlier and name it fillborderholes.m

type fillborderholes.m

function bw_filled = fillborderholes(bw)
%#codegen

% Fill against the left and the top border. The following call to padarray
% adds a column of white pixels on the left and a row of white pixels on
% the top.
bw_a_filled = imfill(bw_a,'holes');

% Now fill against the top and the right border.
bw_b_filled = imfill(bw_b,'holes');

% Next, fill against the right and bottom borders.
bw_c_filled = imfill(bw_c,'holes');

% Finally, fill against the bottom and left borders. (I wouldn't expect
% this to have any effect on this particular image.)
bw_d_filled = imfill(bw_d,'holes');

% The last step is then to "logical OR" all these images together.


Next, I have to create an example input to be used while generating code. C is a static typed language hence an example input is required to constrain the generated code to a particular input size and datatype.

url = 'http://blogs.mathworks.com/images/steve/2013/eye-heart-matlab.png';
imshow(bw)


Now I am ready to generate C code from my function using the codegen command. The '-args' option provides the example input. The '-report' option creates a compilation report which provides links to the generated code. I setup the configuration parameters to turn off run-time checks to make the code run faster and improve its readability. For more details on run-time checks look here

cfg = coder.config('mex');
cfg.IntegrityChecks = false;
cfg.ResponsivenessChecks = false;
codegen fillborderholes -args {bw} -config cfg -report

Code generation successful: To view the report, open('/Volumes/inside_files_dev/ltc/imagemap/ipblog/published/2013/codegen/mex/fillborderholes/html/index.html').


Here's what the generated report looks like:

The default coder configuration setting creates a MEX file in the current directory named fillborderholes_mex. Ok! Let's run the executable

bw_filled_cg = fillborderholes_mex(bw);
imshow(bw_filled_cg);


By default, the code generated for the function is in the file, fillborderholes.c in the folder codegen/mex/fillborderholes/. Alternatively, you can click on the C-code tab in the report and look at the generated code. It includes function calls to imreconstruct which is the morphological reconstruction algorithm used by imfill. The code uses a precompiled platform-specific shared library to implement the algorithm.

Let me focus on a few excerpts of the generated code. Notice how the comments in my MATLAB code appeared in my generated code. The comments help correlate the generated code with the source code. The top-level function, fillborderholes.m for which we generated code takes an image, bw as input and returns bw_filled as output. The input image used in this example at compilation time is of boolean datatype and size 445 x 501 i.e. 222945 elements. As a result, the function declaration is as follows

void fillborderholes(const boolean_T bw[222945], emxArray_boolean_T *bw_filled)

The function expects the first input to be a boolean array of size 222945. A static array is used in the C interface because the input was defined to be fixed size during compilation. The output is dynamically allocated which is represented as a structure type called emxArray specialized for the particular datatype. For more details on emxArray look here

The first operation in fillborderholes.m is to pad the input image using padarray and then call imfill on it. As shown below, this translates to the function call ConstantPad() which pre-pads the input image, bw to create a padded image, bw_a. The call to imfill() creates a filled image, bw_b which is then copied over to bw_a_filled, since bw_b is reused later in the generated code.

The rest of the generated code is a series of calls to ConstantPad() to pad different borders of the image followed by calls to imfill() to fill the padded image.

In the last section, the generated code creates two for-loops to OR the filled results from the intermediate steps and store the final filled result in bw_filled.

As of R2013a, the Image Processing Toolbox using MATLAB Coder can generate code from the following functions: conndef, imcomplement, imfill, imhmax, imhmin, imreconstruct, imregionalmax, imregionalmin, iptcheckconn, and padarray. Which functions would you like to be able to generate code?

Get the MATLAB code

Published with MATLAB® R2013b

## Watershed transform question from tech support

A support call came in this week from a customer trying to use watershed to segment this image:

The complaint was that calling watershed did not produce a good segmentation.

Today I want to show how to use watershed to segment this image. Along the way I'll explain the difference between the watershed transform and watershed segmentation.

First, let's just try calling watershed and see what happens.

url = 'http://blogs.mathworks.com/images/steve/2013/blobs.png';
L = watershed(bw);
Lrgb = label2rgb(L);
imshow(Lrgb)


When I saw that result, I was puzzled at first. Then I realized what was going on. Let me use imfuse to show these two images together, zooming in on one particular blob.

imshow(imfuse(bw,Lrgb))
axis([10 175 15 155])


You see, the watershed transform always gives you one watershed region for every local minimum (or regional minimum) in the image. These little black "noise spots" are local minima, and so there's a watershed region around each one.

Even if we fill these holes, though, just using the watershed transform by itself is never going to produce the segmentation that the customer was seeking. That brings me to my point about the distinction between watershed segmentation and the watershed transform. Watershed segmentation refers to a family of algorithms that are based on the watershed transform. Except for very specific cases, the watershed transform isn't a full segmentation method on its own.

Some years ago, I wrote a MathWorks newsletter article called The Watershed Transform: Strategies for Image Segmentation. It's worth reviewing in order to brush up on the basics. An central concept from the article is this:

The key behind using the watershed transform for segmentation is this: Change your image into another image whose catchment basins are the objects you want to identify.

For an image such as this, consisting of roughly circular, touching blobs. the distance transform can be useful for producing an image whose "catchment basins are the objects you want to identify."

Before go to the distance transform, though, let's clean up the noise a bit. The function bwareaopen can be used to remove very small dots. It removes them in the foreground, though, so we complement the image before and after calling bwareaopen.

bw2 = ~bwareaopen(~bw, 10);
imshow(bw2)

D = -bwdist(~bw);
imshow(D,[])


Now we're starting to get somewhere. Next, compute the watershed transform of D.

Ld = watershed(D);
imshow(label2rgb(Ld))


The watershed ridge lines, in white, correspond to Ld == 0. Let's use these ridge lines to segment the binary image by changing the corresponding pixels into background.

bw2 = bw;
bw2(Ld == 0) = 0;
imshow(bw2)


The "raw" watershed transform is known for its tendency to "oversegment" an image. The reason is something I mentioned above: each local minimum, no matter how small, becomes a catchment basin. A common trick, then, in watershed-based segmentation methods is to filter out tiny local minima using imextendedmin and then modify the distance transform so that no minima occur at the filtered-out locations. This is called "minima imposition" and is implemented via the function imimposemin.

The following call to imextendedmin should ideally just produce small spots that are roughly in the middle of the cells to be segmented. I'll use imshowpair to superimpose the mask on the original image.

mask = imextendedmin(D,2);


Home stretch, now. Modify the distance transform so it only has minima at the desired locations, and then repeat the watershed steps above.

D2 = imimposemin(D,mask);
Ld2 = watershed(D2);
bw3 = bw;
bw3(Ld2 == 0) = 0;
imshow(bw3)


That's it! Read my newsletter article for other thoughts on image segmentation using the watershed transformation.

Get the MATLAB code

Published with MATLAB® R2013b

## Chess and a little text file manipulation

Here's an image of a chess position:

And that's about as close to image processing as today's blog post will come. Because this post is really about text processing.

It seems like a lot of computational tasks in engineering and science involve manipulating data in text files. This weekend I had such a task, although I must admit that it had nothing to do with engineering or science. Even so, I thought the task would be a good illustration of some basic text processing techniques.

I have a text file, tactics.pgn, that is a database of chess tactics puzzles. Here are lines 2,759 through 2,784 from the file.

[Event "?"]
[Site "?"]
[Date "????.??.??"]
[Round "?"]
[White "Diagram 211"]
[Black ""]
[Result "*"]
[EventDate "2006.05.20"]
[FEN "2b1R3/ppk4p/8/2q2p2/2Br2n1/2QP2N1/P4PPP/6K1 w - - 0 1"]
[SetUp "1"]
[SourceDate "2011.02.22"]
1.Rxc8+ Kxc8 2.Be6+ Kd8 3.Qxc5 *
[Event "?"]
[Site "?"]
[Date "????.??.??"]
[Round "?"]
[White "Diagram 212"]
[Black ""]
[Result "*"]
[FEN "3rr1k1/p1p2ppp/Q2b1q2/8/3Np3/4P3/PP3PPP/R1B2RK1 b - - 0 1"]
[SetUp "1"]
[SourceDate "2011.02.22"]
1...Bxh2+ 2.Kxh2 Qxa6 *

These lines store two chess positions. The first is the position I showed above, with White to move. Can you figure out how White wins by capturing the Bishop on c8 with his rook? (If you're interested in how the position is encoded in the text, see the Wikipedia article on Forsyth-Edwards Notation, or FEN.)

I wanted to shuffle the positions in this file randomly. None of the chess software programs that I have can do this, so I decided to tackle it with MATLAB. (Fun thing to do on a Saturday morning, right?)

OK, so here's the basic procedure:

1. Read the entire file into MATLAB.
2. Split the data into chunks, one chunk for each position.
3. Rearrange the chunks randomly.
4. Write out the rearranged chunks to a new text file.

There are a lot of different ways to do this. Here's what I came up with.

First, read in the entire file. The function fileread is just the ticket.

characters = fileread('tactics.pgn');
size(characters)

ans =

1      106394



You can see that there are about 106,000 characters in the file. Let's split the data into lines using strsplit.

lines = strsplit(characters,'\n')';
size(lines)

ans =

4838           1



There are about 4,800 lines of text. But how many positions are there? I'm going to find the starting line of each position by searching for the string "[Event " at the beginning a line. It's time for regexp.

idx = regexp(lines,'^[Event ');
idx(1:15)

ans =

[1]
[]
[]
[]
[]
[]
[]
[]
[]
[]
[]
[1]
[]
[]
[]



This shows us that this string is found twice in the first 15 lines of the file. idx is a cell array, so I'll use cellfun and find to identify all the lines that contain the matching string. Each of these lines is the start of an entry for one chess position.

first_lines = find(~cellfun(@isempty,idx));
first_lines(1:3)

ans =

1
12
23



So there are positions starting on lines 1, 12, and 23.

Next, I'll make a cell array such that that each cell contains all the lines for one position. To make the for-loop work, I'm going to add an "extra" value to the first_lines vector that points to a nonexistent line just past the of the file.

first_lines(end+1) = length(lines) + 1;
for k = 1:length(first_lines)-1
positions{k} = lines(first_lines(k):first_lines(k+1)-1);
end


Let's take a look at what we have now.

size(positions)

ans =

1   421



There are 421 positions in the file. For example:

positions{205}

ans =

'[Event "?"]'
'[Site "?"]'
'[Date "????.??.??"]'
'[Round "?"]'
'[White "Diagram 211"]'
'[Black "Bain"]'
'[Result "*"]'
'[EventDate "2006.05.20"]'
'[FEN "2b1R3/ppk4p/8/2q2p2/2Br2n1/2QP2N1/P4PPP/6K1 w - - 0 1"]'
'[SetUp "1"]'
'[SourceDate "2011.02.22"]'
'1.Rxc8+ Kxc8 2.Be6+ Kd8 3.Qxc5 *'



Again, this is the position shown at top of this post.

Getting near the end, now. It's time to rearrange the positions. Before I do that, though, I'll shuffle the random number generator. I only do this so that if I repeat these steps in a new MATLAB session, I'll be sure to get a different result. After shuffling the random number generator using rng, a quick call to randperm randomly rearranges the positions.

rng shuffle
shuffled_positions = positions(randperm(length(positions)));


We've arrived at the last step: writing out the shuffled positions to a new file.

fid = fopen('shuffled_tactics.pgn','w');
for k = 1:length(shuffled_positions)
position = shuffled_positions{k};
for p = 1:length(position)
fprintf(fid,'%s\n',position{p});
end
end
fclose(fid);


And that's it!

Reading text, manipulating it in some useful way, and writing the results back out -- a common computing task accomplished using several basic MATLAB functions.

Get the MATLAB code

Published with MATLAB® R2013b

## timeit makes it into MATLAB

This is my first blog post with "Published with MATLAB R2013b" at the bottom. The latest MATLAB release shipped earlier in September. And, for the first time in a while, a function that I wrote has made it into MATLAB.

Back in 2008, I spent some time trying to incorporate performance benchmarking lessons I learned from Cleve Moler and Bill McKeeman into a general-purpose function for accurately measuring the "typical" running time of MATLAB functions or expressions. The result was timeit, which I submitted to the File Exchange.

Well, I'm happy to say that timeit has made it into MATLAB in the R2013b release! The MATLAB development team took the code from the File Exchange, made some minor refinements, did some additional testing, wrote some nice doc, and got it into the release. (Thanks, everyone!)

Here are a couple of examples. The first one is just MATLAB, but the second one requires Image Processing Toolbox. (Note that it's helpful to know a little about anonymous functions in order to use timeit.)

How much time does it take to compute sum(A.' .* B, 1), where A is 12000-by-400 and B is 400-by-12000?

A = rand(12000,400);
B = rand(400,12000);
f = @() sum(A.' .* B, 1);
timeit(f)

ans =

0.0397



How much time does it take to dilate the text.png image with a 25-by-25 all-ones structuring element?

bw = imread('text.png');
se = strel(ones(25,25));
g = @() imdilate(bw,se);
timeit(g)

ans =

7.7779e-04



Loren also wrote about timeit recently in her Art of MATLAB blog.

There are other goodies in R2013b that I'll write about soon.

Get the MATLAB code

Published with MATLAB® R2013b

## Defining and filling holes on the border of an image

url = 'http://blogs.mathworks.com/images/steve/2013/eye-heart-matlab.png';
imshow(bw)


Alex tried imfill but it didn't give him exactly what he was looking for.

bw2 = imfill(bw,'holes');
imshow(bw2)
title('Filling holes (first attempt)')


Alex asked how he could fill the objects that are touching the border. At this point, MathWorker Sean de Wolski jumped in and pointed out the problem wasn't well defined; there are multiple interpretations and possible solutions.

It might not be obvious from Alex's image that this is the case, so let me illustrate with a different example. This image has one object, and it touches the border in a way that divides the background into two regions. Topologically, this object is the same as each of the five "unfilled" objects in Alex's image.

bw_eye = imdilate(eye(200),ones(5,5));
imshow(bw_eye)


imfill doesn't do anything to this image:

imshow(imfill(bw_eye,'holes'))
title('Output of imfill')


What does it mean to "fill" this object? Would you expect the upper-right portion of the background to be filled, or the lower-right? This is the point Sean was making.

Well, I stared at Alex's image for a while, trying to come up with a formulation of the problem that would be better defined and that would do what Alex expects.

I think I've hit upon a reasonable formulation. It's not the only one, however, and I'd be interested in what ideas other readers have about this.

Here's my idea: Fill objects against two adjacent borders at a time. By "fill against a border", I mean pad that border with white pixels, do the fill, and then remove the padding.

Let's go back to Alex's image and "fill against" the left and the top border. The following call to padarray adds a column of white pixels on the left and a row of white pixels on the top.

bw_a = padarray(bw,[1 1],1,'pre');
bw_a_filled = imfill(bw_a,'holes');
bw_a_filled = bw_a_filled(2:end,2:end);
imshow(bw_a_filled)


Now fill against the top and the right border.

bw_b = padarray(padarray(bw,[1 0],1,'pre'),[0 1],1,'post');
bw_b_filled = imfill(bw_b,'holes');
bw_b_filled = bw_b_filled(2:end,1:end-1);
imshow(bw_b_filled);


Next, fill against the right and bottom borders.

bw_c = padarray(bw,[1 1],1,'post');
bw_c_filled = imfill(bw_c,'holes');
bw_c_filled = bw_c_filled(1:end-1,1:end-1);
imshow(bw_c_filled)


Finally, fill against the bottom and left borders. (I wouldn't expect this to have any effect on this particular image.)

bw_d = padarray(padarray(bw,[1 0],1,'post'),[0 1],1,'pre');
bw_d_filled = imfill(bw_d,'holes');
bw_d_filled = bw_d_filled(1:end-1,2:end);
imshow(bw_d_filled)


The last step is then to "logical OR" all these images together.

bw_filled = bw_a_filled | bw_b_filled | bw_c_filled | bw_d_filled;
imshow(bw_filled)
title('Final filled result')


Alex, I hope you find this useful.

Get the MATLAB code

Published with MATLAB® R2013a

## Introduction to spatial referencing

Today's blog post was written by Alex Taylor, a developer on the Image Processing Toolbox team. Thanks, Alex!

In MATLAB, an image is typically represented as a numeric matrix. There are some naturally geometric image processing operations where a user may need to care about where an image is located in a coordinate system. The term "spatial referencing" refers to the process of defining where an image is located in some coordinate system. This blog post will explore spatial referencing and the various ways to accomplish the spatial referencing of image data in MATLAB and the Image Processing Toolbox.

### Contents

#### Background: The Intrinsic Coordinate System

There is a consistent, default coordinate system used in MATLAB and the Image Processing Toolbox. This coordinate system is used by default for all naturally geometric operations such as image display, and geometric transformations. You can see the default coordinate system by displaying a small image and observing where the pixels are positioned within the coordinate system of the axes.

I = reshape(9:-1:1,3,3);
figure;
imagesc(I);
grid on;
colormap gray


In the default coordinate system, the center of each pixel falls on an integral value that corresponds to the one-based (column,row) subscripts of the image. For example, the center of the first pixel in I is located at coordinate (1.0,1.0).

In the Image Processing Toolbox, the default coordinate system is referred to as the "intrinsic coordinate system." The intrinsic coordinate system is a continuous extension of the discrete column,row subscripts of an image grid. This relationship between column,row matrix subscripting is where the term "intrinsic coordinates" comes from. The coordinate system is intrinsic to the one based column,row indexing of MATLAB.

#### Spatial referencing objects: imref2d and imref3d

The spatial referencing objects imref2d and imref3d were introduced in the R2013a release of the Image Processing Toolbox to make it easier to work with 2-D and 3-D images in a non-default or "world" coordinate system. This blog post will focus on the 2-D case and imref2d. The following example defines a world coordinate system for the simple image I where the X limits of the image are [3 6] and the Y limits the image are [6 9]. The X and Y world limits can be defined in any consistent unit system. In this example, assume the world limits are defined in units of meters.

RI = imref2d(size(I));
RI.XWorldLimits = [3 6];
RI.YWorldLimits = [6 9];


The equivalent spatial referencing object can be created in one line with a three-argument syntax.

RI = imref2d(size(I),[3 6],[6 9]);


Together, an image, I, and its associated spatial referencing object, RI, describe a spatially referenced image. The imshow function accepts a spatially referenced image.

figure, imshow(I,RI,[],'InitialMagnification','fit');


#### Properties of imref2d

Spatial referencing objects have properties that describe various aspects of spatial referencing. The PixelExtentInWorldX and PixelExtentInWorldY properties define how large each pixel is in the world coordinate system along the X and Y dimensions. Because the XWorldLimits and YWorldLimits properties were defined in meters, the PixelExtentInWorldX and PixelExtentInWorldY properties are in units of meters.

RI.PixelExtentInWorldX
RI.PixelExtentInWorldY

ans =

1

ans =

1



The ImageExtentInWorldX and ImageExtentInWorldY properties define the size of the image grid along the X and Y dimensions in the units of the world coordinate system.

RI.ImageExtentInWorldX
RI.ImageExtentInWorldY

ans =

3

ans =

3



The ImageSize property defines the number of rows and columns in the associated image matrix.

RI.ImageSize

ans =

3     3



#### Converting between coordinate systems using spatial referencing objects

The spatial referencing objects have a set of related functions that can be used to convert from world coordinates to the default intrinsic coordinate system and to discrete pixel subscripts. For example, the worldToIntrinsic function can be used to find the intrinsic coordinates with the image I that correspond to the world (X,Y) coordinate (4.6,7.3).

[intrinsicX,intrinsicY] = worldToIntrinsic(RI,4.6,7.3)

intrinsicX =

2.1000

intrinsicY =

1.8000



The intrinsic coordinates that coorespond to this world coordinate pair fall at intrinsic (X,Y) location (2.1,1.8). This means that the value we are interested in is located between the first and second row and between that second and third column of the image I. The worldToSubscript function can be used to obtain the nearest integral row, column subscript into I at a given world location.

[r,c] = worldToSubscript(RI,4.6,7.3)

r =

2

c =

2



The output (row,column) location can be used to index into I to find the value of I at world (X,Y) location (4.6,7.3).

I(r,c)

ans =

5



The use of worldToSubscript to find the closest integrally valued subscript to a given world coordinate location is equivalent to performing nearest neighbor interpolation within the spatially referenced image grid.

interp2(I,intrinsicX,intrinsicY,'nearest')

ans =

5



You can also use bilinear interpolation to estimate the value of I at intrinsic coordinates that fall between integral row,column indices.

interp2(I,intrinsicX,intrinsicY,'bilinear')

ans =

4.9000



#### imshowpair and spatial referencing

The image display function imshowpair can be used to display overlays of image pairs and is a good example of the importance of spatial referencing in naturally geometric workflows.

Start with an MRI image of my knee after a running injury.

knee = dicomread('knee1.dcm');


Resize the knee image by a factor of 0.5.

smallKnee = imresize(knee,0.5);


Display a falsecolor overlay of each image. Note that because no spatial referencing information is specified to imshowpair, imshowpair assumes each image is in the intrinsic coordinate system. This causes each image to be aligned at the upper left corner (0.5,0.5).

figure, imshowpair(knee,smallKnee);


There is a one input argument syntax of imref2d that allows for easy construction of the default spatial referencing that corresponds to an image of a given size. We can pass these spatial referencing objects to imshowpair to see the same visualization of each image overlayed in the intrinsic coordinate system.

Rknee    = imref2d(size(knee));
RsmallKnee = imref2d(size(smallKnee));
figure, imshowpair(knee,Rknee,smallKnee,RsmallKnee);


In this case, the physical size of the knee in the world is the same, but one of the images is sampled at half the resolution of the other. Adjusting the world limits of the resized image has the effect of simulating a case in which two images of different resolution contain a view of the same scene in the world.

RsmallKnee.XWorldLimits = Rknee.XWorldLimits;
RsmallKnee.YWorldLimits = Rknee.YWorldLimits;
figure, imshowpair(knee,Rknee,smallKnee,RsmallKnee);


Applying a shift to the world limits of one of the images simulates image misalignment due to translation.

RsmallKnee.XWorldLimits = RsmallKnee.XWorldLimits+50;
figure, imshowpair(knee,Rknee,smallKnee,RsmallKnee);


#### The future direction of spatial referencing in the Image Processing Toolbox

Older Image Processing Toolbox functions such as imshow and imtransform use the the MATLAB graphics 'XData' and 'YData' referencing system. In this system, the two element vector XData described the centers of the first and last pixel along the X extent of the image, and the YData described the centers of the first and last pixel along the Y extent of the image. See the following doc link for details: The 'XData' and 'YData' properties of the HG image object.

Newer Image Processing processing functions such as imshowpair, imregister, imregtform, and imwarp use spatial referencing objects instead of 'XData' and 'YData'. Moving forward, we will de-emphasize the use of 'XData' and 'YData' referencing in the Image Processing Toolbox. We will expand the use of spatial refeferencing objects in new functions and find ways to integrate spatial referencing into existing functions that are part of naturally geometric workflows.

Are there any Image Processing Toolbox functions in which you would like to see the introduction of spatial referencing? If so, please let us know, we would appreciate your feedback.

Get the MATLAB code

Published with MATLAB® R2013a

## TIFF, BigTIFF, and blockproc

I'd like to welcome back guest blogger Ashish Uthama for today's post. Ashish, a developer on the Image Processing Toolbox team, posted here previously about the mystery of the failing TIFF append loop.

### Contents

#### Brief introduction to the TIFF format

A TIFF file is, well, a Tagged Image File Format (TIFF). As the name implies, the file itself consists of a whole bunch of 'tags' and corresponding values. For example, all TIFF files are expected to have the ImageWidth and ImageLength tags, whose values are the width and length (height) of the image stored within.

A core feature of this format lies in how image data is organized in a physical file. TIFF allows for the image data to be stored in chunks as either 'Tiles' or 'Strips'. You could think of a strip as a tile whose width is the same as the image. The format stores offsets for each of these strips or tiles in an image header. This allows for efficient random access of any chunk in the image. The original TIFF format specification calls for use of 32 bit file offset values. I guess you can imagine where I am going with that bit of information. 32 bits implies a maximum value of 2^32 for any offset. Note that these are offsets from the very beginning of the file. In effect, the 32 bit requirement limits the size of the largest possible TIFF file to be under 4 gigabytes. This was fine for a long time (The earliest TIFF release was in 1986), but we now have more data producers hitting this limit.

#### Enter BigTIFF

The libtiff community, a loosely organized collection of volunteer software developers, maintains and develops the LibTIFF library. LibTIFF implements routines to manipulate TIFF files. Since version 4.0.0 of LibTIFF, they have included a variant of the TIFF format which uses 64 bits for the offsets instead of 32 bits. As you can imagine, this allows for really big TIFF files. This variant is appropriately called BigTIFF . Do note that this underlying change causes BigTIFF to an incompatible variant to the 32 bit offset TIFF format (now called the classic format). The first stable version of libtiff with BigTIFF support came out just in time for us to be able to qualify and test it to go out with MATLAB R2012b.

#### BigTIFF support in MATLAB

As I mentioned in the introduction, TIFF files have tags and values. This approach gives the format a lot of flexibility in the type of data it can hold, and yes, with great flexibility comes great complexity.

Reading BigTIFF file in MATLAB is no different from reading classic (32bit offset) version TIFF files. imread and all its TIFF related options work with BigTIFF. With BigTIFF files, the most relevant option is the 'PixelRegion' parameter. If a BigTIFF image file is too large to fit in system memory, you could use this parameter to load a smaller part of the image. imwrite can easily write out classic TIFF files. However, writing out BigTIFF files is a bit more involved. The key assumption with using imwrite is that the entire image data is available as a matrix in memory, something not necessarily true with BigTIFF. Tiff images can be complex and there are hundreds of possible combinations of the various tag values. In addition, the format enforces rules (see Tables 2, 3, 4 and 5) for legal combinations of tag values.

MATLAB provides a gateway to the LibTIFF library routines through the Tiff class. This interface supports over 60 tags, allowing MATLAB users to work with a wide range of TIFF files. Moreover, this interface also supports working with BigTIFF files. Here is a code snippet which uses the Tiff class to write a (albeit small) BigTIFF file:

imdata = imread('example.tif');
% 'w'  will create a classic TIFF file
% 'w8' will create a BigTIFF file
% This option is the only differentiator when writing to these two formats.
bt     = Tiff('myfile.tif','w8');


tags.ImageLength         = size(imdata,1);
tags.ImageWidth          = size(imdata,2);
tags.Photometric         = Tiff.Photometric.RGB;
tags.BitsPerSample       = 8;
tags.SamplesPerPixel     = size(imdata,3);
tags.TileWidth           = 128;
tags.TileLength          = 128;
tags.Compression         = Tiff.Compression.JPEG;
tags.PlanarConfiguration = Tiff.PlanarConfiguration.Chunky;
tags.Software            = 'MATLAB';

setTag(bt, tags);
write(bt,  imdata);
close(bt);


#### Using BigTIFF with BLOCKPROC (Image Processing Toolbox)

type bigTiffWriter.m

classdef bigTiffWriter < ImageAdapter
%BIGTIFFWRITER - A basic image adapter to write Big TIFF files.
%
% object for use with BLOCKPROC.
% - Tile dimensions must be multiples of 16.
% - Only uint8, RGB input image data is supported.
%
% Based on "Working with Data in Unsupported Formats"
% http://www.mathworks.com/help/toolbox/images/f7-12726.html#bse_q4y-1
%
% Example:
%    %%
%    % Set file names and obtain size information from input file.
%    inFile        = 'example.tif';
%    inFileInfo    = imfinfo(inFile);
%    outFile       = 'out.tif';
%    %%
%    % Create an output TIFF file with tile size of 128x128
%    tileSize      = [128, 128]; % has to be a multiple of 16.
%    outFileWriter = bigTiffWriter(outFile, inFileInfo(1).Height, inFileInfo(2).Width, tileSize(1), tileSize(2));
%    %%
%    % Now call blockproc to rearrange the color channels.
%    blockproc(inFile, tileSize, @(b) flipdim(b.data,3), 'Destination', outFileWriter);
%    outFileWriter.close();
%
%

%   Copyright 2013 The MathWorks, Inc.

properties(GetAccess = public, SetAccess = private)
Filename;
TiffObject;
TileLength;
TileWidth;
end

methods

function obj = bigTiffWriter(fname, imageLength, imageWidth, tileLength, tileWidth)
% Constructor

validateattributes(fname,       {'char'},   {'row'});
validateattributes(imageLength, {'numeric'},{'scalar'});
validateattributes(imageWidth,  {'numeric'},{'scalar'});
validateattributes(tileLength,  {'numeric'},{'scalar'});
validateattributes(tileWidth,   {'numeric'},{'scalar'});

if(mod(tileLength,16)~=0 || mod(tileWidth,16)~=0)
error('bigTiffWriter:invalidTileSize',...
'Tile size must be a multiple of 16');
end

obj.Filename   = fname;
obj.ImageSize  = [imageLength, imageWidth, 1];
obj.TileLength = tileLength;
obj.TileWidth  = tileWidth;

% Create the Tiff object.
obj.TiffObject = Tiff(obj.Filename, 'w8');

% Setup the tiff file properties
% See "Exporting Image Data and Metadata to TIFF files
% http://www.mathworks.com/help/techdoc/import_export/f5-123068.html#br_c_iz-1
%
obj.TiffObject.setTag('ImageLength',   obj.ImageSize(1));
obj.TiffObject.setTag('ImageWidth',    obj.ImageSize(2));
obj.TiffObject.setTag('TileLength',    obj.TileLength);
obj.TiffObject.setTag('TileWidth',     obj.TileWidth);
obj.TiffObject.setTag('Photometric',   Tiff.Photometric.RGB);
obj.TiffObject.setTag('BitsPerSample', 8);
obj.TiffObject.setTag('SampleFormat',  Tiff.SampleFormat.UInt);
obj.TiffObject.setTag('SamplesPerPixel', 3);
obj.TiffObject.setTag('PlanarConfiguration', Tiff.PlanarConfiguration.Chunky);
end

function [] = writeRegion(obj, region_start, region_data)
% Write a block of data to the tiff file.

% Map region_start to a tile number.
tile_number = obj.TiffObject.computeTile(region_start);

% If region_data is greater than tile size, this function
% warns, else it will silently pad with 0s.
obj.TiffObject.writeEncodedTile(tile_number, region_data);

end

% Not implemented.
end

function close(obj)
% Close the tiff file
obj.TiffObject.close();
end

end
end


Let us run the included example:

Set file names and obtain size information from input file.

inFile        = 'example.tif';
inFileInfo    = imfinfo(inFile);
outFile       = 'out.tif';


Create an output TIFF file with tile size of 128x128

tileSize      = [128, 128]; % has to be a multiple of 16.
outFileWriter = bigTiffWriter(outFile, inFileInfo(1).Height, inFileInfo(2).Width, tileSize(1), tileSize(2));


Now call blockproc to rearrange the color channels.

blockproc(inFile, tileSize, @(b) flipdim(b.data,3), 'Destination', outFileWriter);
outFileWriter.close();


This adapter can serve as a starting point to write more complicated writer adapters for your particular tag combination of the Tiff format. Do you work with large image data? Do you, or would you use BigTIFF for your work? Do let us know!

Get the MATLAB code

Published with MATLAB® R2013a

## Homomorphic filtering – part 2

I'd like to welcome back guest blogger Spandan Tiwari for the second post in his two-part series on homomorphic filtering.

Last time we looked at how to apply a simple homomorphic filter. Today we continue our discussion on homomorphic filtering. First I'll load the variables I, H, and Ihmf that I computed last time.

load homomorphic_part1


In homomorphic filtering we apply a high-pass filter to the log-transformed image. The high-pass filtering step provides us with an opportunity to simultaneously apply other enhancements to the image. Consider a modified version of the high-pass filter $H(u,v)$ that we used last time.

$$H_{e}(u,v) = \alpha + \beta \ H(u,v)$$

We added an offset and a scaling factor for the Gaussian high-pass filter. If $\alpha < 1$ and $\beta > 1$, this filter will amplify the high-frequency components more than the low-frequency components. This filter is called high-frequency emphasis filter. The resulting image, typically, is sharper and also has better contrast. We choose $\alpha = 0.5$ and $\beta = 1.5$, and formulate the high-frequency emphasis filter.

alpha = 0.5;
beta = 1.5;
Hemphasis = alpha + beta*H;


Let's compare the original high-pass filter and the high-frequency emphasis filter by looking at their cross-sections.

plot(1:30,H(1,1:30),'r',1:30,Hemphasis(1,1:30),'b','LineWidth',2);
grid on;
legend('High-pass Filter','High-frequency Emphasis Filter','Location','best');


Now let's apply the filter and look at the result of homomorphic filtering. The image below shows the original (on the left) and the homomorphic filtered (on the right) images together. If you compare the two images you can see that the gradual change in illumination in the left image has been corrected to a large extent in the image on the right.

If = fft2(I, M, N);
Iout = real(ifft2(Hemphasis.*If));
Iout = Iout(1:size(I,1),1:size(I,2));

Ihmf_2 = exp(Iout) - 1;

imshowpair(I, Ihmf_2, 'montage')


The non-uniform illumination has largely been corrected. Now let's compare our earlier result of homomorphic filtering with regular high-pass filter (below, left) and the result with high-frequency emphasis filter (below, right). We see that the latter seems to have better non-uniform illumination compensation of the two.

imshowpair(Ihmf, Ihmf_2, 'montage')


Also, looking at these two images side-by-side highlights an interesting effect. In the image on the left there seems to be a bright halo type artifact on the borders. This can be seen more clearly if we increase the contrast by applying histogram equalization on the image on the left using the histeq function. Note that I am normalizing the image using mat2gray before passing it to histeq.

imshow(histeq(mat2gray(Ihmf)))


The halos on the borders can be seen clearly now. But if you look closely it is not just the image on the left that has the halo effect on the borders. The output from the high-frequency emphasis filter (image on the right) also has a similar, but less-pronounced, halo effect. Let's look at the histogram equalized version of both these images together to make this more apparent.

imshowpair(histeq(mat2gray(Ihmf)), histeq(mat2gray(Ihmf_2)), 'montage')


paddedI = padarray(I,ceil(size(I)/2)+1,'replicate');


Let's look at the result with the replicate style padding. I will show the results only for the high-frequency emphasis filter. We apply the high-frequency emphasis filter to the padded image.

If = fft2(paddedI);
Iout = real(ifft2(Hemphasis.*If));


Note that we padded the image on both sides in each dimension. This is unlike fft2 which appends trailing zeros along each dimension. Consequently, the output will also be padded on both sides.

imshow(Iout)


While cropping the image back to the original size we have to mindful of this and crop the image around the center of the larger output image.

Iout = Iout(ceil(M/2)-size(I,1)/2+1:ceil(M/2)+size(I,1)/2, ...
ceil(N/2)-size(I,2)/2+1:ceil(N/2)+size(I,2)/2);


Let's apply the exponential and look at the result. The image on the left is the output with zero padding that we computed earlier. The one of the right is the output with replicate style padding. We can see that the halo effect has been mitigated to a large extent.

Ihmf_3 = exp(Iout) - 1;
imshowpair(Ihmf_2, Ihmf_3, 'montage')


Let's use histogram equalization to get a better a better look.

imshowpair(histeq(mat2gray(Ihmf_2)), histeq(mat2gray(Ihmf_3)), 'montage')


The halo effect at the border has been mitigated significantly. But we had to do some code mechanics to get homomorphic filtering to play well in the FFT-domain. This is where the spatial domain filtering might offer a better alternative. It is relatively simpler to do the replicate style padding in the spatial domain. Let's create a simple spatial domain high-pass filter. Although we can create a spatial kernel exactly equivalent to the frequency-domain Gaussian high-pass filter we used earlier, we won't get into that to keep things simple. Here we make a simple ideal high-pass filter.

filterRadius = sigma;
hLowpass = fspecial('average', filterSize);
hImpulse = zeros(filterSize);
hHighpass = hImpulse - hLowpass;


Now we apply this high-pass filter in the spatial domain. To get the replicate style padding, we simply specify the 'replicate' option in imfilter. Then we apply the exponential to get the homomorphic filtered image.

Ihmf_spatial = imfilter(I, hHighpass, 'replicate');

Ihmf_spatial = exp(Ihmf_spatial) - 1;


Here's the filtered image (right) juxtaposed with the original input image (left).

imshowpair(I, Ihmf_spatial, 'montage')


Again, we see that the non-uniform illumination has been corrected to a large extent and there aren't any noticeable border artifacts. To show the absence of the border artifacts let's look at the histogram equalized version of the output. For comparison we will also enhance the output from frequency-domain high-emphasis filtered image (with zero-padding), which is shown on the left. Again, we can see that there aren't any noticeable border artifacts after filtering in the spatial domain (image on the right).

imshowpair(histeq(mat2gray(Ihmf_2)), histeq(mat2gray(Ihmf_spatial)), 'montage');


Both frequency domain and spatial domain filtering offer some practical advantages of their own. Spatial domain filtering offers a straightforward way to get the padding right. Also, it might be faster for small kernel sizes. But the frequency domain offers an easier and more intuitive way to get the filter with desired frequency characteristics.

In conclusion, homomorphic filtering is a useful tool to have in your quiver of enhancement techniques. It may not be applicable as a generic enhancement technique, but it works well on a certain set of problems.

Which applications have you used homomorphic filtering for? Were there any specific problems you faced while using it? Let us know.

Get the MATLAB code

Published with MATLAB® R2013a

## Homomorphic filtering – part 1

I'd like to welcome back guest blogger Spandan Tiwari for today's post. Spandan, a developer on the Image Processing Toolbox team, posted here previously about circle finding.

Recently I attended a friend's wedding at an interesting wedding venue, the Charles River Museum of Industry and Innovation in the town of Waltham, not very far from the Mathworks headquarters here in Natick. While waiting for the guests to arrive, we all got a chance to look around the museum. It's a nice little place with several quirky old machines on display, some dating back to 1812. Hanging on a wall in a dimly lit corner of the museum was a framed copy of an engineer's decadic logarithm (and anti-logarithm) tables. A friend who was looking around with me commented, "Ah, logs...haven't seen these since school days." That got me thinking. When was the last time I used logarithms? Well, just in the morning I had used the log function for scaling the Fourier magnitude spectrum of an image for better visualization. It took me a moment though to make the mental connection between the printed log tables and log's functional form, but as soon as I did, I realized that I used log quite often.

It also reminded me of a neat image enhancement technique from the olden days of image processing called homomorphic filtering, which makes clever use of the log function. There are several image denoising/enhancement techniques in literature that assume an additive noise model, i.e. $\tilde I(x,y) = I(x,y) + n(x,y)$, where $n$ is the noise signal. But there are relatively few techniques that work with other noise models, such as the multiplicative model $\tilde I(x,y) = I(x,y)\ n(x,y)$.

Homomorphic filtering is one such technique for removing multiplicative noise that has certain characteristics.

Homomorphic filtering is most commonly used for correcting non-uniform illumination in images. The illumination-reflectance model of image formation says that the intensity at any pixel, which is the amount of light reflected by a point on the object, is the product of the illumination of the scene and the reflectance of the object(s) in the scene, i.e.,

$$I(x,y) = L(x,y)\ R(x,y)$$

where $I$ is the image, $L$ is scene illumination, and $R$ is the scene reflectance. Reflectance $R$ arises from the properties of the scene objects themselves, but illumination $L$ results from the lighting conditions at the time of image capture. To compensate for the non-uniform illumination, the key is to remove the illumination component $L$ and keep only the reflectance component $R$. If we consider illumination as the noise signal (which we want to remove), this model is similar to the multiplicative noise model shown earlier.

Illumination typically varies slowly across the image as compared to reflectance which can change quite abruptly at object edges. This difference is the key to separating out the illumination component from the reflectance component. In homomorphic filtering we first transform the multiplicative components to additive components by moving to the log domain.

$$\ln(I(x,y)) = \ln(L(x,y)\ R(x,y))$$

$$\ln(I(x,y)) = \ln(L(x,y)) + \ln(R(x,y))$$

Then we use a high-pass filter in the log domain to remove the low-frequency illumination component while preserving the high-frequency reflectance component. The basic steps in homomorphic filtering are shown in the diagram below:

For a working example I will use an image from the Image Processing Toolbox.

I = imread('AT3_1m4_01.tif');
imshow(I)


In this image the background illumination changes gradually from the top-left corner to the bottom-right corner of the image. Let's use homomorphic filtering to correct this non-uniform illumination.

The first step is to convert the input image to the log domain. Before that we will also convert the image to floating-point type.

I = im2double(I);
I = log(1 + I);


The next step is to do high-pass filtering. We can do high-pass filtering in either the spatial or the spectral domain. Although they are both exactly equivalent, each domain offers some practical advantages of its own. We will do both it ways. Let's start with frequency-domain filtering. In frequency domain the homomorphic filtering process looks like:

First we will construct a frequency-domain high-pass filter. There are different types of high-pass filters you can construct, such as Gaussian, Butterworth, and Chebychev filters. We will construct a simple Gaussian high-pass filter directly in the frequency domain. In frequency domain filtering we have to careful about the wraparound error which comes from the fact that Discrete Fourier Transform treats a finite-length signal (such as the image) as an infinite-length periodic signal where the original finite-length signal represents one period of the signal. Therefore, there is interference from the non-zero parts of the adjacent copies of the signal. To avoid this, we will pad the image with zeros. Consequently, the size of the filter will also increase to match the size of the image.

M = 2*size(I,1) + 1;
N = 2*size(I,2) + 1;


Note that we can make the size of the filter (M,N) even numbers to speed-up the FFT computation, but we will skip that step for the sake of simplicity. Next, we choose a standard deviation for the Gaussian which determines the bandwidth of low-frequency band that will be filtered out.

sigma = 10;


And create the high-pass filter...

[X, Y] = meshgrid(1:N,1:M);
centerX = ceil(N/2);
centerY = ceil(M/2);
gaussianNumerator = (X - centerX).^2 + (Y - centerY).^2;
H = exp(-gaussianNumerator./(2*sigma.^2));
H = 1 - H;


Couple of things to note here. First, we formulate a low-pass filter and then subtracted it from 1 to get the high-pass filter. Second, this is a centered filter in that the zero-frequency is at the center.

imshow(H,'InitialMagnification',25)


We can rearrange the filter in the uncentered format using fftshift.

H = fftshift(H);


Next, we high-pass filter the log-transformed image in the frequency domain. First we compute the FFT of the log-transformed image with zero-padding using the fft2 syntax that allows us to simply pass in size of the padded image. Then we apply the high-pass filter and compute the inverse-FFT. Finally, we crop the image back to the original unpadded size.

If = fft2(I, M, N);
Iout = real(ifft2(H.*If));
Iout = Iout(1:size(I,1),1:size(I,2));


The last step is to apply the exponential function to invert the log-transform and get the homomorphic filtered image.

Ihmf = exp(Iout) - 1;


Let's look at the original and the homomorphic-filtered images together. The original image is on the left and the filtered image is on the right. If you compare the two images you can see that the gradual change in illumination in the left image has been corrected to a large extent in the image on the right.

imshowpair(I, Ihmf, 'montage')


Next time we will explore homomorphic filtering some more. We will see a slightly modified form of homomorphic filtering and also discuss why we might want to do filtering in the spatial domain.

Get the MATLAB code

Published with MATLAB® R2013a

## R2013a – image processing and computer vision

I've mentioned the R2013a release in two previous posts (15-May-2013 and 12-Mar-2013). Today I want to point out that R2013a is a pretty significant release in terms of new features related to image processing and computer vision.

Here's the quick summary of what's new in Image Processing Toolbox, Computer Vision System Toolbox, and Image Acquisition Toolbox.

Image Processing Toolbox version 8.2

• Image segmentation using active contours
• Classes and functions for representing and applying 2-D and 3-D geometric transformations
• Classes for defining the world coordinate system of an image
• Code generation for conndef, imcomplement, imfill, imhmax, imhmin, imreconstruct, imregionalmax, imregionalmin, iptcheckconn, and padarray functions (using MATLAB Coder)
• GPU acceleration for imrotate, imfilter, imdilate, imerode, imopen, imclose, imtophat, imbothat, imshow, padarray, and bwlookup functions (using Parallel Computing Toolbox)

Computer Vision System Toolbox version 5.2

• Cascade object detector training using Haar, Histogram of Oriented Gradients (HOG), and Local Binary Patterns (LBP) features
• Fast Retina Keypoint (FREAK) algorithm for feature extraction
• Hamming distance method for matching features
• Multicore support in matchFeatures function and Foreground Detector System object
• Functions for corner detection, geometric transformation estimation, and text and graphics overlay, augmenting similar System objects

Image Acquisition Toolbox version 4.5

• Kinect for Windows sensor support for acquiring images, depth maps, skeleton data, and related metadata

I've written extensively here in the past about geometric transformations (the spatial transformations category includes 21 different posts over a three-year period). Now there is much more to say. The new geometric transformation functionality represents a ground-up redesign, based partially on customer response and feedback to the functionality (such as imtransform) added way back in 2001. I hope to recruit a developer on the Image Processing Toolbox team to post here soon about the changes.

In the meantime, here's an example from the documentation for the new function imsharpen.

a = imread('hestain.png');
imshow(a)
title('Original Image');

b = imsharpen(a);
imshow(b)
title('Sharpened Image');


Get the MATLAB code

Published with MATLAB® R2013a

Steve Eddins is a software development manager in the MATLAB and image processing areas at MathWorks. Steve coauthored Digital Image Processing Using MATLAB. He writes here about image processing concepts, algorithm implementations, and MATLAB.

These postings are the author's and don't necessarily represent the opinions of MathWorks.