Package gov.nih.mipav.model.algorithms.filters

Class AlgorithmNLNoiseReduction

java.lang.Object
java.lang.Thread
gov.nih.mipav.model.algorithms.AlgorithmBase
gov.nih.mipav.model.algorithms.filters.AlgorithmNLNoiseReduction
All Implemented Interfaces:
ActionListener, WindowListener, Runnable, EventListener
public class AlgorithmNLNoiseReduction extends AlgorithmBase

This is a port of the SUSAN Nonlinear Noise reduction program. SUSAN noise reduction uses nonlinear filtering to reduce noise in an image (2D or 3D) while preserving the underlying structure. It does this by only averaging a voxel with local voxels which have similar intensity. The original code could only handle integer data going from 0 to 255. 2 main changes were made to allow any type of black and white data. In setupBrightness temp = ((float)k)/(float)brightThres); was changed to: temp=(float)((k*(srcMax-srcMin))/(bpSize * brightThres)); In smoothing and smoothing3d bp[Math.abs(center-brightness)] was changed to: bp[(int)Math.abs(256*(center-brightness)/(srcMax - srcMin))]; Also note that using the Math.abs allowed changing the bp allocation from bpSize*2 + 1 to bpSize + 1. In the dialog the default value of brightThres was set equal to 0.1 * (srcMax - srcMin).

The area of a mask excluding the center which has a same or similar intensity value as the center of the mask is known as the USAN. The SUSAN filter works by taking an average of the pixels in the locality which lie in the USAN. The SUSAN filter uses a Gaussian in both the brightness and spatial domains. If useMedian is true, the sums taken over the local neighborhod do not include the center pixel itself. This allows good reduction of impulse noise. If the USAN area is zero and useMedian is true, then the median of the pixel's 8 nearest neighbors in 2D or 26 nearest neighbors in 3D is used to estimate the pixel's correct value. If useMedian is false, then the sums include the center pixel. If useMedian is false and the USAN is zero, then the pixel remains unchanged.

For useMedian == true the complete equation for the SUSAN filter is: J(x,y) = num/denom, where num = sum(i,j)!=(0,0)I(x+i,y+j)*exp[-(r*r/2*sigma*sigma) - (I(x+i,y+j) - I(x,y))**2/thresh*thresh] denom = sum(i,j)!=(0,0) exp[-(r*r/2*sigma*sigma) - (I(x+i,y+j) - I(x,y))**2/thresh*thresh] with r = sqrt(i*i + j*j) and the above replacement rule is used when the denominator or USAN is zero. If useMedian == false, the restriction (i,j) != (0,0) in the above is removed.

The SUSAN filter tends to preserve edges. Pixels are pulled toward the neighboring region to which they are closest in value.

SUSAN Version 2i_medx - nonlinear 2D/3D smoothing

Oxford center for Functional Magnetic Resonance Imaging of the Brain, Department of Clinical Neurology, Oxford University, Oxford, UK (Previously in Computer Vision and Image Processing Group - now Computer Vision and Electro Optics Group - DERA Chertsey, UK) Email: steve@fmrib.ox.ac.uk WWW: http://www.fmrib.ox.ac.uk/~steve

(C) Crown Copyright (1995-1999), Defence Evaluation and Research Agency, Farnborough, Hampshire, GU14 6TD, UK DERA WWW site: http://www.dera.gov.uk/ DERA Computer Vision and Electro Optics Group WWW site: http://www.dera.gov.uk/imageprocessing/dera/group_home.html DERA Computer Vision and Electro Optics Group point of contact: Dr. John Savage, jtsavage@dera.gov.uk, +44 1344 633203

A UK patent has been granted: "Method for digitally processing images to determine the position of edges and/or corners therein for guidance of unmanned vehicle", UK Patent 2272285. Proprietor: Secretary of State for Defence, UK. 15 January 1997

This code is issued for research purposes only and remains the property of the UK Secretary of State for Defence. This code must not be passed on without this header information being kept intact. This code must not be sold.

