Lab 4.3 Image analysis: Regions of interest

Introduction

Now that we have seen how image data is represented digitally, it's time to do some image analysis.

One common task in image analysis is to identify interesting objects in a field and to calculate and perform statistics of some parameter, such as the object's brightness or size.

Manual selection of regions of interest

The simplest way to identify a region of interest (abbreviated ROI) is to use your own visual system and the mouse to select a region. In this example, we'll keep the raw image separate from the image that we will use for display so that we do not alter the underlying data values when we perform analysis. We will the function roipoly to obtain a series of mouse clicks that define the object (see help roipoly). When you are finished selecting the object, either double-click or click on an existing point to close the polygon.

As you are working through this example, pick out a cell in the cellimage.tiff to draw.

im_raw = double(imread('cell_image.tiff'));

im_disp = dasw.math.rescale(im_raw, [min(im_raw(:)) max(im_raw(:))],[0 255]);

figure;

colormap(gray(256));

image(im_disp);

% at this point, you may want to zoom in so you can see a cell well

zoom on

% be sure you are done zooming before you call roipoly below

[BW,xi,yi] = roipoly;

hold on

plot(xi,yi,'b-');

You may notice that the points xi and yi define a contour around the region of interest. Now let's examine the mask image that is returned in variable BW. BW is an image of 0's and 1's that is the same size as im_disp; it is 1 for points that are inside the region of interest, and 0 for points outside.

size(BW)

max(BW(:))

figure;

colormap(gray(256))

image(BW*255); % let's make the brightest pixel white (255)

In Matlab, one generally specifies a region of interest as a collection of image matrix index values

We have seen 2 ways of representing a region of interest: a contour and a mask image. However, to actually examine the brightness values in pixels inside the region of interest, we need to identify the matrix index values that correspond to these pixels. (Note that here I am talking about the matrix index values, not the index values of a color look-up table.) We can use the find command for this purpose:

[myroi_rowinds, myroi_colinds] = find(BW); % find all non-zero points

% let's examine the contents of this variable

myroi_rowinds

myroi_colinds

We can use these index values to reference the brightness values of pixels within the region of interest:

im_raw(myroi_rowinds, myroi_colinds)

Matrix index values in 2-dimensional matrices (or n-dimensional matrices)

In Matlab, we have several options for referring to the elements of a 2-dimensional matrix. Let's examine this in more detail.

Let's create a matrix A:

A = [ 13 14 15 ; 16 17 18 ; 19 20 21 ; 22 23 24]

We can refer to a portion of A by providing the 2-dimensional matrix index values that we want to access A(rows, columns):

A(1:2,2:3)

However, we also have the ability to reference the ith element of A, where i runs from 1 to the number of elements of A, which is equal to the number of rows times the number of columns (4 * 3 = 12 in the case of our matrix here). Let's see how this works:

A(1)

A(2)

A(3)

A(4)

A(5)

A(6)

A(7)

You can see that we are working our way through the columns of the matrix.

Suppose we want to identify the pixel that has the value 20. We can do it in 2 ways:

[rowinds,colinds] = find(A==20)

[inds] = find(A==20)

Both will find the correct pixel:

A(rowinds,colinds)

A(inds)

In imaging processing, it is sometimes more convenient to use the 1-dimensional version, and sometimes it is more convenient to use the 2-dimensional version. We will use both.

Examining simple properties of a region of interest

We are now in a position to examine some simple features of our ROI. For example, we can obtain the average brightness. For this purpose, it is easier to use the 1-dimensional index values:

myroi_inds = find(BW); % find non-zero elements

avg_brightness = mean(im_raw(myroi_inds))

We can examine the size of the ROI in terms of number of pixels:

disp(['Number of pixels is ' int2str(length(myroi_inds)) '.']);

We can find the maximum and minimum brightness values:

mx = max(im_raw(myroi_inds));

mn = min(im_raw(myroi_inds));

disp(['The maximum value is ' num2str(mx) '.']);

disp(['The minimum value is ' num2str(mn) '.']);

Automatic selection of regions of interest

While manual drawing has the advantage that the user can carefully inspect the ROI as it is being drawn, it is often desirable to have some means of automatically identifying interesting objects within an image.

