Function Reference

Phase congruency of an image using monogenic filters.

This code is considerably faster than phasecong3() but you may prefer the output from phasecong3()'s oriented filters.

There are potentially many arguments, here is the full usage:

  (PC, or, ft, T) =
          phasecongmono(img; nscale, minwavelength, mult,
                        sigmaonf, k, cutoff, g, deviationgain, noisemethod)

However, apart from the image, all parameters have defaults and the
usage can be as simple as:

   (PC,) = phasecongmono(img)   # Use (PC,) so that PC is not a tuple of all
                                # the returned values

More typically you will pass the image followed by a series of keyword
arguments that you wish to set, leaving the remaining parameters set to
their defaults, for example:

   (PC,) = phasecongmono(img, nscale = 5, minwavelength = 3, k = 2.5)

Keyword arguments:
             Default values      Description

   nscale           4    - Number of wavelet scales, try values 3-6
                           A lower value will reveal more fine scale
                           features. A larger value will highlight 'major'
                           features.
   minwavelength    3    - Wavelength of smallest scale filter.
   mult             2.1  - Scaling factor between successive filters.
   sigmaonf         0.55 - Ratio of the standard deviation of the Gaussian
                           describing the log Gabor filter's transfer function
                           in the frequency domain to the filter center frequency.
   k                3.0  - No of standard deviations of the noise energy beyond
                           the mean at which we set the noise threshold point.
                           You may want to vary this up to a value of 10 or
                           20 for noisy images
   cutoff           0.5  - The fractional measure of frequency spread
                           below which phase congruency values get penalized.
   g                10   - Controls the sharpness of the transition in
                           the sigmoid function used to weight phase
                           congruency for frequency spread.
   deviationgain    1.5  - Amplification to apply to the calculated phase
                           deviation result. Increasing this sharpens the
                           edge responses, but can also attenuate their
                           magnitude if the gain is too large.  Sensible
                           values to use lie in the range 1-2.
   noisemethod      -1   - Parameter specifies method used to determine
                           noise statistics.
                             -1 use median of smallest scale filter responses
                             -2 use mode of smallest scale filter responses
                              0+ use noiseMethod value as the fixed noise threshold
                           A value of 0 will turn off all noise compensation.

Returned values:
   PC         - Phase congruency indicating edge significance
   or         - Orientation image in radians -pi/2 to pi/2,  +ve anticlockwise.
                0 corresponds to a vertical edge, pi/2 is horizontal.
   ft         - Local weighted mean phase angle at every point in the
                image.  A value of pi/2 corresponds to a bright line, 0
                corresponds to a step and -pi/2 is a dark line.
   T          - Calculated noise threshold (can be useful for
                diagnosing noise characteristics of images).  Once you know
                this you can then specify fixed thresholds and save some
                computation time.

The convolutions are done via the FFT. Many of the parameters relate to the specification of the filters in the frequency plane. The values do not seem to be very critical and the defaults are usually fine. You may want to experiment with the values of nscales and k, the noise compensation factor.

Typical sequence of operations to obtain an edge image:

 > (PC, or) = phasecongmono(img)
 > nm = nonmaxsup(PC, or, 1.5)   # nonmaxima suppression
 > bw = hysthresh(nm, 0.1, 0.3)  # hysteresis thresholding 0.1 - 0.3

Notes on filter settings to obtain even coverage of the spectrum
sigmaonf       .85   mult 1.3
sigmaonf       .75   mult 1.6     (filter bandwidth ~1 octave)
sigmaonf       .65   mult 2.1
sigmaonf       .55   mult 3       (filter bandwidth ~2 octaves)

Note that better results are generally achieved using the large bandwidth filters. I typically use a sigmaOnf value of 0.55 or even smaller.

See also: phasecong3, phasesymmono, gaborconvolve, filtergrid

Phase Preserving Dynamic Range Compression

Generates a series of dynamic range compressed images at different scales. This function is designed to reveal subtle features within high dynamic range images such as aeromagnetic and other potential field grids. Often this kind of data is presented using histogram equalisation in conjunction with a rainbow colourmap. A problem with histogram equalisation is that the contrast amplification of a feature depends on how commonly its data value occurs, rather than on the amplitude of the feature itself.

Phase Preserving Dynamic Range Compression allows subtle features to be revealed without these distortions. Perceptually important phase information is preserved and the contrast amplification of anomalies in the signal is purely a function of their amplitude. It operates as follows: first a highpass filter is applied to the data, this controls the desired scale of analysis. The 2D analytic signal of the data is then computed to obtain local phase and amplitude at each point in the image. The amplitude is attenuated by adding 1 and then taking its logarithm, the signal is then reconstructed using the original phase values.

Usage: dimg = ppdrc(img, wavelength; clip, n)

Arguments:     img - Image to be processed. A 2D array of Real or Gray elements.
        wavelength - Scalar value, or Vector, of wavelengths, in pixels, of
                     the cut-in frequencies to be used when forming the highpass
                     versions of the image.  Try a range of values starting
                     with, say, a wavelength corresponding to half the size
                     of the image and working down to something like 50
                     grid units.