This 3D version derived from 2d version SUSAN Version 2i SUSAN = Smallest Univalue Segment Assimilating Nucleus Email: steve@fmrib.ox.ac.uk WWW: http://www.fmrib.ox.ac.uk/~steve Related paper: article{Smith97, author = "Smith, S.M. and Brady, J.M.", title = "{SUSAN} - A New Approach to Low Level Image Processing", journal = "Int. Journal of Computer Vision", pages = "45--78", volume = "23", number = "1", month = "May", year = 1997}

To be registered for automatic (bug) updates of SUSAN, send an email.

Note: FDT is the image data type

Note: edge and corner finding have been taken out of the 3D version of SUSAN - it is not thought that they would be of great value.

See following section for different machine information. Please report any bugs (and fixes). There are a few optional changes that can be made in the "defines" section which follows shortly.

This code is written using an emacs folding mode, making moving around the different sections very easy. This is why there are various marks within comments and why comments are indented.

SPATIAL CONTROL: d

In SUSAN smoothing d controls the size of the Gaussian mask; its default is 4.0. Increasing d gives more smoothing. In edge finding, a fixed flat mask is used, either 37 pixels arranged in a "circle" (default), or a 3 by 3 mask which gives finer detail. In corner finding, only the larger 37 pixel mask is used; d is not variable. In smoothing, the flat 3 by 3 mask can be used instead of a larger Gaussian mask; this gives low smoothing and fast operation.

BRIGHTNESS CONTROL: t

In all three algorithms, t can be varied (default=20); this is the main threshold to be varied. It determines the maximum difference in greylevels between two pixels which allows them to be considered part of the same "region" in the image. Thus it can be reduced to give more edges or corners, i.e. to be more sensitive, and vice versa. In smoothing, reducing t gives less smoothing, and vice versa. Set t=10 for the test image available from the SUSAN web page.

ITERATIONS:

With SUSAN smoothing, more smoothing can also be obtained by iterating the algorithm several times. This has a different effect from varying d or t.

BRIGHTNESS FUNCTION LUT IMPLEMENTATION: (Only read this if you are interested in the C implementation)

The SUSAN brightness function is implemented as a LUT (Look-Up-Table) for speed. The resulting pointer-based code is a little hard to follow, so here is a brief explanation. In setupBrightness() the LUT is setup. This mallocs enough space for *bp and then repositions the pointer to the center of the malloced space. The SUSAN function e^-(x^6) or e^-(x^2) is calculated and converted to a unsigned char in the range 0-100, for all possible image brightness differences (including negative ones). Thus bp[23] is the output for a brightness difference of 23 greylevels. In the SUSAN algorithms this LUT is used as follows:

p=in + (i-3)*x_size + j - 1; p points to the first image pixel in the circular mask surrounding point (x,y).

cp=bp + in[i*x_size+j]; cp points to a position in the LUT corresponding to the brightness of the center pixel (x,y).

now for every pixel within the mask surrounding (x,y), n+=*(cp-*p++); the brightness difference function is found by moving the cp pointer down by an amount equal to the value of the pixel pointed to by p, thus subtracting the two brightness values and performing the exponential function. This value is added to n, the running USAN area.

