imagefft.java

Digital Image Processing: A Practical Introduction Using Java Nick Efford

/***************************************************************************

  ImageFFT.java

  Written by Nick Efford.

  Copyright (c) 2000, Pearson Education Ltd.  All rights reserved.

  THIS SOFTWARE IS PROVIDED BY THE AUTHOR "AS IS" AND ANY EXPRESS OR
  IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
  WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
  ARE DISCLAIMED.  IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE
  LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
  CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
  SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR
  BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
  LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
  NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
  SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

***************************************************************************/


package com.pearsoneduc.ip.op;


import java.awt.image.*;
import com.pearsoneduc.ip.util.Complex;


/**
 * Computes the FFT of an image, and the inverse FFT of its frequency
 * domain representation.  The transformed data can be inspected or
 * filtered before performing the inverse transform.
 *
 * @author Nick Efford
 * @version 1.3 [1999/08/09]
 * @see FFTException
 * @see java.awt.image.BufferedImage
 */

public class ImageFFT {


  ///////////////////////////// CLASS CONSTANTS ////////////////////////////


  public static final int NO_WINDOW = 1;
  public static final int BARTLETT_WINDOW = 2;
  public static final int HAMMING_WINDOW = 3;
  public static final int HANNING_WINDOW = 4;
  private static final String NO_DATA = "no spectral data available";
  private static final String INVALID_PARAMS = "invalid filter parameters";
  private static final double TWO_PI = 2.0*Math.PI;


  /////////////////////////// INSTANCE VARIABLES ///////////////////////////


  /** Complex storage for results of FFT. */
  private Complex[] data;

  /** base-2 logarithm of transform width. */
  private int log2w;

  /** base-2 logarithm of transform height. */
  private int log2h;

  /** Width of transform. */
  private int width;

  /** Height of transform. */
  private int height;

  /** Windowing function applied to image data. */
  private int window;

  /** Indicates whether we have spectral or spatial data. */
  private boolean spectral = false;


  ///////////////////////////// PUBLIC METHODS /////////////////////////////


  /**
   * Computes one half of a radial Bartlett windowing function.
   * @param r distance from centre of data
   * @param rmax maximum distance
   * @return function value.
   */

  public static final double bartlettWindow(double r, double rmax) {
    return 1.0 - Math.min(r, rmax)/rmax;
  }


  /**
   * Computes one half of a radial Hamming windowing function.
   * @param r distance from centre of data
   * @param rmax maximum distance
   * @return function value.
   */

  public static final double hammingWindow(double r, double rmax) {
    double f = (rmax - Math.min(r, rmax)) / rmax;
    return 0.54 - 0.46*Math.cos(f*Math.PI);
  }


  /**
   * Computes one half of a radial Hanning windowing function.
   * @param r distance from centre of data
   * @param rmax maximum distance
   * @return function value.
   */

  public static final double hanningWindow(double r, double rmax) {
    double f = (rmax - Math.min(r, rmax)) / rmax;
    return 0.5 - 0.5*Math.cos(f*Math.PI);
  }


  /**
   * Computes the transfer function for a Butterworth low pass filter.
   * @param n order of filter
   * @param radius filter radius
   * @param r distance from centre of spectrum
   * @return transfer function value.
   */

  public static final double butterworthLowPassFunction(
   int n, double radius, double r) {
    double p = Math.pow(r/radius, 2.0*n);
    return 1.0/(1.0 + p);
  }


  /**
   * Computes the transfer function for a Butterworth high pass filter.
   * @param n order of filter
   * @param radius filter radius
   * @param r distance from centre of spectrum
   * @return transfer function value.
   */

  public static final double butterworthHighPassFunction(
   int n, double radius, double r) {
    try {
      double p = Math.pow(radius/r, 2.0*n);
      return 1.0/(1.0 + p);
    }
    catch (ArithmeticException e) {
      return 0.0;
    }
  }


  /**
   * Computes the transfer function for a Butterworth band pass filter.
   * @param n order of filter
   * @param radius filter radius
   * @param delta band width
   * @param r distance from centre of spectrum
   * @return transfer function value.
   */

  public static final double butterworthBandPassFunction(
   int n, double radius, double delta, double r) {
    return 1.0-butterworthBandStopFunction(n, radius, delta, r);
  }


  /**
   * Computes the transfer function of a Butterworth band stop filter.
   * @param n order of filter
   * @param radius filter radius
   * @param delta band width
   * @param r distance from centre of spectrum
   * @return transfer function value.
   */

  public static final double butterworthBandStopFunction(
   int n, double radius, double delta, double r) {
    try {
      double p = Math.pow(delta*radius/(r*r - radius*radius), 2.0*n);
      return 1.0/(1.0 + p);
    }
    catch (ArithmeticException e) {
      return 0.0;
    }
  }