Keyword arguments:
              clip - Percentage of output image histogram to clip.  Only a
                     very small value should be used, say 0.01 or 0.02, but
                     this can be beneficial.  Defaults to 0.01%
                 n - Order of the Butterworth high pass filter.  Defaults
                     to 2

Returns:      dimg - Array of the dynamic range reduced images.  If only
                     one wavelength is specified the image is returned
                     directly, and not as a one element array of image arrays.

Important: Scaling of the image affects the results. If your image has values of order 1 or less it is useful to scale the image up a few orders of magnitude. The reason is that when the frequency amplitudes are attenuated we add one before taking the log to avoid obtaining negative results for values less than one. Thus if v is small log(1 + v) will not be a good approximation to log(v). However, if you scale the image by say, 1000 then log(1 + 1000*v) will be a reasonable approximation to log(1000*v).

When specifying the array wavelength it is suggested that you use wavelengths that increase in a geometric series. You can use the function geoseries() to conveniently do this

Example using geoseries() to generate a set of wavelengths that increase geometrically in 10 steps from 50 to 800.

   dimg = ppdrc(img, geoseries((50 800), 10))

See also: highpassmonogenic, geoseries

Compute phase and amplitude in highpass images via monogenic filters.

Usage: (phase, orient, E) = highpassmonogenic(img, maxwavelength, n)

Arguments:           img - Image to be processed.  A 2D array of Real or Gray elements.
           maxwavelength - Wavelength(s) in pixels of the  cut-in frequency(ies)
                           of the Butterworth highpass filter.
                       n - The order of the Butterworth filter. This is an
                           integer >= 1.  The higher the value the sharper
                           the cutoff.

Returns:           phase - The local phase. Values are between -pi/2 and pi/2
                  orient - The local orientation. Values between -pi and pi.
                           Note that where the local phase is close to
                           +-pi/2 the orientation will be poorly defined.
                       E - Local energy, or amplitude, of the signal.

Note that maxwavelength can be an array in which case the outputs will be an array of output images of length nscales, where nscales = length(maxwavelength).

See also: bandpassmonogenic, ppdrc, monofilt

Compute phase and amplitude in bandpass images via monogenic filters.

Usage: (phase, orient, E) = bandpassmonogenic(img, minwavelength, maxwavelength, n)

Arguments:           img - Image to be processed.  A 2D array of Real or Gray elements.
           minwavelength - } Wavelength(s) in pixels of the cut-in and cut-out frequency(ies)
           maxwavelength - } of the Butterworth bandpass filter(s).
                       n - The order of the Butterworth filter. This is an
                           integer >= 1.  The higher the value the sharper
                           the cutoff.

Returns:           phase - The local phase. Values are between -pi/2 and pi/2
                  orient - The local orientation. Values between -pi and pi.
                           Note that where the local phase is close to
                           +-pi/2 the orientation will be poorly defined.
                       E - Local energy, or amplitude, of the signal.

Note that minwavelength and maxwavelength can be (equal length) arrays in which case the outputs will be an array of output images of length nscales, where nscales = length(maxwavelength).

See also: highpassmonogenic, ppdrc, monofilt

Phase symmetry of an image using monogenic filters.

This function calculates the phase symmetry of points in an image. This is a contrast invariant measure of symmetry. This function can be used as a line and blob detector. The greyscale polarity of the lines that you want to find can be specified.

This code is considerably faster than phasesym() but you may prefer the output from phasesym()'s oriented filters.

There are potentially many arguments, here is the full usage:

  (phSym, symmetryEnergy, T) =
               phasesymmono(img; nscale, minwaveLength, mult,
                            sigmaonf, k, polarity, noisemethod)

However, apart from the image, all parameters have defaults and the usage can be as simple as:

   (phSym,) = phasesymmono(img)

Keyword arguments:
             Default values      Description

   nscale           5    - Number of wavelet scales, try values 3-6
   minwaveLength    3    - Wavelength of smallest scale filter.
   mult             2.1  - Scaling factor between successive filters.
   sigmaonf         0.55 - Ratio of the standard deviation of the Gaussian
                           describing the log Gabor filter's transfer function
                           in the frequency domain to the filter center frequency.
   k                2.0  - No of standard deviations of the noise energy beyond
                           the mean at which we set the noise threshold point.
                           You may want to vary this up to a value of 10 or
                           20 for noisy images
   polarity         0    - Controls 'polarity' of symmetry features to find.
                            1 - just return 'bright' points
                           -1 - just return 'dark' points
                            0 - return bright and dark points.
   noisemethod      -1   - Parameter specifies method used to determine
                           noise statistics.
                             -1 use median of smallest scale filter responses
                             -2 use mode of smallest scale filter responses
                              0+ use noiseMethod value as the fixed noise threshold
                           A value of 0 will turn off all noise compensation.