in SUSAN smoothing, the variable height mask is implemented by multiplying the above by the moving mask pointer, reset for each new center pixel. tmp = *dpt++ * *(cp-brightness);

  • Field Details

    • brightnessTable

      private byte[] brightnessTable
      DOCUMENT ME!

    • brightThres

      private final double brightThres
      DOCUMENT ME!

    • bTableSize

      private int bTableSize
      DOCUMENT ME!

    • dimension

      private int dimension
      DOCUMENT ME!

    • gaussMask

      private byte[] gaussMask
      DOCUMENT ME!

    • inVolume

      private float[][] inVolume
      DOCUMENT ME!

    • inVolume3d

      private float[] inVolume3d
      DOCUMENT ME!

    • maskSize

      private int maskSize
      DOCUMENT ME!

    • xMaskSize

      private int xMaskSize
      DOCUMENT ME!

    • yMaskSize

      private int yMaskSize
      DOCUMENT ME!

    • zMaskSize

      private int zMaskSize
      DOCUMENT ME!

    • maskStdDev

      private final float maskStdDev
      DOCUMENT ME!

    • medianVal

      private float[] medianVal
      DOCUMENT ME!

    • newProgressValue

      private int newProgressValue
      DOCUMENT ME!

    • progressValue

      private int progressValue
      DOCUMENT ME!

    • sliceSize

      private int sliceSize
      DOCUMENT ME!

    • smallSliceSize

      private int smallSliceSize
      DOCUMENT ME!

    • srcMin

      private double srcMin
      DOCUMENT ME!

    • srcMax

      private double srcMax
      DOCUMENT ME!

    • threeByThree

      private boolean threeByThree
      DOCUMENT ME!

    • tmpImage

      private float[] tmpImage
      DOCUMENT ME!

    • useMedian

      private boolean useMedian
      DOCUMENT ME!

    • volSize

      private int volSize
      DOCUMENT ME!

    • xSize

      private int xSize
      DOCUMENT ME!

    • ySize

      private int ySize
      DOCUMENT ME!

    • zSize

      private int zSize
      DOCUMENT ME!

    • xStdDev

      private float xStdDev
      DOCUMENT ME!

    • yStdDev

      private float yStdDev
      DOCUMENT ME!

    • zStdDev

      private float zStdDev
      DOCUMENT ME!

  • Constructor Details

    • AlgorithmNLNoiseReduction

      public AlgorithmNLNoiseReduction(ModelImage srcImg, double bt, float maskSD, boolean useMedian, boolean img25D)
      AlgorithmNLNoiseReduction - Constructor.
      Parameters:
      srcImg - source image model
      bt - Brightness threshold: this allows discrimination between noise and the underlying image. Ideally, the value should be set greater than the noise level and less than the contrast of the underlying image. Edges of contrast smaller than this threshold will be blurred whereas those of greater contrast will not be. Reducing bt gives less smoothing. bt has a default value = 0.1 * (srcMax - srcMin).
      maskSD - This determines the spatial extent of the smoothing. The mask is basically Gaussian with standard deviation (in image units - e.g. mm.) set by the user. However, for a small, fast, flat response with a 3x3 or 3x3x3 voxel mask, set maskSD to zero. maskSD has a default value equal to the x resolution.
      useMedian - If true, the center pixel is not included in the sums. When the local neighborhood of similar brightness voxels is empty, a local median filter is used. This allows the correction of impulse("salt-and-pepper")noise. If false, this feature is turned off. In this case, the center pixel is included in the sums. When no neighborhood is found, the original intensity of the voxel of interest remains unchanged.
      img25D - Flag, if true, indicates that each slice of the 3D volume should be processed independently. 2D images disregard this flag.

    • AlgorithmNLNoiseReduction

      public AlgorithmNLNoiseReduction(ModelImage destImg, ModelImage srcImg, double bt, float maskSD, boolean useMedian, boolean img25D)
      AlgorithmNLNoiseReduction - Constructor.
      Parameters:
      destImg - image model where result image is to stored
      srcImg - source image model
      bt - Brightness threshold: this allows discrimination between noise and the underlying image. Ideally, the value should be set greater than the noise level and less than the contrast of the underlying image. Edges of contrast smaller than this threshold will be blurred whereas those of greater contrast will not be. Reducing bt gives less smoothing. bt has a default value = 0.1 * (srcMax - srcMin).
      maskSD - This determines the spatial extent of the smoothing. The mask is basically Gaussian with standard deviation (in image units - e.g. mm.) set by the user. However, for a small, fast, flat response with a 3x3 or 3x3x3 voxel mask, set maskSD to zero. maskSD has a default value equal to the x resolution.
      useMedian - If true, the center pixel is not included in the sums. When the local neighborhood of similar brightness voxels is empty, a local median filter is used. This allows the correction of impulse("salt-and-pepper")noise. If false, this feature is turned off. In this case, the center pixel is included in the sums. When no neighborhood is found, the original intensity of the voxel of interest remains unchanged.
      img25D - Flag, if true, indicates that each slice of the 3D volume should be processed independently. 2D images disregard this flag.

  • Method Details

    • finalize

      public void finalize()
      Prepares this class for destruction.
      Overrides:
      finalize in class AlgorithmBase

    • runAlgorithm

      public void runAlgorithm()
      Starts the program.
      Specified by:
      runAlgorithm in class AlgorithmBase

    • addBorder

      private void addBorder(float[] in, float[] tmpEnImage, int border)
      DOCUMENT ME!
      Parameters:
      in - DOCUMENT ME!
      tmpEnImage - DOCUMENT ME!
      border - DOCUMENT ME!

    • addBorder3D

      private void addBorder3D(float[] data, float[] tmpEnImage, int border)
      This method adds a border to the input data so that borders can be dealt with easily.
      Parameters:
      data - DOCUMENT ME!
      tmpEnImage - DOCUMENT ME!
      border - DOCUMENT ME!

    • calcMinMaxSlice

      private void calcMinMaxSlice(float[] data)
      Calculates the min and max values for the image array, so that the image is displayed properly.
      Parameters:
      data - DOCUMENT ME!

    • median2D

      private float median2D(float[] data, int x, int y)
      Simple 3x3 median filter.
      Parameters:
      data - image data
      x - width index into image data
      y - height index into image data
      Returns:
      median value

    • median3D

      private float median3D(float[] data, int x, int y, int z)
      Simple 3x3x3 median filter.
      Parameters:
      data - image data
      x - width index into image data
      y - height index into image data
      z - image slice index into image data
      Returns:
      median value

    • setupBrightness

      private void setupBrightness()
      The SUSAN brightness function is implemented as a LUT (Look-Up-Table) for speed. In setupBrightness() the LUT is setup. This creates a 257 entry brightness table. The SUSAN function e^-(x^2) is calculated and converted to a byte in the range 0-100, for all possible image brightness differences (including negative ones). Thus brightnessTable[23] is the output for a brightness difference of 23*(srcMax - srcMin)/256. In the smoothing and smoothing3D algorithms this LUT is used as follows: ip = ((i-yMaskSize)*xSize) + j - xMaskSize; ip = (i-1)*x_size + j - 1; ip points to the first image pixel in the circular mask surrounding point (x,y). ip = (k-zMaskSize)*xSize*ySize + (i-yMaskSize)*xSize + j-xMaskSize; ip = (k-1)*xSize*ySize + (i-1)*xSize + j-1; ip points to the first image voxel in the spherical mask surrounding point (x,y,z). brightness is the intensity value at this mask point i*x_size+j points to the center pixel (x,y). k*xSize*ySize+i*xSize+j points to the center voxel (x,y,z). center is the intensity value of this center point tmp = gaussMask[dpt++] * brightnessTable[(int)Math.abs(256*(center-brightness)/(srcMax - srcMin))]; brightnessTable gives the SUSAN brightness function of the absolute value of the difference of the 2 intensity values. dpt gives the current position inside the mask, so the SUSAN brightness function is multiplied by a Gaussian weighting function. In area += tmp, tmp is added to the area and in total += tmp * brightness, the total area*brightness is calculated. After summing over all mask values and performing median replacement if requested, total is divided by area to obtain the filtered intensity value.

    • smoothing2D

      private void smoothing2D(int z)
      Performs 2D and 2.5 nonlinear smoothing.
      Parameters:
      z - The z slice

    • smoothing3D

      private void smoothing3D()
      Performs 3D nonlinear smoothing.