BME355 Lab Listing: Image Processing
Top navigation banner.


Image
Processing



Preparation

Lab Outline

Procedure

Learning more


Sample Images

GOALS

To understand the basics of image processing, including:
  • sampling, quantization and the effect of noise on two-dimensional signals (images).
  • image enhancement and restoration techniques
  • image reconstruction from projections

    Introduction

    Biomedical images arise from a variety of modalities. The radiologic modalities for imaging the body include: X-rays (planar projections and computed tomography), ultrasound, magnetic resonance, SPECT (single photon emission computed tomography) and PET (positron emission tomography). Microscopes, including confocal and electron, are a source of imagery at the cellular level and smaller. Conventional light cameras are also used for imaging the surface of the body and the interior of the body using endoscopes (with fiber optics). Infrared cameras can be used for thermography (measuring heat distribution on the body).

    Image processing is the natural extension of signal processing (1D signals) to images (2D, 3D, 4D signals).

    Quality

    Image quality can be characterized in a number of ways: brightness, contrast, sharpness, noisiness. Contrast is the relative change in intensity between two regions, typically neighboring. The sharper an image, the steeper the transition at a boundary. Noise can be characterized by the variance of the gray levels within a uniform area.

    Quantization

    Images are typically quantized to 256 levels (8 bits) although biomedical images may have many more levels. While the human eye can distinguish fewer than 100 gray levels in an image, with image processing we can take advantage of additional quantization. Some medical images may have up to 65536 gray levels (16 bits). When there are too few quantization levels in a smooth, gradually shaded region, you will see "false contours" at the arbitrary curve where there is a transition from one gray level to the next.

    Fourier Transform

    The Fourier transform of a two-dimensional signal
    F(u,v) = F { f(x,y) } is the weight or magnitude of each component frequency and is computed by this double integral:

    F(u,v) = F { f(x,y)) } =




    - 



    - 
    f(x,y) e-2pi (x u + v y) dx dy

    The zero frequency value is just the average value of the signal.

    Inverse Fourier Transform

    We can recover the original signal from its frequency domain representation using:

    f(x,y) = F-1 { F(u,v) } =




    - 



    - 
    F(u,v) e2pi (x u + y v) du dv

    f is a linear combination of elementary periodic patterns
    e2pi (xu + yv) (sinusoids). F is a measure of the relative contributions of each sinusoid.




    Lines corresponding to maxima of cos2 p(u x + v y), with the angle and the spacing determined by u and v. This corresponds to spatial frequencies (u,v). The period is
    1 / [(u2 + v2)].



    The image is decomposed into frequencies in both the x and the y direction. High and low frequencies correspond to sinusoidal variation in a particular direction through the image. Sharp transitions of gray level in an image require high frequency variation.

    When you have a stack of two dimensional images (such as slices through a body), you have a three-dimensional image, which would have a three dimensional Fourier transform. You can even have a four dimensional image when you have a time sequence of three dimensional images, such as of the beating heart.

    Sampling

    Sampling issues are the same as for one dimensional signals except now you must consider two sampling intervals (usually these are equal) and the frequency content in both x and y. Each sample in a image is called a pixel (short for picture element).

    Enhancement

    One way to enhance images is to use linear filtering, as with one-dimensional signals. Lowpass filtering can be used to remove high frequency noise. The resulting image will look smoother or blurrier because of the attenuation of the high frequency content of the image. Highpass filtering can be used to sharpen an image. If you take an image and add a highpass filtered image to it, the result will be sharper because the transitions (edges) in the image have been strengthened. These filters can often be implemented simply in the discrete space domain by convolving with small kernels. This is equivalent to multiplying in the Fourier domain.

    Discrete convolution is computed by centering the flipped filter at each pixel in the image and replacing the pixel value with the sum of the products of the corresponding elements.

    The simplest low pass filter is just the average of each pixel in a 3x3 neighborhood:

    flow_pass (i,j) = f(i,j) ** 1
    9




    1
    1
    1
    1
    1
    1
    1
    1
    1




    The simplest high pass filter is the negative of the Laplacian (the sum of the second derivatives in x and y):

    fhigh_pass (i,j) = f(i,j) **



    0
    - 1
    0
    -1
    4
    -1
    0
    -1
    0




    This is like taking the image and subtracting a blurred or low pass version of it (averaging the four neighbors). Since we are subtracting the low frequency content, the high frequency content remains.

    Enhancing the edges in images sometimes improves their quality. Edges are sharp transitions in gray levels and are measured by differentiation. Vertical edges correspond to derivatives in the x direction; horizontal edges are derivatives in the y direction. These can be computed by discrete approximation with simple filters:

    fvert_edges(i,j) = f(i,j) ** 1
    6




    -1
    0
    1
    -1
    0
    1
    -1
    0
    1




    fhorz_edges(i,j) = f(i,j) ** 1
    6




    -1
    -1
    -1
    0
    0
    0
    1
    1
    1




    The gradient magnitude (edge strength) is then just:
    grad_mag =

     

    fhorz_edges2 + fvert_edges2
     

    Median Filter

    Instead of linear filtering, there are many ways to do nonlinear filtering for enhancement. One example is the median filter. At each point in the image, replace the value with the median computed from a neighborhood around the point (e.g. 3x3 or 5x5). This will reduce or remove "salt and pepper" noise, that is, isolated pixels set to black or white, due to transmission or other errors, without excessive blurring.

    Histogram Equalization

    The image histogram is a histogram or count of the number of pixels at each gray level. If the histogram is concentrated in one portion of the gray level range, the image contrast will not be good. If concentrated at the high end, the image will appear overexposed; at the low end, underexposed. In terms of histograms, the best overall contrast is when the histogram is flat or equalized. This means that the full contrast range is best utilized. The gray levels can be transformed to result in a flat histogram for image enhancement. Note that unless the pixel values are divided up (i.e. assigned to different output levels) equalization cannot result in a perfectly flat histogram.

    Noise

    Gaussian noise is typically modeled as additive (although sometimes multiplicative) and independent of the signal.

    For Poisson noise, such as results from nuclear events, the signal is governed by p(k,m) = [(e-mmk)/k!]. The average number of events observed equals m = pq = prob of event ×time. The average equals m, i.e. k k p(k,m) = m. The variance equals m, i.e. s2 = k (k-m)2 p(k,m) = m. Thus, the noise is signal dependent.

    Salt and pepper noise occurs due to transmission or other errors that result in either a white or black pixel.

    Restoration

    The goal of image restoration is to remove degradations. Unlike enhancement, we explicitly model the degradation. A simple model of degradation is a convolution with a blurring function plus additive noise:

    degraded_image = blur **   true_image + noise

    The blurring function must be known either from the physics of the camera or imager or could be determined empirically. For some techniques, statistical properties of the noise must be known.

    Inverse Filtering

    The simplest restoration technique is inverse filtering. Here, we ignore noise and simply try to undo the blurring effect. Consider,

    g(x,y) = f(x,y) **  h(x,y)

    where g is the degraded image, f is the true image and h is the blurring or point spread function. In the Fourier domain, this would be:

    G(u,v) = F(u,v) H(u,v)

    where G, F and H are the Fourier transforms of g, f and h, respectively. It is natural to determine an estimate [F'] of F from G using:

    [F'](u,v) = G(u,v) [1/H(u,v)]

    Thus, the image restored by this inverse filter would be:

    [f'](x,y) = F { G(u,v) [1/H(u,v)] }

    Note that where H is zero, this is not defined. Where H is zero, the corresponding frequencies have been completely lost and the inverse filter cannot restore these frequencies. Noise can be a real problem, especially where H is small, because the inverse filter will make it even worse. Thus, in practice, the inverse filter is often set to zero when H is less than some threshold. More sophisticated restorations account for noise explicitly, e.g. Wiener Filtering.

    Reconstruction from Projections

    Images of slices through the body can be produced using tomographic imaging by reconstructing the image from projections. The primary examples in medical imaging are CT (computed tomography) which uses X-rays, and nuclear medicine, including SPECT (single photon emission computed tomography) and PET (positron emission tomography), which use injected radioactive compounds (like the gamma camera).

    As X-rays pass through the body, they are absorbed according to the local density and composition. The resulting intensity is a function of the integral of this absorption (m(x,y)) along the straight line path through the body:


    m(x,y) dS = log I0
    Id
    (1)
    where I0 is the incident X-ray intensity and Id is the detected intensity.

    If a set of projections of an image f(x,y) is taken at the same angle f with varying position t, the Fourier transform of this function, Pf(t) would be:

    Sf(w) =
    Pf(t) e-2pi wt dt
    (2)
    Radon showed that this was equal to:
    Sf(w) =

    f(x,y) e-2pi w(x cosf+ y sinf) dx dy
    (3)
    which is just the Fourier transform of the image along the line of the projection, F(wcosf, wsinf).


    In principle, with projections along many angles, we can reconstruct the Fourier transform (and therefore the image itself) of the slice. In practice, other algorithms, using this relationship, can be developed.

    The most popular reconstruction algorithm is filtered back-projection. By manipulating the above, we can write the image as:

    f(x,y) =
    p

    0 




     
    Sf(w) |w| e2pi wt dw
    d f
    (4)
    where t = x cosf+ y sinf. The term in brackets, called Qf(t), is the filtered version of the projection data (Fourier transformed, multiplied by filter, inverse Fourier transformed). The remaining operation:
    f(x,y) =
    p

    0 
    [ Qf (t) ] d f
    (5)
    is called backprojection because it is equivalent to smearing each Qf(t) at angle f across the image and adding them up.

    File translated from TEX by TTH, version 2.58. On 8 Dec 1999, 11:05.
    Bottom navigation banner.