Return values:
   phSym                 - Phase symmetry image (values between 0 and 1).
   symmetryEnergy        - Un-normalised raw symmetry energy which may be
                           more to your liking.
   T                     - Calculated noise threshold (can be useful for
                           diagnosing noise characteristics of images)

The convolutions are done via the FFT. Many of the parameters relate to the specification of the filters in the frequency plane. The values do not seem to be very critical and the defaults are usually fine. You may want to experiment with the values of nscales and k, the noise compensation factor.

Notes on filter settings to obtain even coverage of the spectrum

sigmaonf       .85   mult 1.3
sigmaonf       .75   mult 1.6     (filter bandwidth ~1 octave)
sigmaonf       .65   mult 2.1
sigmaonf       .55   mult 3       (filter bandwidth ~2 octaves)

See Also: phasesym, phasecongmono

Apply monogenic filters to an image to obtain 2D analytic signal.

This is an implementation of Felsberg's monogenic filters

Usage: (f, h1f, h2f, A, theta, psi) =
            monofilt(img, nscale, minWaveLength, mult, sigmaOnf, orientWrap)
                             3         4           2     0.65    true/false
Arguments:
The convolutions are done via the FFT.  Many of the parameters relate
to the specification of the filters in the frequency plane.

  Variable       Suggested   Description
  name           value
 ----------------------------------------------------------
   img                       Image to be convolved. An Array of Real or Gray.
   nscale          = 3       Number of filter scales.
   minWaveLength   = 4       Wavelength of smallest scale filter.
   mult            = 2       Scaling factor between successive filters.
   sigmaonf        = 0.65    Ratio of the standard deviation of the
                             Gaussian describing the log Gabor filter's
                             transfer function in the frequency domain
                             to the filter center frequency.
   orientWrap       false    Optional Boolean flag  to turn on/off
                             'wrapping' of orientation data from a
                             range of -pi .. pi to the range 0 .. pi.
                             This affects the interpretation of the
                             phase angle - see note below. Defaults to false.
Returns:
       f  - vector of bandpass filter responses with respect to scale.
     h1f  - vector of bandpass h1 filter responses wrt scale.
     h2f  - vector of bandpass h2 filter responses.
       A  - vector of monogenic energy responses.
   theta  - vector of phase orientation responses.
     psi  - vector of phase angle responses.

If orientWrap is true theta will be returned in the range 0 .. pi

Experimentation with sigmaonf can be useful depending on your application. I have found values as low as 0.2 (a filter with a very large bandwidth) to be useful on some occasions.

See also: gaborconvolve

Convolve an image with a bank of log-Gabor filters.

Usage: (EO, BP) = gaborconvolve(img,  nscale, norient, minWaveLength, mult,
                                 sigmaOnf, dThetaOnSigma, Lnorm)

Arguments:
The convolutions are done via the FFT.  Many of the parameters relate
to the specification of the filters in the frequency plane.

  Variable       Suggested   Description
  name           value
 ----------------------------------------------------------
   img                       Image to be convolved.
   nscale          = 4       Number of wavelet scales.
   norient         = 6       Number of filter orientations.
   minWaveLength   = 3       Wavelength of smallest scale filter.
   mult            = 1.7     Scaling factor between successive filters.
   sigmaOnf        = 0.65    Ratio of the standard deviation of the
                             Gaussian describing the log Gabor filter's
                             transfer function in the frequency domain
                             to the filter center frequency.
   dThetaOnSigma   = 1.3     Ratio of angular interval between filter
                             orientations and the standard deviation of
                             the angular Gaussian function used to
                             construct filters in the freq. plane.
   Lnorm            0        Optional integer indicating what norm the
                             filters should be normalized to.  A value of 1
                             will produce filters with the same L1 norm, 2
                             will produce filters with matching L2
                             norm. the default value of 0 results in no
                             normalization (the filters have unit height
                             Gaussian transfer functions on a log frequency
                             scale)
Returns:

  EO - 2D array of arrays of complex valued convolution results
       EO[s,o] = convolution result for scale s and orientation o.
       The real part is the result of convolving with the even
       symmetric filter, the imaginary part is the result from
       convolution with the odd symmetric filter.

       Hence:
       abs.(EO[s,o]) returns the magnitude of the convolution over the
                    image at scale s and orientation o.
       angle.(EO[s,o]) returns the phase angles.

  BP - Array of bandpass images corresponding to each scale s.

Notes on filter settings to obtain even coverage of the spectrum energy

dThetaOnSigma 1.2 - 1.3
sigmaOnf  .90   mult 1.15
sigmaOnf  .85   mult 1.2
sigmaOnf  .75   mult 1.4       (bandwidth ~1 octave)
sigmaOnf  .65   mult 1.7
sigmaOnf  .55   mult 2.2       (bandwidth ~2 octaves)

The determination of mult given sigmaOnf is entirely empirical. What I do is plot out the sum of the squared filter amplitudes in the frequency domain and see how even the coverage of the spectrum is. If there are concentric 'gaps' in the spectrum one needs to reduce mult and/or reduce sigmaOnf (which increases filter bandwidth)