  /**
   * Creates an ImageFFT for the specified image.  There will be
   * no windowing of image data.
   * @param image input image
   * @exception FFTException if the image is not 8-bit greyscale.
   */

  public ImageFFT(BufferedImage image) throws FFTException {
    this(image, NO_WINDOW);
  }

  /**
   * Creates an ImageFFT for the specified image, applying the
   * specified windowing function to the data.
   * @param image input image
   * @param win windowing function
   * @exception FFTException if the image is not 8-bit greyscale.
   */

  public ImageFFT(BufferedImage image, int win) throws FFTException {

    if (image.getType() != BufferedImage.TYPE_BYTE_GRAY)
      throw new FFTException("image must be 8-bit greyscale");

    // Compute dimensions, allowing for zero padding

    log2w = powerOfTwo(image.getWidth());
    log2h = powerOfTwo(image.getHeight());
    width = 1 << log2w;
    height = 1 << log2h;
    window = win;

    // Allocate storage for results of FFT

    data = new Complex[width*height];
    for (int i = 0; i < data.length; ++i)
      data[i] = new Complex();

    Raster raster = image.getRaster();
    double xc = image.getWidth()/2.0, yc = image.getHeight()/2.0;
    double r, rmax = Math.min(xc, yc);
    switch (window) {

      case HAMMING_WINDOW:
        for (int y = 0; y < image.getHeight(); ++y)
          for (int x = 0; x < image.getWidth(); ++x) {
            r = Math.sqrt((x-xc)*(x-xc) + (y-yc)*(y-yc));
            data[y*width+x].re =
             (float) (hammingWindow(r, rmax)*raster.getSample(x, y, 0));
          }
        break;

      case HANNING_WINDOW:
        for (int y = 0; y < image.getHeight(); ++y)
          for (int x = 0; x < image.getWidth(); ++x) {
            r = Math.sqrt((x-xc)*(x-xc) + (y-yc)*(y-yc));
            data[y*width+x].re =
             (float) (hanningWindow(r, rmax)*raster.getSample(x, y, 0));
          }
        break;

      case BARTLETT_WINDOW:
        for (int y = 0; y < image.getHeight(); ++y)
          for (int x = 0; x < image.getWidth(); ++x) {
            r = Math.sqrt((x-xc)*(x-xc) + (y-yc)*(y-yc));
            data[y*width+x].re =
             (float) (bartlettWindow(r, rmax)*raster.getSample(x, y, 0));
          }
        break;

      default:  // NO_WINDOW
        for (int y = 0; y < image.getHeight(); ++y)
          for (int x = 0; x < image.getWidth(); ++x)
            data[y*width+x].re = raster.getSample(x, y, 0);
        break;

    }

  }


  /**
   * @return width of FFT.
   */

  public int getWidth() {
    return width;
  }


  /**
   * @return height of FFT.
   */

  public int getHeight() {
    return height;
  }


  /**
   * @return current windowing function.
   */

  public int getWindow() {
    return window;
  }


  /**
   * @return true if data are spectral, false if data are spatial.
   */

  public boolean isSpectral() {
    return spectral;
  }


  /**
   * @return information string for an ImageFFT object.
   */

  public String toString() {
    String s = new String(getClass().getName() + ": " + width + "x" + height +
     (spectral ? ", frequency domain" : ", spatial domain"));
    return s;
  }


  /**
   * Transforms data via a forward or inverse FFT, as appropriate.
   * An inverse transform is computed if the previous transform was
   * in the forward direction; otherwise, the forward transform
   * is computed.
   */

  public void transform() {

    int x, y, i;
    Complex[] row = new Complex[width];
    for (x = 0; x < width; ++x)
      row[x] = new Complex();
    Complex[] column = new Complex[height];
    for (y = 0; y < height; ++y)
      column[y] = new Complex();

    int direction;
    if (spectral)
      direction = -1;   // inverse transform
    else
      direction = 1;    // forward transform

    // Perform FFT on each row

    for (y = 0; y < height; ++y) {
      for (i = y*width, x = 0; x < width; ++x, ++i) {
        row[x].re = data[i].re;
        row[x].im = data[i].im;
      }
      reorder(row, width);
      fft(row, width, log2w, direction);
      for (i = y*width, x = 0; x < width; ++x, ++i) {
        data[i].re = row[x].re;
        data[i].im = row[x].im;
      }
    }

    // Perform FFT on each column

    for (x = 0; x < width; ++x) {
      for (i = x, y = 0; y < height; ++y, i += width) {
        column[y].re = data[i].re;
        column[y].im = data[i].im;
      }
      reorder(column, height);
      fft(column, height, log2h, direction);
      for (i = x, y = 0; y < height; ++y, i += width) {
        data[i].re = column[y].re;
        data[i].im = column[y].im;
      }
    }

    if (spectral)
      spectral = false;
    else
      spectral = true;

  }


