imagefft.java

Digital Image Processing: A Practical Introduction Using Java Nick Efford

    if (!spectral)
      throw new FFTException(NO_DATA);
    if (u >= 0 && u < width && v >= 0 && v < height)
      return data[v*width+u].getMagnitude();
    else
      return 0.0f;
  }


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

  public float getPhase(int u, int v) throws FFTException {
    if (!spectral)
      throw new FFTException(NO_DATA);
    if (u >= 0 && u < width && v >= 0 && v < height)
      return data[v*width+u].getPhase();
    else
      return 0.0f;
  }


  /**
   * Modifies magnitude at any point in the spectrum.
   * @param u horizontal spatial frequency
   * @param v vertical spatial frequency
   * @param mag new magnitiude for specified point
   * @exception FFTException if spectral data are not available.
   */

  public void setMagnitude(int u, int v, float mag) throws FFTException {
    if (!spectral)
      throw new FFTException(NO_DATA);
    if (u >= 0 && u < width && v >= 0 && v < height) {
      int i = v*width+u;
      data[i].setPolar(mag, data[i].getPhase());
    }
  }


  /**
   * Modifies phase at any point in the spectrum.
   * @param u horizontal spatial frequency
   * @param v vertical spatial frequency
   * @param phase new phase for specified point
   * @exception FFTException if spectral data are not available.
   */

  public void setPhase(int u, int v, float phase) throws FFTException {
    if (!spectral)
      throw new FFTException(NO_DATA);
    if (u >= 0 && u < width && v >= 0 && v < height) {
      int i = v*width+u;
      data[i].setPolar(data[i].getMagnitude(), phase);
    }
  }


  /**
   * Performs ideal low pass filtering on the spectrum.
   * @param radius filter radius
   * @exception FFTException if spectral data are not available
   *  or filter radius is invalid.
   */

  public void idealLowPassFilter(double radius) throws FFTException {

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

    if (radius < 0.0 || radius > 1.0)
      throw new FFTException("invalid filter radius");

    int u2 = width/2;
    int v2 = height/2;
    int su, sv, i = 0;
    double r, rmax = Math.min(u2, v2);

    for (int v = 0; v < height; ++v) {
      sv = shift(v, v2) - v2;
      for (int u = 0; u < width; ++u, ++i) {
        su = shift(u, u2) - u2;
        r = Math.sqrt(su*su + sv*sv) / rmax;
        if (r > radius)
          data[i].re = data[i].im = 0.0f;
      }
    }

  }


  /**
   * Performs ideal high pass filtering on the spectrum.
   * @param radius filter radius
   * @exception FFTException if spectral data are not available
   *  or filter radius is invalid.
   */

  public void idealHighPassFilter(double radius) throws FFTException {

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

    if (radius < 0.0 || radius > 1.0)
      throw new FFTException("invalid filter radius");

    int u2 = width/2;
    int v2 = height/2;
    int su, sv, i = 0;
    double r, rmax = Math.min(u2, v2);

    for (int v = 0; v < height; ++v) {
      sv = shift(v, v2) - v2;
      for (int u = 0; u < width; ++u, ++i) {
        su = shift(u, u2) - u2;
        r = Math.sqrt(su*su + sv*sv) / rmax;
        if (r < radius)
          data[i].re = data[i].im = 0.0f;
      }
    }

  }


  /**
   * Performs ideal band pass filtering on the spectrum.
   * @param radius filter radius
   * @param delta band width
   * @exception FFTException if spectral data are not available
   *  or filter parameters are invalid.
   */

  public void idealBandPassFilter(double radius, double delta)
   throws FFTException {

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

    double delta2 = delta/2.0;
    double r1 = radius - delta2;
    double r2 = radius + delta2;
    if (r1 < 0.0 || r2 > 1.0)
      throw new FFTException(INVALID_PARAMS);

    int u2 = width/2;
    int v2 = height/2;
    int su, sv, i = 0;
    double r, rmax = Math.min(u2, v2);

    for (int v = 0; v < height; ++v) {
      sv = shift(v, v2) - v2;
      for (int u = 0; u < width; ++u, ++i) {
        su = shift(u, u2) - u2;
        r = Math.sqrt(su*su + sv*sv) / rmax;
        if (r < r1 || r > r2)
          data[i].re = data[i].im = 0.0f;
      }
    }

  }


  /**
   * Performs ideal band stop filtering on the spectrum.
   * @param radius filter radius
   * @param delta band width
   * @exception FFTException if spectral data are not available
   *  or filter parameters are invalid.
   */

  public void idealBandStopFilter(double radius, double delta)
   throws FFTException {

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

    double delta2 = delta/2.0;
    double r1 = radius - delta2;
    double r2 = radius + delta2;
    if (r1 < 0.0 || r2 > 1.0)
      throw new FFTException(INVALID_PARAMS);

    int u2 = width/2;
    int v2 = height/2;
    int su, sv, i = 0;
    double r, rmax = Math.min(u2, v2);

    for (int v = 0; v < height; ++v) {
      sv = shift(v, v2) - v2;
      for (int u = 0; u < width; ++u, ++i) {
        su = shift(u, u2) - u2;
        r = Math.sqrt(su*su + sv*sv) / rmax;
        if (r >= r1 && r <= r2)
          data[i].re = data[i].im = 0.0f;
      }
    }

  }


  /**
   * Performs order-1 Butterworth low pass filtering of the spectrum.
   * @param radius filter radius
   * @exception FFTException if spectral data are not available
   *  or filter parameters are invalid.
   */

  public void butterworthLowPassFilter(double radius) throws FFTException {
    butterworthLowPassFilter(1, radius);
  }

  /**
   * Performs Butterworth low pass filtering of the spectrum.
   * @param n order of filter
   * @param radius filter radius
   * @exception FFTException if spectral data are not available or
   *  filter parameters are invalid.
   */

  public void butterworthLowPassFilter(int n, double radius)
   throws FFTException {

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

    if (n < 1 || radius <= 0.0 || radius > 1.0)
      throw new FFTException(INVALID_PARAMS);

    int u2 = width/2;
    int v2 = height/2;
    int su, sv, i = 0;
    double mag, r, rmax = Math.min(u2, v2);

    for (int v = 0; v < height; ++v) {
      sv = shift(v, v2) - v2;
      for (int u = 0; u < width; ++u, ++i) {
        su = shift(u, u2) - u2;
        r = Math.sqrt(su*su + sv*sv) / rmax;
        mag = butterworthLowPassFunction(n, radius, r)*data[i].getMagnitude();
        data[i].setPolar(mag, data[i].getPhase());
      }
    }

  }


  /**
   * Performs order-1 Butterworth high pass filtering of the spectrum.
   * @param radius filter radius
   * @exception FFTException if spectral data are not available or
   *  filter parameters are invalid.
   */

  public void butterworthHighPassFilter(double radius) throws FFTException {
    butterworthHighPassFilter(1, radius);
  }

  /**
   * Performs Butterworth high pass filtering of the spectrum.
   * @param n order of filter
   * @param radius filter radius
   * @exception FFTException if spectral data are not available or
   *  filter parameters are invalid.
   */

  public void butterworthHighPassFilter(int n, double radius)
   throws FFTException {

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

    if (n < 1 || radius <= 0.0 || radius > 1.0)
      throw new FFTException(INVALID_PARAMS);

    int u2 = width/2;
    int v2 = height/2;
    int su, sv, i = 0;
    double mag, r, rmax = Math.min(u2, v2);

    for (int v = 0; v < height; ++v) {
      sv = shift(v, v2) - v2;
      for (int u = 0; u < width; ++u, ++i) {
        su = shift(u, u2) - u2;
        r = Math.sqrt(su*su + sv*sv) / rmax;
        mag = butterworthHighPassFunction(n, radius, r)*data[i].getMagnitude();
        data[i].setPolar(mag, data[i].getPhase());
      }
    }

  }


  /**
   * Performs order-1 Butterworth band pass filtering of the spectrum.
   * @param radius filter radius
   * @param delta band width
   * @exception FFTException if spectral data are not available or
   *  filter parameters are invalid.
   */

  public void butterworthBandPassFilter(double radius, double delta)
   throws FFTException {
    butterworthBandPassFilter(1, radius, delta);
  }

  /**
   * Performs Butterworth band pass filtering of the spectrum.
   * @param n order of filter
   * @param radius filter radius
   * @param delta band width
   * @exception FFTException if spectral data are not available or
   *  filter parameters are invalid.
   */

  public void butterworthBandPassFilter(int n, double radius, double delta)
   throws FFTException {

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

    double delta2 = delta/2.0;
    if (n < 1 || radius-delta2 <= 0.0 || radius+delta2 > 1.0)
      throw new FFTException(INVALID_PARAMS);

    int u2 = width/2;
    int v2 = height/2;
    int su, sv, i = 0;
    double mag, r, rmax = Math.min(u2, v2);

    for (int v = 0; v < height; ++v) {
      sv = shift(v, v2) - v2;
      for (int u = 0; u < width; ++u, ++i) {
        su = shift(u, u2) - u2;
        r = Math.sqrt(su*su + sv*sv) / rmax;
        mag = butterworthBandPassFunction(n, radius, delta, r)
         * data[i].getMagnitude();
        data[i].setPolar(mag, data[i].getPhase());
      }
    }

  }


  /**
   * Performs order-1 Butterworth band stop filtering of the spectrum.
   * @param radius filter radius
   * @param delta band width
   * @exception FFTException if spectral data are not available or
   *  filter parameters are invalid.
   */

  public void butterworthBandStopFilter(double radius, double delta)
   throws FFTException {
    butterworthBandStopFilter(1, radius, delta);
  }

  /**
   * Performs Butterworth band stop filtering of the spectrum.
   * @param n order of filter
   * @param radius filter radius
   * @param delta band width
   * @exception FFTException if spectral data are not available or
   *  filter parameters are invalid.
   */

  public void butterworthBandStopFilter(int n, double radius, double delta)
   throws FFTException {

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

    double delta2 = delta/2.0;
    if (n < 1 || radius-delta2<= 0.0 || radius+delta2 > 1.0)
      throw new FFTException(INVALID_PARAMS);

    int u2 = width/2;
    int v2 = height/2;
    int su, sv, i = 0;
    double mag, r, rmax = Math.min(u2, v2);

    for (int v = 0; v < height; ++v) {
      sv = shift(v, v2) - v2;
      for (int u = 0; u < width; ++u, ++i) {
        su = shift(u, u2) - u2;
        r = Math.sqrt(su*su + sv*sv) / rmax;
        mag = butterworthBandStopFunction(n, radius, delta, r)
         * data[i].getMagnitude();
        data[i].setPolar(mag, data[i].getPhase());
      }
    }

  }


  ///////////////////////////// PRIVATE METHODS ////////////////////////////


  /**
   * Computes the power of two for which the corresponding value
   * equals or exceeds the specified integer.
   * @param n integer value
   * @return power of two
   */

  private static int powerOfTwo(int n) {
   int i = 0, m = 1;
   while (m < n) {
     m <<= 1;
     ++i;
   }
   return i;
  }


  /**
   * Reorders an array of data to prepare for an FFT.  The element at
   * index i is swapped with the element at an index given by the
   * bit-reversed value of i.
   * @param data array of complex values
   * @param n number of values
   */

  private static void reorder(Complex[] data, int n) {
    int j = 0, m;
    for (int i = 0; i < n; ++i) {
      if (i > j)
        data[i].swapWith(data[j]);
      m = n >> 1;
      while ((j >= m) && (m >= 2)) {
        j -= m;
        m >>= 1;
      }
      j += m;
    }
  }


  /**
   * Performs a one-dimensional FFT on the specified data.
   * @param data input data, already reordered by bit-reversal
   * @param size number of data elements
   * @param log2n base-2 logarithm of number of elements
   * @param dir direction of transform (1 = forward, -1 = inverse)
   */

  private static void fft(Complex[] data, int size, int log2n, int dir) {

    double angle, wtmp, wpr, wpi, wr, wi, tmpr, tmpi;
    int n = 1, n2;
    for (int k = 0; k < log2n; ++k) {

      n2 = n;
      n <<= 1;
      angle = (-TWO_PI/n) * dir;
      wtmp = Math.sin(0.5*angle);
      wpr = -2.0*wtmp*wtmp;
      wpi = Math.sin(angle);
      wr = 1.0;
      wi = 0.0;

      for (int m = 0; m < n2; ++m) {
        for (int i = m; i < size; i += n) {
          int j = i + n2;
          tmpr = wr*data[j].re - wi*data[j].im;
          tmpi = wi*data[j].re + wr*data[j].im;
          data[j].re = (float) (data[i].re - tmpr);
          data[i].re += (float) tmpr;
          data[j].im = (float) (data[i].im - tmpi);
          data[i].im += (float) tmpi;
        }
        wtmp = wr;
        wr = wtmp*wpr - wi*wpi + wr;
        wi = wi*wpr + wtmp*wpi + wi;
      }

    }

    // Rescale results of inverse transform

    if (dir == -1)
      for (int i = 0; i < size; ++i) {
        data[i].re /= size;
        data[i].im /= size;
      }

  }


  /**
   * Shifts a coordinate relative to a centre point.
   * @param d coordinate
   * @param d2 centre point
   * @return shifted coordinate.
   */

  private static final int shift(int d, int d2) {
    return (d >= d2 ? d-d2 : d+d2);
  }


  /**
   * Collects magnitudes from spectral data, storing them
   * in the specified array.
   * @param mag destination for magnitudes
   * @return maximum magnitude.
   */

  private float calculateMagnitudes(float[] mag) {
    float maximum = 0.0f;
    for (int i = 0; i < data.length; ++i) {
      mag[i] = data[i].getMagnitude();
      if (mag[i] > maximum)
        maximum = mag[i];
    }
    return maximum;
  }


}