If there are 'gaps' radiating outwards then one needs to reduce dthetaOnSigma (increasing angular bandwidth of the filters)

Computes edge and corner phase congruency in an image via log-Gabor filters.

There are potentially many arguments, here is the full usage:

  (M, m, or, ft, EO, T) = phasecong3(img; nscale, norient, minwavelength,
                          mult, sigmaonf, k, cutoff, g, noisemethod)

However, apart from the image, all parameters have defaults and the usage can be as simple as:

    (M,) = phasecong3(img)

Keyword Arguments:
             Default values      Description

   nscale           4    - Number of wavelet scales, try values 3-6
   norient          6    - Number of filter orientations.
   minwavelength    3    - Wavelength of smallest scale filter.
   mult             2.1  - Scaling factor between successive filters.
   sigmaonf         0.55 - Ratio of the standard deviation of the Gaussian
                           describing the log Gabor filter's transfer function
                           in the frequency domain to the filter center frequency.
   k                2.0  - No of standard deviations of the noise energy beyond
                           the mean at which we set the noise threshold point.
                           You may want to vary this up to a value of 10 or
                           20 for noisy images
   cutoff           0.5  - The fractional measure of frequency spread
                           below which phase congruency values get penalized.
   g                10   - Controls the sharpness of the transition in
                           the sigmoid function used to weight phase
                           congruency for frequency spread.
   noisemethod      -1   - Parameter specifies method used to determine
                           noise statistics.
                             -1 use median of smallest scale filter responses
                             -2 use mode of smallest scale filter responses
                              0+ use noisemethod value as the fixed noise threshold

Returned values:
   M          - Maximum moment of phase congruency covariance.
                This is used as a indicator of edge strength.
   m          - Minimum moment of phase congruency covariance.
                This is used as a indicator of corner strength.
   or         - Orientation image in radians -pi/2 to pi/2,  +ve anticlockwise.
                0 corresponds to a vertical edge, pi/2 is horizontal.
   ft         - Local weighted mean phase angle at every point in the
                image.  A value of pi/2 corresponds to a bright line, 0
                corresponds to a step and -pi/2 is a dark line.
   EO         - A 2D array of complex valued convolution results for each scale
                and orientation
   T          - Calculated noise threshold (can be useful for
                diagnosing noise characteristics of images).  Once you know
                this you can then specify fixed thresholds and save some
                computation time.

EO[s,o] = convolution result for scale s and orientation o. The real part is the result of convolving with the even symmetric filter, the imaginary part is the result from convolution with the odd symmetric filter.

Hence: abs.(EO[s,o]) returns the magnitude of the convolution over the image at scale s and orientation o, angle.(EO[s,o]) returns the phase angles.

The convolutions are done via the FFT. Many of the parameters relate to the specification of the filters in the frequency plane. The values do not seem to be very critical and the defaults are usually fine. You may want to experiment with the values of nscales and k, the noise compensation factor.

Some filter parameters to obtain even coverage of the spectrum

sigmaonf       .85   mult 1.3
sigmaonf       .75   mult 1.6     (filter bandwidth ~1 octave)
sigmaonf       .65   mult 2.1
sigmaonf       .55   mult 3       (filter bandwidth ~2 octaves)

See also: phasesym, gaborconvolve

Compute phase symmetry on an image via log-Gabor filters.

This function calculates the phase symmetry of points in an image. This is a contrast invariant measure of symmetry. This function can be used as a line and blob detector. The greyscale polarity of the lines that you want to find can be specified.

Usage:   (phSym, orientation, totalEnergy, T) =
                phasesym(img; nscale = 5, norient = 6, minwavelength = 3, mult = 2.1,
                         sigmaonf = 0.55, k = 2, polarity = 0, noisemethod = -1)

However, apart from the image, all parameters have defaults and the
usage can be as simple as:

    (phSym,) = phasesym(img)

Argument:
                    img  - Image to be processed. 2D Array of Real or Gray

Keyword Arguments:
              Default values      Description
    nscale           5    - Number of wavelet scales, try values 3-6
    norient          6    - Number of filter orientations.
    minwavelength    3    - Wavelength of smallest scale filter.
    mult             2.1  - Scaling factor between successive filters.
    sigmaonf         0.55 - Ratio of the standard deviation of the Gaussian
                            describing the log Gabor filter's transfer function
                            in the frequency domain to the filter center frequency.
    k                2.0  - No of standard deviations of the noise energy beyond
                            the mean at which we set the noise threshold point.
                            You may want to vary this up to a value of 10 or
                            20 for noisy images
    polarity         0    - Controls 'polarity' of symmetry features to find.
                             1 - just return 'bright' points
                            -1 - just return 'dark' points
                             0 - return bright and dark points.
    noisemethod      -1   - Parameter specifies method used to determine
                            noise statistics.
                              -1 use median of smallest scale filter responses
                              -2 use mode of smallest scale filter responses
                               0+ use noiseMethod value as the fixed noise threshold.