  /**
   * Converts stored data into an image.
   * @param image destination image, or null
   * @return FFT data as an image.
   * @exception FFTException if the data are in spectral form; an
   *  image can be created only from data in the spatial domain.
   */

  public BufferedImage toImage(BufferedImage image) throws FFTException {
    return toImage(image, 0);
  }

  /**
   * Converts stored data into an image.
   * @param image destination image, or null
   * @param bias constant value added to data
   * @return FFT data as an image.
   * @exception FFTException if the data are in spectral form; an
   *  image can be created only from data in the spatial domain.
   */

  public BufferedImage toImage(BufferedImage image, int bias)
   throws FFTException {

    if (spectral)
      throw new FFTException("cannot transfer spectral data to an image");

    if (image == null)
      image = new BufferedImage(width, height, BufferedImage.TYPE_BYTE_GRAY);
    WritableRaster raster = image.getRaster();

    int w = Math.min(image.getWidth(), width);
    int h = Math.min(image.getHeight(), height);

    // If destination image is bigger, zero it

    if (w < image.getWidth() || h < image.getHeight())
      for (int y = 0; y < image.getHeight(); ++y)
        for (int x = 0; x < image.getWidth(); ++x)
          raster.setSample(x, y, 0, 0);

    // Copy real component of data to destination image

    int i = 0, value;
    for (int y = 0; y < height; ++y)
      for (int x = 0; x < width; ++x, ++i) {
        value = Math.max(0, Math.min(255, bias + Math.round(data[i].re)));
        raster.setSample(x, y, 0, value);
      }

    return image;

  }


  /**
   * Returns the amplitude spectrum of an image, as another image.
   * The data are shifted such that the DC component is at the image
   * centre, and scaled logarithmically so that low-amplitude
   * detail is visible.
   * @return shifted spectrum, as an image.
   * @exception FFTException if spectral data are not available
   * (e.g. because last transform was in the inverse direction).
   */

  public BufferedImage getSpectrum() throws FFTException {

    if (!spectral)
      throw new FFTException(NO_DATA);

    // Collect magnitudes and find maximum

    float[] magData = new float[width*height];
    float maximum = calculateMagnitudes(magData);
    BufferedImage image =
     new BufferedImage(width, height, BufferedImage.TYPE_BYTE_GRAY);
    WritableRaster raster = image.getRaster();

    // Shift, rescale and copy to image

    double scale = 255.0 / Math.log(maximum + 1.0);
    int x2 = width/2, y2 = height/2;
    int sx, sy, value;
    for (int y = 0; y < height; ++y) {
      sy = shift(y, y2);
      for (int x = 0; x < width; ++x) {
        sx = shift(x, x2);
        value = (int) Math.round(scale*Math.log(magData[sy*width+sx]+1.0));
        raster.setSample(x, y, 0, value);
      }
    }

    return image;

  }


  /**
   * Returns the amplitude spectrum of an image, as another image.
   * The data are unshifted and are scaled logarithmically so that
   * low-amplitude detail is visible.
   * @return unshifted spectrum, as an image.
   * @exception FFTException if spectral data are not available
   * e.g. because last transform was in the inverse direction).
   */

  public BufferedImage getUnshiftedSpectrum() throws FFTException {

    if (!spectral)
      throw new FFTException(NO_DATA);

    // Collect magnitudes and find maximum

    float[] magData = new float[width*height];
    float maximum = calculateMagnitudes(magData);
    BufferedImage image =
     new BufferedImage(width, height, BufferedImage.TYPE_BYTE_GRAY);
    WritableRaster raster = image.getRaster();

    // Rescale and copy to image

    double scale = 255.0 / Math.log(maximum + 1.0);
    int i = 0, value;
    for (int y = 0; y < height; ++y)
      for (int x = 0; x < width; ++x, ++i) {
        value = (int) Math.round(scale*Math.log(magData[i]+1.0));
        raster.setSample(x, y, 0, value);
      }

    return image;

  }


  /**
   * Computes magnitude for any point in the spectrum.
   * @param u horizontal spatial frequency
   * @param v vertical spatial frequency
   * @return magnitude at the specified point, or zero if point
   * does not exist.
   * @exception FFTException if spectral data are not available.
   */

  public float getMagnitude(int u, int v) throws FFTException {