Thresholding

Let's consider the following images from Sue Paradis' lab. Suppose we want to automatically detect the locations of synaptic proteins PSD-95 and synapsin, which have been labeled with fluorescent antibodies in these images. The cell of interest is labeled with green fluorescent protein (GFP).

First, let's look at the images. They have a little twist to them that we haven't seen yet; they contain multiple frames that we can read out as follows:

img_psd95 = imread('CELL10_DSP.tif',1);

img_cell = imread('CELL10_DSP.tif',2);

img_synapsin = imread('CELL10_DSP.tif',3);

figure(100);

colormap(gray(256));

subplot(2,2,1);

image(img_psd95);

subplot(2,2,2);

image(img_cell);

subplot(2,2,3);

image(img_synapsin);

% now let's get fancy and make a full color combined image

fullcolor = cat(3,img_psd95, img_cell, img_synapsin);

subplot(2,2,4);

image(fullcolor);

(When you're done, leave this figure open!)

Zoom around the images, particularly the combined image. Note that there are some spots that only exhibit PSD-95 staining (red), others that only exhibit synapsin staining (blue), and other areas that exhibit both staining (purple, or, when this location overlaps the green cell, white).

Suppose we want to detect all of the PSD-95 spots. Let's look carefully at the PSD-95 image above. To my eye, it has some very bright spots, and then some dimmer areas that could be very weak staining, or an imaging artifact. We'd really like to detect PSD-95 spots, but ignore the dimmer areas in the image. To do this, we can apply a threshold to the image and only examine pixels that exceed a certain brightness.

To determine the value of the threshold we should pick, we should examine the intensity histogram:

figure;

hist(img_psd95(:), 1:255);

Zoom around this figure and look carefully. Because so much of the image is completely dark, the largest bin by far is the 0 bin. If you zoom in and look carefully, there are several thousand pixels that have the value 255. However, there is a huge gap between intensity values 255 and about value 51, where we start to see some pixels again. Zoom in on the region between 1 and 51. There are a few pixels with these values that we'd like to consider as if they were 0. To do this, we apply a threshold:

img_psd95_thresh = img_psd95;

myinds = find(img_psd95 < 60);

img_psd95_thresh(myinds) = 0;

figure;

colormap(gray);

image(img_psd95_thresh);

Detecting spots in "binary" images

The thresholding process has created what is called a binary image where pixels that are determined to contain "signal" are non-zero and those that are 0 are said to contain non-signal. Now we can move on to detect the spots.

Because we might want to both plot and examine the brightness of our ROIs, we want procedures that will return both contours and the matrix index values for our spots. Fortunately, Matlab has these routines to do both built into the images toolbox (help images):

In order to understand how these routines work, we'll run them on a toy example. (You can cut and paste this:)

THRESH = ([1 1 1 0 0 0 0 0

1 1 1 0 1 1 0 0

1 1 1 0 1 1 0 0

1 1 1 0 0 0 1 0

1 1 1 0 0 0 1 0

1 1 1 0 0 0 1 0

1 1 1 0 0 1 1 0

1 1 1 0 0 0 0 0]);

The first thing we will do is to detect the objects and trace their boundaries:

[BW, L] = bwboundaries(THRESH,4);

Let's examine the outputs. L is a "labeled" matrix where each pixel is labeled with its corresponding object number (there are 3 objects in this matrix):

L,

BW is a cell list of X/Y pairs that define the contours around these objects. Look at

BW{1}

BW{2}

Now let's draw it:

figure;

colormap(gray(256));

subplot(2,2,1);

image(255*THRESH);

subplot(2,2,2);

image(128*L);

subplot(2,2,1);

hold on

plot(BW{1}(:,2),BW{1}(:,1),'-r');

plot(BW{2}(:,2),BW{2}(:,1),'-b');

plot(BW{3}(:,2),BW{3}(:,1),'-g');

We can use another routine called regionprops to identify the pixel index values that correspond to the labeled regions in L:

stats = regionprops(L,'PixelIdxList')

Let's examine this:

stats(1)

stats(1).PixelIdxList

So it does contain the pixel index values for the first object. Now let's create a new image where we give these pixel values specific gray values to show it really corresponds:

newimg = zeros(size(L));

subplot(2,2,3);

image(newimg);

subplot(2,2,4);

newimg(stats(1).PixelIdxList) = 255;

newimg(stats(2).PixelIdxList) = 128;

newimg(stats(3).PixelIdxList) = 64;

image(newimg);

Now let's apply this to the real image data. Create a new m-file called spotdetector.m. This function has some general purpose bells and whistles to make data analysis easier. Rather than returning just the contours and index values, it returns a structure that you can add on to (put this .m file in a new package called +dasw/+roi).

Now let's run this function on our image:

[rois_psd95,L] = dasw.roi.spotdetector(img_psd95_thresh,4,'CELL10',1,{'psd95'});

Next let's write a simple .m file to plot the ROIs (place in +dasw/+roi package):

Let's plot it on the raw data in figure 100 from above:

figure(100);

subplot(2,2,1);

hold on;

[h_lines,h_text]=dasw.roi.plot_rois(rois_psd95,9,[1 0 0]);

Zoom around and check this out. It will probably make your computer's brain hurt a little bit because there are a lot of ROIs. After checking it out we can delete these plots to save your computer some trouble:

delete(h_lines)

delete(h_text)

Examining properties of regions of interest using the Matlab image toolbox

We snuck in a call to Matlab's regionprops command in our function spotdetector. It returned a bunch of statistics on these spots. Let's look at just a couple of these.

To see the list, check out help regionprops. You can also look at the first ROI's data structure:

rois_psd95(1)

rois_psd95(1).stats

Now just for fun let's plot the distribution of areas:

stats=[rois_psd95.stats]; % this pulls stats fields into separate variable...

area = [stats.Area]; % so we can do this in one more line to get areas

Now let's plot a histogram of the areas:

figure;

hist(area,100);

xlabel('Area (pixels^2)');

ylabel('Number of spots');

title('Spots of PSD-95');

Simple colocalization analysis

Suppose we wanted to identify the subset of PSD-95 that overlapped the cell that is expressing GFP (to make an attempt to distinguish synapses that might be on the cell, or at least very close to it). There are many ways to do this, but the simplest might be to check, for each ROI, if there is any overlap in a thresholded green image.

First let's write a little function roi_overlap.m to check for overlap (put it in +dasw/+roi)

Now we can get the overlaps and plot overlapping ROIs on the green image (which we will threshold above 10). In order to do this, you will need to put the function label_rois.m from the zip file below into +dasw/+roi to create the package function dasw.roi.label_rois.m.

overlaps = dasw.roi.roi_overlap(rois_psd95, (img_cell>10));

inds = find(overlaps>0);

% let's use the provided function roi_label to label these as GFP+

rois_psd95(inds) = dasw.roi.label_rois(rois_psd95(inds), 'GFP');

% now let's plot this subset on the GFP image

figure(100);

subplot(2,2,2);

[h_lines,h_text]=dasw.roi.plot_rois(rois_psd95(inds),9,[1 0 0]);

disp(['Number of overlapping ROIs is ' int2str(length(inds)) '.']);

Now just for exercise, let's plot the (cartesian, not dendritic) distances of all of these ROIs from the soma. Get ready to click on the soma in the white image (upper right):

figure(100);

subplot(2,2,2);

delete(h_lines);

delete(h_text);

somapt = ginput(1) % you'll click here

stats = [rois_psd95(inds).stats];

centroids = [stats.Centroid]; % this is now 1-d, got to convert to 2-d

centroids = reshape(centroids,2,length(centroids)/2)';

diffs = centroids-repmat(somapt,length(centroids),1);

dists = sqrt(sum(diffs.*diffs,2));

figure;

hist(dists,100);

xlabel('Distances (pixels)');

ylabel('Number of examples');

References:

The function spotdetector is loosely based on a function from Sri Raghavachari.

Function reference:

Matlab functions introduced

(type help functionname for help):

• bwboundaries

• regionprops

User functions introduced:

• dasw.roi.spotdetector

• dasw.roi.plot_rois

• dasw.roi.roi_overlap

• dasw.roi.label_roi

User functions used:

• dasw.math.rescale - from Lab4.2

Labs are copyright 2011-2012, Stephen D. Van Hooser, all rights reserved.