Return values:
    phSym                 - Phase symmetry image (values between 0 and 1).
    orientation           - Orientation image. Orientation in which local
                            symmetry energy is a maximum, in radians
                            (-pi/2 - pi/2), angles positive anti-clockwise. Note
                            the orientation info is quantized by the number
                            of orientations
    totalEnergy           - Un-normalised raw symmetry energy which may be
                            more to your liking.
    T                     - Calculated noise threshold (can be useful for
                            diagnosing noise characteristics of images).  Once you know
                            this you can then specify fixed thresholds and save some
                            computation time.

The convolutions are done via the FFT. Many of the parameters relate to the specification of the filters in the frequency plane. The values do not seem to be very critical and the defaults are usually fine. You may want to experiment with the values of nscales and k, the noise compensation factor.

Notes on filter settings to obtain even coverage of the spectrum

sigmaonf       .85   mult 1.3
sigmaonf       .75   mult 1.6     (filter bandwidth ~1 octave)
sigmaonf       .65   mult 2.1
sigmaonf       .55   mult 3       (filter bandwidth ~2 octaves)

See also: phasesymmono, phasecong3

Phase preserving wavelet image denoising.

Usage: cleanimage = ppdenoise(img,  nscale = 5, norient = 6,
                              mult = 2.5, minwavelength = 2, sigmaonf = 0.55,
                              dthetaonsigma = 1.0, k = 3, softness = 1.0)
Argument:
          img - Image to be processed (greyscale)

Keyword arguments:
        nscale - No of filter scales to use (5-7) - the more scales used
                 the more low frequencies are covered.
       norient - No of orientations to use (6)
          mult - Multiplying factor between successive scales  (2.5-3)
 minwavelength - Wavelength of smallest scale filter (2)
      sigmaonf - Ratio of the standard deviation of the Gaussian
                 describing the log Gabor filter's transfer function
                 in the frequency domain to the filter center frequency (0.55)
 dthetaonsigma - Ratio of angular interval between filter orientations
                 and the standard deviation of the angular Gaussian (1)
                 function used to construct filters in the freq. plane.
             k - No of standard deviations of noise to reject 2-3
      softness - Degree of soft thresholding (0-hard to 1-soft)

The convolutions are done via the FFT. Many of the parameters relate to the specification of the filters in the frequency plane. Most arguments do not need to be changed from the defaults and are mostly not that critical. The main parameter that you may wish to play with is k, the number of standard deviations of noise to reject.

Generate grids for constructing frequency domain filters.

Usage:  (f, fx, fy) = filtergrids(rows, cols)
        (f, fx, fy) = filtergrids((rows, cols))

Arguments:  rows, cols - Size of image/filter

Returns:             f - Grid of size (rows, cols) containing frequency
                         values from 0 to 0.5,  where f =
                         sqrt(fx^2 + fy^2). The grid is quadrant
                         shifted so that 0 frequency is at f[1,1]

                fx, fy - Grids containing normalised frequency values
                         ranging from -0.5 to 0.5 in x and y directions
                         respectively. fx and fy are quadrant shifted.

See also: filtergrid where you are only needing radius

Generate grid for constructing frequency domain filters.

Usage:  f = filtergrid(rows, cols)
        f = filtergrid((rows, cols))

Arguments:  rows, cols - Size of image/filter

Returns:             f - Grid of size (rows, cols) containing normalised
                         frequency values from 0 to 0.5.  Grid is quadrant
                         shifted so that 0 frequency is at f[1,1]

Used by phasecongmono, phasecong3, etc etc

See also: filtergrids if you also want normalized frequency grids in the x and y directions as well.

Generate monogenic filter grids.

Usage: (H1, H2, f) = monogenicfilters(rows, cols)
       (H1, H2, f) = monogenicfilters((rows, cols))

Arguments:  rows,cols - Size of filters to generate

Returns: H1, H2 - The two monogenic filters.
              f - Frequency grid corresponding to the filters.

where:
       H1 = i*fx./f
       H2 = i*fy./f

Note that H1, H2, and f and quadrant shifted to that the 0 frequency value is at coordinate [1,1].

See also: packedmonogenicfilters

Monogenic filter where both filters are packed in the one Complex grid.

Usage: (H, f) = packedmonogenicfilters(rows, cols)
       (H, f) = packedmonogenicfilters((rows, cols))

Arguments:  rows,cols - Size of filters to generate

Returns:      H - The two monogenic filters packed into the
                  one Complex64 grid.
              f - Frequency grid corresponding to the filter.

The two monogenic filters are defined as

       H1 = i*fx./f
       H2 = i*fy./f

However the two filters can be packed together as a complex valued matrix, one in the real part and one in the imaginary part. Do this by multiplying H2 by i and then adding it to H1. When the convolution is performed via the fft the real part of the result will correspond to the convolution with H1 and the imaginary part with H2. This allows the two convolutions to be done as one in the frequency domain, saving time and memory.

Note that H and f and quadrant shifted to that the 0 frequency value is at coordinate [1,1].

See also: monogenicfilters

Construct a low-pass Butterworth filter.

Usage: f = lowpassfilter(sze, cutoff, n)

where: sze    is a two element tuple specifying the size of filter
              to construct (rows, cols).
       cutoff is the cutoff frequency of the filter 0 - 0.5
       n      is the order of the filter, the higher n is the sharper
              the transition is. (n must be an integer >= 1).
              Note that n is doubled so that it is always an even integer.

                      1
      f =    --------------------
                              2n
              1.0 + (w/cutoff)

The frequency origin of the returned filter is at the corners.

See also: highpassfilter, highboostfilter, bandpassfilter

Construct a band-pass Butterworth filter.

Usage: f = bandpassfilter(sze, cutin, cutoff, n)

Arguments:
             sze - A 2 element tuple specifying the size of filter
                   to construct (rows, cols).
   cutin, cutoff - The frequencies defining the band pass 0 - 0.5
               n - The order of the filter, the higher n is the sharper
                   the transition is. (n must be an integer >= 1).
Returns:
               f - Frequency domain filter of size==sze, the frequency
                   origin is at the corners.

See also: lowpassfilter, highpassfilter, highboostfilter

Construct a high-boost Butterworth filter.

Usage: f = highboostfilter(sze, cutoff, n, boost)

Arguments:
         sze - A 2 element tuple specifying the size of filter
               to construct (rows, cols).
      cutoff - The cutoff frequency of the filter 0 - 0.5
           n - The order of the filter, the higher n is the sharper
               the transition is. (n must be an integer >= 1).
       boost - The ratio that high frequency values are boosted
               relative to the low frequency values.  If boost is less
               than one then a 'lowboost' filter is generated
Returns:
           f - Frequency domain filter of size==sze, the frequency
               origin is at the corners.

See also: lowpassfilter, highpassfilter, bandpassfilter

Construct a high-pass Butterworth filter.

Usage: f = highpassfilter(sze, cutoff, n)

         sze - A 2 element tuple specifying the size of filter
               to construct (rows, cols).
      cutoff - The cutoff frequency of the filter 0 - 0.5
           n - The order of the filter, the higher n is the sharper
               the transition is. (n must be an integer >= 1).
Returns:
           f - Frequency domain filter of size==sze, the frequency
               origin is at the corners.

See also: lowpassfilter, highboostfilter, bandpassfilter

The logarithmic Gabor function in the frequency domain.

Usage: v = loggabor(f::Real, fo::Real, sigmaOnf::Real)

Arguments:
            f - Frequency to evaluate the function at.
           fo - Centre frequency of filter.
     sigmaOnf - Ratio of the standard deviation of the Gaussian
                describing the log Gabor filter's transfer function
                in the frequency domain to the filter center frequency.

sigmaOnf = 0.75 gives a filter bandwidth of about 1 octave.
sigmaOnf = 0.55 gives a filter bandwidth of about 2 octaves.

Generate arrays of filter grid angles.

Usage: (sintheta, costheta) = gridangles(freq, fx, fy)

Arguments: freq, fx, fy - The output of filtergrids()

Returns:       sintheta - The sine and cosine of the angles in the filtergrid
               costheta

See also filtergrids

Orientation selective frequency domain filter with cosine windowing function.

Usage: filter = cosineangularfilter(angl, wavelen, sintheta, costheta)

Arguments:
               angl - Orientation of the filter (radians)
            wavelen - Wavelength of the angular cosine window function.
 sintheta, costheta - Grids as returned by gridangles()

See also: gaussianangularfilter, filtergrids

Orientation selective frequency domain filter with Gaussian windowing function.

Usage: filter = gaussianangularfilter(angl, thetaSigma, sintheta, costheta)

Arguments:
               angl - Orientation of the filter (radians)
         thetasigma - Standard deviation of angular Gaussian window function.
 sintheta, costheta - Grids as returned by gridangles()

See also: cosineangularfilter, gridangles, filtergrids

2D Fourier transform of Moisan's periodic image component.

Usage: (P, S, p, s) = perfft2(img)

Argument: img - Image to be transformed
Returns:    P - 2D fft of periodic image component
            S - 2D fft of smooth component
            p - Periodic component (spatial domain)
            s - Smooth component (spatial domain)

Moisan's "Periodic plus Smooth Image Decomposition" decomposes an image into two components

    img = p + s

where s is the 'smooth' component with mean 0 and p is the 'periodic' component which has no sharp discontinuities when one moves cyclically across the image boundaries.

This decomposition is very useful when one wants to obtain an FFT of an image with minimal artifacts introduced from the boundary discontinuities. The image p gathers most of the image information but avoids periodization artifacts.

The typical use of this function is to obtain a 'periodic only' fft of an image

  P = perfft2(img)

Displaying the amplitude spectrum of P will yield a clean spectrum without the typical vertical-horizontal 'cross' arising from the image boundaries that you would normally see.

Note if you are using the function to perform filtering in the frequency domain you may want to retain s (the smooth component in the spatial domain) and add it back to the filtered result at the end.

The computational cost of obtaining the 'periodic only' FFT involves taking an additional FFT.

Generate geometric series.

Useful for generating geometrically scaled wavelengths for specifying filter banks.

Usage 1: s = geoseries(s1, mult, n)

Arguments:      s1 - The starting value in the series.
              mult - The scaling factor between succesive values.
                 n - The desired number of elements in the series.

Usage 2: s = geoseries((s1, sn), n)

Arguments: (s1, sn) - Tuple specifying the 1st and last values
                      in the the series.
                  n - The desired number of elements in the series.

Example:

      s = geoseries(0.5, 2, 4)
      s =  [0.5000,    1.0000,    2.0000,    4.0000]

Alternatively obtain the same series using

           s = geoseries((0.5, 4), 4)

A phase congruent test image that interpolates from a step to a line.

Generates a test image where the feature type changes from a step edge to a line feature from top to bottom. Gradient based edge detectors will only correctly mark the step-like feature towards the top of the image and incorrectly mark two features towards the bottom of the image whereas phase congruency will correctly mark a single feature from top to bottom. In general, natural images contain a roughly uniform distribution of the full continuum of feature types from step to line.

Usage:
    img = step2line(sze; nscales=50, ampexponent=-1, ncycles=1.5, phasecycles=0.25)

Arguments:
      sze::Integer - Number of rows in test image, defaults to 512.

Keyword Arguments:
  nscales::Integer - No of Fourier components used to construct the signal.
                     Defaults to 50.
 ampexponent::Real - Decay exponent of amplitude with frequency.
                     A value of -1 will produce amplitude inversely
                     proportional to frequency (corresponds to step feature).
                     A value of -2 will result in the line feature
                     appearing as a triangular waveform. Defaults to -1.
     ncycles::Real - Number of wave cycles across the width of the image.
                     Defaults to 1.5
 phasecycles::Real - Number of feature type phase cycles going vertically
                     down the image. Defaults to 0.25 giving a sequence of feature
                     phase congruency angle varying from 0 to pi/2.
Returns:
   img::Array{Float64,2} - The test image.


Examples of use:
  > img = step2line()                              # Default pattern
  > img = step2line(ncycles=3, ampexponent=-1.5);  # 3 cycles, 'soft' step to line
  > img = step2line(ncycles=3, ampexponent=-1.5, phasecycles = 3);

See also: circsine, starsine

Generate a phase congruent circular sine wave grating.

Useful for testing the isotropy of response of a feature dectector.

Usage:    img = circsine(sze; wavelength = 40, nscales = 50, ampexponent = -1,
                          offset = 0, p = 2, trim = false)

Arguments:
   sze::Integer  - The size of the square image to be produced. Defaults to 512.

Keyword arguments:
 wavelength::Real - The wavelength in pixels of the sine wave. Defaults to 40.
 nscales::Integer - No of Fourier components used to construct the
                    signal. This is typically 1, if you want a simple sine
                    wave, or >50 if you want to build a phase congruent
                    waveform. Defaults to 50.
ampexponent::Real - Decay exponent of amplitude with frequency.
                    A value of -1 will produce amplitude inversely
                    proportional to frequency (this will produce a step
                    feature if offset is 0)
                    A value of -2 with an offset of pi/2 will result in a
                    triangular waveform.  Defaults to -1;
     offset::Real - Angle of phase congruency at which the features of the
                    star pattern will be generated at. This controls the feature type.
                    0 for a step-like feature, pi/2 for a line/triangular-like feature.
                    Defaults to 0. If nscales = 1 use pi/2 to get continuity
                    at the centre.
      p::Integer  - Optional parameter specifying the norm to use in
                    calculating the radius from the centre. This defaults to
                    2, resulting in a circular pattern.  Large values gives
                    a square pattern
      trim::Bool  - Optional boolean flag indicating whether you want the
                    circular pattern trimmed from the corners leaving
                    only complete circles. Defaults to false.
Returns:
   img::Array{Float64,2} - The test image.

Examples:
> circsine(nscales = 1) - A simple circular sine wave pattern
> circsine(nscales = 50, ampexponent = -1, offset =  0)     - Square waveform
> circsine(nscales = 50, ampexponent = -2, offset = pi/2)   - Triangular waveform
> circsine(nscales = 50, ampexponent = -1.5, offset = pi/4) - Something in between
                                                              square and triangular
> circsine(nscales = 50, ampexponent = -1.5, offset = 0)    - Looks like a square but is not.

See also: starsine, step2line

Generate a phase congruent star shaped sine wave grating.

Useful for testing the behaviour of feature detectors at line junctions.

Usage:    img = starsine(sze; ncycles=10, nscales=50, ampexponent=-1, offset=0)

Argument:
     sze::Integer - The size of the square image to be produced. Defaults to 512.

Keyword arguments:
    ncycles::Real - The number of sine wave cycles around centre point.
                    Typically an integer, but any value can be used.
 nscales::Integer - No of fourier components used to construct the
                    signal. This is typically 1, if you want a simple sine
                    wave, or >50 if you want to build a phase congruent
                    waveform.  Defaults to 50.
ampexponent::Real - Decay exponent of amplitude with frequency.
                    A value of -1 will produce amplitude inversely
                    proportional to frequency (this will produce a step
                    feature if offset is 0)
                    A value of -2 with an offset of pi/2 will result in a
                    triangular waveform.
     offset::Real - Angle of phase congruency at which the features of the
                    star pattern will be generated at. This controls the feature type.
                    0 for a step-like feature, pi/2 for a line/triangular-like feature.
Returns:
   img::Array{Float64,2} - The test image.

Examples:
> starsine(nscales = 1) - A simple sine wave pattern radiating out
                          from the centre. Use 'offset' if you wish to
                          rotate it a bit.
> starsine(nscales = 50, ampexponent = -1, offset =  0)     - Square waveform
> starsine(nscales = 50, ampexponent = -2, offset = pi/2)   - Triangular waveform
> starsine(nscales = 50, ampexponent = -1.5, offset = pi/4) - Something in between
                                                              square and triangular
> starsine(nscales = 50, ampexponent = -1.5, offset = 0)    - Looks like a square but is not.

See also: circsine, step2line

Create $1/f^p$ spectrum noise images.

When displayed as a surface these images also generate great landscape terrain.

Usage: img = noiseonf(sze, p)

Arguments:
    sze::Tuple{Integer, Integer} or ::Integer
              - A tuple (rows, cols) or single value specifying size of
                image to produce.
      p::Real - Exponent of spectrum decay = 1/(f^p)

Returns:
   img::Array{Float64,2} - The noise image with specified spectrum.


Reference values for p:
            p = 0   - raw Gaussian noise image.
              = 1   - gives the supposedly 1/f 'standard' drop-off for
                      'natural' images.
              = 1.5 - seems to give the most interesting 'cloud patterns'.
              = > 2 - produces 'blobby' images.

Randomize image phase leaving amplitude spectrum unchanged.

Usage:   newimg = nophase(img)

Argument:       img::AbstractArray{T,2} where T <: Real - Input image

Returns:     newimg::Array{Float64,2} - Image with randomized phase

In general most images will be destroyed by this transform. However, some textures are reproduced in an 'amplitude only' image quite well. Typically these are textures which have an amplitude spectrum that have a limited number of isolated peaks. That is, a texture made up from a limited number of strong harmonics.

See also: noiseonf, quantizephase, swapphase

Quantize phase values in an image.

Usage:  qimg = quantizephase(img, N)

Arguments: img::Array{T,2} where T <: Real - Image to be processed
                     N::Integer - Desired number of quantized phase values

Returns:  qimg::Array{Float64,2} - Phase quantized image

Phase values in an image are important. However, despite this, they can be quantized very heavily with little perceptual loss. The value of N can be as low as 4, or even 3! Using N = 2 is also worth a look.

See also: swapphase

Demonstrates phase - amplitude swapping between images.

Usage:   (newimg1, newimg2) = swapphase(img1, img2)

Arguments:
 img1, img2::Array{<:Real,2} - Two images of same size to be used as input

Returns:
    newimg1::Array{Float64,2} - Image obtained from the phase of img1
                                and the magnitude of img2.
    newimg2::Array{Float64,2} - Phase of img2, magnitude of img1.

See also: quantizephase, nophase

Fill NaN values in an image with closest non NaN value.

This can be used as a crude (but quick) 'inpainting' function to allow a FFT to be computed on an image containing NaN values. While the 'inpainting' is crude it is typically good enough to remove most of the edge effects one might get at the boundaries of the NaN regions. The NaN regions should then be remasked out of the final processed image.

Usage:  (newim, mask) = fillnan(img)

  Argument:  img   - Image to be 'filled'.
  Returns:   newim - Filled image.
             mask  - Binary image indicating the valid, non-NaN, regions in
                     the original image.

See also: replacenan

Replace NaNs in an array with a specified value.

Usage: (newimg, mask) = replacenan(img, defaultval=0)

Arguments:
   img        - The Array containing NaN values.
   defaultval - The default value to replace NaNs.

Returns:
        newim - Filled image,
        mask  - Boolean image indicating non-NaN regions in the original
                image.

See also: fillnan

Hysteresis thresholding of an image.

Usage: bw = hysthresh(img, T1, T2)

Arguments:
           img  - Image to be thresholded
        T1, T2  - Upper and lower threshold values.  T1 and T2
                  can be entered in any order, the larger of the
                  two values is used as the upper threshold.
Returns:
            bw  - The binary thresholded image as a BitArray

All pixels with values above threshold T1 are marked as edges. All pixels that are connected to points that have been marked as edges and with values above threshold T2 are also marked as edges. Eight connectivity is used.