# ----------------------------------------------------------------------
# Numenta Platform for Intelligent Computing (NuPIC)
# Copyright (C) 2013, Numenta, Inc.  Unless you have an agreement
# with Numenta, Inc., for a separate license for this software code, the
# following terms and conditions apply:
#
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License version 3 as
# published by the Free Software Foundation.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
# See the GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program.  If not, see http://www.gnu.org/licenses.
#
# http://numenta.org/licenses/
# ----------------------------------------------------------------------

"""
## @file
This file defines GaborFilter, an ImageSensor filter that applies
one or more gabor filter masks to incoming images.
"""

import os
import shutil
import math

import numpy
from PIL import (Image,
                 ImageChops)
from nupic.regions.ImageSensorFilters.BaseFilter import BaseFilter
from nupic.math import GetNTAReal


# Separable convolution is not supported in the public release
try:
  from nupic.bindings.algorithms import (Float32SeparableConvolution2D,
                                         Float32Rotation45)
except:
  pass

dtype = GetNTAReal()


#+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+
# GaborFilter
#+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+
class GaborFilter(BaseFilter):

  """
  Apply a bank of Gabor filters to the original image, and
  return one or more images (same dimensions as the original)
  containing the Gabor responses.
  """

  #+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+
  # Public API methods
  #+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+

  def __init__(self, gaborBankParams=None,
                     debugMode=False,
                     debugOutputDir='gabor.d',
                     suppressLobes=False,
                     wipeOutsideMask=False,
                     convolutionMethod='1D'):
    """
    @param gaborBankParams -- A list of sub-lists that specify the parameters of
            the individual Gabor filters that comprise the filter bank.
            Each sub-list contains one or more dicts; each such dict specifies
            a particular Gabor filter.  The responses of each Gabor filter in
            a particular sub-list are combined together additively (using an
            L2 norm) to generate a single response image.

            Thus, the total number of response images generated by GaborFilter
            will equal the length of the main 'gaborBankParams' list.  The total
            number of individual Gabor filters used to generate these responses
            will equal to the total number of dicts contained in all of the
            sub-lists of 'gaborBankParams'.

            Each of the dicts is expected to contain the following
            required keys:

            'scale': an integer between 3 and 18 that specifies the 'radius'
                    of the gabor filter.  For example, if 'scale' is set to 5,
                    then the resulting filter will have size 11x11 pixels.

            'orient': the orientation of the filter, in degrees, with zero
                    corresponding to vertical, and increasing in a clockwise
                    manner (e.g., a value of 90 produces a horizontally-oriented
                    filter; a value of 45 degrees produces a diagonally-oriented
                    filter along the southwest-to-northeast axis, etc.)
                    Legal values are: 0, 45, 90, and 135.

            'phase': the phase of the filter, in degrees.  Legal values are: 0,
                    90, 180, and 270.  Phases of 0 and 180 produce single-lobed
                    symmetric gabor filters that are generally selective to lines.
                    Phases of 90 and 270 produce dual-lobed (asymmetric)
                    gabor filters that are generally selective to edges.

            'mode': specifies the post-processing mode to apply to the raw
                    filter output.  Legal values are:

                    'raw'       Perform no post-processing. Outputs will be in (-1.0, +1.0)
                    'abs'       Output is absolute value of raw response.  Output will
                                lie in range (0.0, +1.0)
                    'positive'  Clip negative raw response values to zero.  Output will lie
                                in range (0.0, +1.0)
                    'rectified' Output is square of input.  Output will lie in the
                                range (0.0, +1.0)
                    'power'     Positive raw responses are squared; negative raw responses
                                are clipped to zero.  Output will be in (0.0, +1.0)
                    'hardFull'  Apply a hard threshold to the input; if raw response is
                                >= 'threshold', then output is 1.0, else output is -1.0.
                                Output lies in (-1.0, +1.0)
                    'hardHalf'  Apply a hard threshold to the input; if raw response is
                                >= 'threshold', then output is 1.0, else output is 0.0.
                                Output lies in (0.0, +1.0)
                    'sigmoidFull': Apply a sigmoidal threshold function to the input using
                                'threshold' and 'steepness'.  Output lies in (-1.0, +1.0).
                    'sigmoidHalf': Apply a sigmoidal threshold function to the input using
                                'threshold' and 'steepness'.  Output lies in (0.0, +1.0).

            'name': a human-meaningful name (primarily used for debugging purposes.)

            In addition, the following key is required if 'mode' is set to 'hardFull',
            'hardHalf', 'sigmoidFull', or 'sigmoidHalf'

            'threshold': the threshold to use for either the hard ('hardFull' or 'hardHalf')
            or soft ('sigmoidFull' or 'sigmoidHalf') thresholding post-processing modes.

            In addition, the following key is required if 'mode' is set to 'sigmoidFull'
            or 'sigmoidHalf':

            'steepness': controls the steepness of the sigmoidal function that performs
            post-processing on the raw response.

            If 'gaborBankParams' is None, then a default set of filter
            bank parameters will be used.  This default set will consist of 36
            separate responses, each computed from a single underlying Gabor filter.

    @param debugMode -- A boolean flag indicating whether or not the filter
            should output debugging information to a file (default is False.)

    @param debugOutputDir -- The name of the output directory that will be
            created to store gabor response images in the event that the
            debugging flag 'debugMode' is set to True.  Defaults to 'gabor.d'

    @param suppressLobes -- A boolean flag indicating whether or not the secondary
            lobes of the gabor filters should be zeroed out and thus removed
            from the mask.  Defaults to False.

    @param convolutionMethod -- Method to use for convolving filters with images.
            '2D' is straightforward convolution of 2D filter with image.  Other options
            are 'FFT' for convolution by multiplication in Fourier space, and '1D' for
            convolution with two 1D filters formed from a separable 2D filter. '1D'
            convolution is not yet available in our public release.
    """

    BaseFilter.__init__(self)

    # Debugging
    self._debugMode = debugMode
    self._debugOutputDir = debugOutputDir
    self._suppressLobes = suppressLobes
    self._wipeOutsideMask = wipeOutsideMask
    self._convolutionMethod = convolutionMethod

    # Separable convolution is not supported in the public release
    if self._convolutionMethod == '1D':
      try:
        Float32SeparableConvolution2D()
        Float32Rotation45()
      except:
        self._convolutionMethod = 'FFT'

    if self._convolutionMethod == '1D':
      # Create an instance of the rotation class
      self._rotation = Float32Rotation45()

    self._debugCompositeDir = os.path.join(self._debugOutputDir, 'composite')

    # Store the parameter set
    if gaborBankParams is None:
      gaborBankParams = self.getDefaultParamSet()
    self._gaborBankParams = gaborBankParams

    self._defAspectRatio = 1.30

    if self._debugMode:
      print 'Using Gabor bank parameters:'
      print self._gaborBankParams

    # Prepare the filters
    self._gaborBank = self._makeGaborBank()

    if self._debugMode:
      print 'Gabor Bank:'
      for f in self._gaborBank:
        print '+=+=+=+=+=+=+=+=+=+=+=+=+=+=+'
        print 'Scale:  %d' % f[0][1]['scale']
        print 'Orient: %d' % f[0][1]['orient']
        print 'Phase:  %d' % int(f[0][1]['phase'])
        self._printGaborFilter(f[0][0])

    if self._debugMode:
      print 'Gabor bank generated with %d filter pairs.' % len(self._gaborBankParams)

    # Initialize
    self._imageCounter = 0

  def getOutputCount(self):
    """
    Return the number of images returned by each call to process().

    If the filter creates multiple simultaneous outputs, return a tuple:
    (outputCount, simultaneousOutputCount).
    """

    return (1, len(self._gaborBank))

  def process(self, origImage):
    """
    Perform Gabor filtering on the input image and return one or more
    response maps having dimension identical to the input image.

    @param origImage -- The image to process.

    Returns a list containing a list of the response maps.
    """

    if self.mode != 'gray':
      raise RuntimeError("GaborFilter only supports grayscale images.")

    BaseFilter.process(self, origImage)
    responses = []

    if self._debugMode:
      print 'GaborFilter: process()'
      print 'GaborFilter: origImage size:', origImage.size

    # Convert image data from PIL to numpy array
    imageData = numpy.asarray(origImage.split()[0], dtype=dtype)
    # Save the mask to put back later
    mask = origImage.split()[1]
    if self._wipeOutsideMask:
      maskData = numpy.asarray(mask, dtype=dtype)
    else:
      maskData = None

    # Perform the actual Gabor filtering
    responseSet = self._doProcessing(imageData, maskData)

    if self._debugMode:
      print 'Responses generated: %d' % len(responseSet)

     # Convert from numpy to PIL
    imageSet = self._convertToPIL(responseSet)

    # Dump the gabor responses to disk
    if self._debugMode:
      self._saveDebugImages(imageSet)

    # Add the mask back
    for image in imageSet:
      image.putalpha(mask)

    self._imageCounter += 1
    return [imageSet]

  def processArray(self, origImageData):
    """
    Perform Gabor filtering on the input array and return one or more
    response maps having dimension identical to the input image array.

    @param origImageData -- two-dimensional numpy array representing the image.

    Returns a list containing the response maps. Each map is a 2D numpy array.
    """


    if self._debugMode:
      print 'GaborFilter: processArray()'

    # Perform the actual Gabor filtering
    responseSet = self._doProcessing(origImageData)

    if self._debugMode:
      print 'Responses generated: %d' % len(responseSet)

    self._imageCounter += 1
    return responseSet


  def getDefaultParamSet(self):
    """
    Provide a default set of Gabor filter bank parameters in the event
    that none were provided by the application.

    Note: this method is classified as a public method so as to allow
    the application to query the default parameter set.
    """

    postProcMode = 'sigmoidHalf'
    scale = {
             'small':   3,
             'medium':  6,
             'large':  12,
            }
    orient = {
              'vert':   0,
              'swne':  45,
              'horz':  90,
              'nwse': 135,
             }
    phase = {
              'line-bright':           0,
              'edge-bright-to-dark':  90,
              'line-dark':           180,
              'edge-dark-to-bright': 270,
             }

    return [

        #+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+
        # Small vertical lines and edges
        #+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+

        # Small vertical lines
        [
          # Phase-0 (detects bright line against dark background)
          {
            'scale': scale['small'],
            'orient': orient['vert'],
            'phase': phase['line-bright'],
            'mode': postProcMode,
            'threshold': 0.400, #0.300,
            'steepness': 8,
            'name': 'small-vert-line-bright',
           },
        ],
        # [
        #   # Phase-180 (detects dark line against bright background)
        #   {
        #    'scale': scale['small'],
        #    'orient': orient['vert'],
        #    'phase': phase['line-dark'],
        #    'mode': postProcMode,
        #    'threshold': 0.400, #0.300,
        #    'steepness': 8,
        #    'name': 'small-vert-line-dark',
        #   },
        # ],

        # Small vertical edges (bright to dark)
        [
          # Phase-90 (detects bright to dark edges)
          {
            'scale': scale['small'],
            'orient': orient['vert'],
            'phase': phase['edge-bright-to-dark'],
            'mode': postProcMode,
            'threshold': 0.400, #0.300,
            'steepness': 8,
            'name': 'small-vert-edge-bright',
           },
         ],

        # Small vertical edges (dark to bright)
        [
          # Phase-270 (detects dark to bright edges)
          {
            'scale': scale['small'],
            'orient': orient['vert'],
            'phase': phase['edge-dark-to-bright'],
            'mode': postProcMode,
            'threshold': 0.400, #0.300,
            'steepness': 8,
            'name': 'small-vert-edge-dark',
           },
         ],

        #+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+
        # Small SW-NE lines and edges
        #+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+

        # Small SW-NE lines
        [
          # Phase-0 (detects bright line against dark background)
          {
            'scale': scale['small'],
            'orient': orient['swne'],
            'phase': phase['line-bright'],
            'mode': postProcMode,
            'threshold': 0.400,
            'steepness': 8,
            'name': 'small-swne-line-bright',
           },
        ],
        # [
        #   # Phase-180 (detects dark line against bright background)
        #   {
        #    'scale': scale['small'],
        #    'orient': orient['swne'],
        #    'phase': phase['line-dark'],
        #    'mode': postProcMode,
        #    'threshold': 0.400,
        #    'steepness': 8,
        #    'name': 'small-swne-line-dark',
        #   },
        # ],

        # Small SW-NE edges (bright to dark)
        [
          # Phase-90 (detects bright to dark edges)
          {
            'scale': scale['small'],
            'orient': orient['swne'],
            'phase': phase['edge-bright-to-dark'],
            'mode': postProcMode,
            'threshold': 0.400,
            'steepness': 8,
            'name': 'small-swne-edge-bright',
           },
         ],

        # Small SW-NE edges (dark to bright)
        [
          # Phase-270 (detects dark to bright edges)
          {
            'scale': scale['small'],
            'orient': orient['swne'],
            'phase': phase['edge-dark-to-bright'],
            'mode': postProcMode,
            'threshold': 0.400,
            'steepness': 8,
            'name': 'small-swne-edge-dark',
           },
         ],

        #+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+
        # Small horizontal lines and edges
        #+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+

        # Small horizontal lines
        [
          # Phase-0 (detects bright line against dark background)
          {
            'scale': scale['small'],
            'orient': orient['horz'],
            'phase': phase['line-bright'],
            'mode': postProcMode,
            'threshold': 0.400, #0.300,
            'steepness': 8,
            'name': 'small-horz-line-bright',
           },
        ],
        # [
        #   # Phase-180 (detects dark line against bright background)
        #   {
        #    'scale': scale['small'],
        #    'orient': orient['horz'],
        #    'phase': phase['line-dark'],
        #    'mode': postProcMode,
        #    'threshold': 0.400,
        #    'steepness': 8,
        #    'name': 'small-horz-line-dark',
        #   },
        # ],

        # Small horizontal edges (bright to dark)
        [
          # Phase-90 (detects bright to dark edges)
          {
            'scale': scale['small'],
            'orient': orient['horz'],
            'phase': phase['edge-bright-to-dark'],
            'mode': postProcMode,
            'threshold': 0.400, #0.300,
            'steepness': 8,
            'name': 'small-horz-edge-bright',
           },
         ],

        # Small horizontal edges (dark to bright)
        [
          # Phase-270 (detects dark to bright edges)
          {
            'scale': scale['small'],
            'orient': orient['horz'],
            'phase': phase['edge-dark-to-bright'],
            'mode': postProcMode,
            'threshold': 0.400, #0.300,
            'steepness': 8,
            'name': 'small-horz-edge-dark',
           },
         ],

        #+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+
        # Small SE-NW lines and edges
        #+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+

        # Small SE-NW lines
        [
          # Phase-0 (detects bright line against dark background)
          {
            'scale': scale['small'],
            'orient': orient['nwse'],
            'phase': phase['line-bright'],
            'mode': postProcMode,
            'threshold': 0.400,
            'steepness': 8,
            'name': 'small-nwse-line-bright',
           },
        ],
        # [
        #   # Phase-180 (detects dark line against bright background)
        #   {
        #    'scale': scale['small'],
        #    'orient': orient['nwse'],
        #    'phase': phase['line-dark'],
        #    'mode': postProcMode,
        #    'threshold': 0.400,
        #    'steepness': 8,
        #    'name': 'small-nwse-line-dark',
        #   },
        # ],

        # Small SE-NW edges (bright to dark)
        [
          # Phase-90 (detects bright to dark edges)
          {
            'scale': scale['small'],
            'orient': orient['nwse'],
            'phase': phase['edge-bright-to-dark'],
            'mode': postProcMode,
            'threshold': 0.400,
            'steepness': 8,
            'name': 'small-nwse-edge-bright',
           },
         ],

        # Small SE-NW edges (dark to bright)
        [
          # Phase-270 (detects dark to bright edges)
          {
            'scale': scale['small'],
            'orient': orient['nwse'],
            'phase': phase['edge-dark-to-bright'],
            'mode': postProcMode,
            'threshold': 0.400,
            'steepness': 8,
            'name': 'small-nwse-edge-dark',
           },
         ],


        #+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+
        # Medium vertical lines and edges
        #+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+

        # Medium vertical lines
        [
          # Phase-0 (detects bright line against dark background)
          {
            'scale': scale['medium'],
            'orient': orient['vert'],
            'phase': phase['line-bright'],
            'mode': postProcMode,
            'threshold': 0.500, #0.450,
            'steepness': 6,
            'name': 'medium-vert-line-bright',
           },
        ],
        # [
        #   # Phase-180 (detects dark line against bright background)
        #   {
        #    'scale': scale['medium'],
        #    'orient': orient['vert'],
        #    'phase': phase['line-dark'],
        #    'mode': postProcMode,
        #    'threshold': 0.500,
        #    'steepness': 6,
        #    'name': 'medium-vert-line-dark',
        #   },
        # ],

        # Medium vertical edges (bright to dark)
        [
          # Phase-90 (detects bright to dark edges)
          {
            'scale': scale['medium'],
            'orient': orient['vert'],
            'phase': phase['edge-bright-to-dark'],
            'mode': postProcMode,
            'threshold': 0.500, #0.450,
            'steepness': 6,
            'name': 'medium-vert-edge-bright',
           },
         ],

        # Medium vertical edges (dark to bright)
        [
          # Phase-270 (detects dark to bright edges)
          {
            'scale': scale['medium'],
            'orient': orient['vert'],
            'phase': phase['edge-dark-to-bright'],
            'mode': postProcMode,
            'threshold': 0.500, #0.450,
            'steepness': 6,
            'name': 'medium-vert-edge-dark',
           },
         ],

        #+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+
        # Medium SW-NE lines and edges
        #+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+

        # Medium SW-NE lines
        [
          # Phase-0 (detects bright line against dark background)
          {
            'scale': scale['medium'],
            'orient': orient['swne'],
            'phase': phase['line-bright'],
            'mode': postProcMode,
            'threshold': 0.600,
            'steepness': 6,
            'name': 'medium-swne-line-bright',
           },
        ],
        # [
        #   # Phase-180 (detects dark line against bright background)
        #   {
        #    'scale': scale['medium'],
        #    'orient': orient['swne'],
        #    'phase': phase['line-dark'],
        #    'mode': postProcMode,
        #    'threshold': 0.600,
        #    'steepness': 6,
        #    'name': 'medium-swne-line-dark',
        #   },
        # ],

        # Medium SW-NE edges (bright to dark)
        [
          # Phase-90 (detects bright to dark edges)
          {
            'scale': scale['medium'],
            'orient': orient['swne'],
            'phase': phase['edge-bright-to-dark'],
            'mode': postProcMode,
            'threshold': 0.600,
            'steepness': 6,
            'name': 'medium-swne-edge-bright',
           },
         ],

        # Medium SW-NE edges (dark to bright)
        [
          # Phase-270 (detects dark to bright edges)
          {
            'scale': scale['medium'],
            'orient': orient['swne'],
            'phase': phase['edge-dark-to-bright'],
            'mode': postProcMode,
            'threshold': 0.600,
            'steepness': 6,
            'name': 'medium-swne-edge-dark',
           },
         ],

        #+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+
        # Medium horizontal lines and edges
        #+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+
        [
          # Phase-0 (detects bright line against dark background)
          {
            'scale': scale['medium'],
            'orient': orient['horz'],
            'phase': phase['line-bright'],
            'mode': postProcMode,
            'threshold': 0.500, #0.450,
            'steepness': 6,
            'name': 'medium-horz-line-bright',
           },
        ],
        # [
        #   # Phase-180 (detects dark line against bright background)
        #   {
        #    'scale': scale['medium'],
        #    'orient': orient['horz'],
        #    'phase': phase['line-dark'],
        #    'mode': postProcMode,
        #    'threshold': 0.500,
        #    'steepness': 6,
        #    'name': 'medium-horz-line-dark',
        #   },
        # ],

        # Medium horizontal edges (bright to dark)
        [
          # Phase-90 (detects bright to dark edges)
          {
            'scale': scale['medium'],
            'orient': orient['horz'],
            'phase': phase['edge-bright-to-dark'],
            'mode': postProcMode,
            'threshold': 0.500, #0.450,
            'steepness': 6,
            'name': 'medium-horz-edge-bright',
           },
         ],

        # Medium horizontal edges (dark to bright)
        [
          # Phase-270 (detects dark to bright edges)
          {
            'scale': scale['medium'],
            'orient': orient['horz'],
            'phase': phase['edge-dark-to-bright'],
            'mode': postProcMode,
            'threshold': 0.500, #0.450,
            'steepness': 6,
            'name': 'medium-horz-edge-dark',
           },
         ],

        #+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+
        # Medium SE-NW lines and edges
        #+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+

        # Medium SE-NW lines
        [
          # Phase-0 (detects bright line against dark background)
          {
            'scale': scale['medium'],
            'orient': orient['nwse'],
            'phase': phase['line-bright'],
            'mode': postProcMode,
            'threshold': 0.600,
            'steepness': 6,
            'name': 'medium-nwse-line-bright',
           },
        ],
        # [
        #   # Phase-180 (detects dark line against bright background)
        #   {
        #    'scale': scale['medium'],
        #    'orient': orient['nwse'],
        #    'phase': phase['line-dark'],
        #    'mode': postProcMode,
        #    'threshold': 0.600,
        #    'steepness': 6,
        #    'name': 'medium-nwse-line-dark',
        #   },
        # ],

        # Medium SE-NW edges (bright to dark)
        [
          # Phase-90 (detects bright to dark edges)
          {
            'scale': scale['medium'],
            'orient': orient['nwse'],
            'phase': phase['edge-bright-to-dark'],
            'mode': postProcMode,
            'threshold': 0.600,
            'steepness': 6,
            'name': 'medium-nwse-edge-bright',
           },
         ],

        # Medium SE-NW edges (dark to bright)
        [
          # Phase-270 (detects dark to bright edges)
          {
            'scale': scale['medium'],
            'orient': orient['nwse'],
            'phase': phase['edge-dark-to-bright'],
            'mode': postProcMode,
            'threshold': 0.600,
            'steepness': 6,
            'name': 'medium-nwse-edge-dark',
           },
         ],


        #+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+
        # Large vertical lines and edges
        #+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+

        # Large vertical lines
        [
          # Phase-0 (detects bright line against dark background)
          {
            'scale': scale['large'],
            'orient': orient['vert'],
            'phase': phase['line-bright'],
            'mode': postProcMode,
            'threshold': 0.600,
            'steepness': 6,
            'name': 'large-vert-line-bright',
           },
        ],
        # [
        #   # Phase-180 (detects dark line against bright background)
        #   {
        #    'scale': scale['large'],
        #    'orient': orient['vert'],
        #    'phase': phase['line-dark'],
        #    'mode': postProcMode,
        #    'threshold': 0.600,
        #    'steepness': 6,
        #    'name': 'large-vert-line-dark',
        #   },
        # ],

        # Large vertical edges (bright to dark)
        [
          # Phase-90 (detects bright to dark edges)
          {
            'scale': scale['large'],
            'orient': orient['vert'],
            'phase': phase['edge-bright-to-dark'],
            'mode': postProcMode,
            'threshold': 0.600,
            'steepness': 6,
            'name': 'large-vert-edge-bright',
           },
         ],

        # Large vertical edges (dark to bright)
        [
          # Phase-270 (detects dark to bright edges)
          {
            'scale': scale['large'],
            'orient': orient['vert'],
            'phase': phase['edge-dark-to-bright'],
            'mode': postProcMode,
            'threshold': 0.600,
            'steepness': 6,
            'name': 'large-vert-edge-dark',
           },
         ],

        #+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+
        # Large SW-NE lines and edges
        #+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+

        # Large SW-NE lines
        [
          # Phase-0 (detects bright line against dark background)
          {
            'scale': scale['large'],
            'orient': orient['swne'],
            'phase': phase['line-bright'],
            'mode': postProcMode,
            'threshold': 0.600,
            'steepness': 6,
            'name': 'large-swne-line-bright',
           },
        ],
        # [
        #   # Phase-180 (detects dark line against bright background)
        #   {
        #    'scale': scale['large'],
        #    'orient': orient['swne'],
        #    'phase': phase['line-dark'],
        #    'mode': postProcMode,
        #    'threshold': 0.600,
        #    'steepness': 6,
        #    'name': 'large-swne-line-dark',
        #   },
        # ],

        # Large SW-NE edges (bright to dark)
        [
          # Phase-90 (detects bright to dark edges)
          {
            'scale': scale['large'],
            'orient': orient['swne'],
            'phase': phase['edge-bright-to-dark'],
            'mode': postProcMode,
            'threshold': 0.600,
            'steepness': 6,
            'name': 'large-swne-edge-bright',
           },
         ],

        # Large SW-NE edges (dark to bright)
        [
          # Phase-270 (detects dark to bright edges)
          {
            'scale': scale['large'],
            'orient': orient['swne'],
            'phase': phase['edge-dark-to-bright'],
            'mode': postProcMode,
            'threshold': 0.600,
            'steepness': 6,
            'name': 'large-swne-edge-dark',
           },
         ],

        #+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+
        # Large horizontal lines and edges
        #+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+

        # Large horizontal lines
        [
          # Phase-0 (detects bright line against dark background)
          {
            'scale': scale['large'],
            'orient': orient['horz'],
            'phase': phase['line-bright'],
            'mode': postProcMode,
            'threshold': 0.600,
            'steepness': 6,
            'name': 'large-horz-line-bright',
           },
        ],
        # [
        #   # Phase-180 (detects dark line against bright background)
        #   {
        #    'scale': scale['large'],
        #    'orient': orient['horz'],
        #    'phase': phase['line-dark'],
        #    'mode': postProcMode,
        #    'threshold': 0.600,
        #    'steepness': 6,
        #    'name': 'large-horz-line-dark',
        #   },
        # ],

        # Large horizontal edges (bright to dark)
        [
          # Phase-90 (detects bright to dark edges)
          {
            'scale': scale['large'],
            'orient': orient['horz'],
            'phase': phase['edge-bright-to-dark'],
            'mode': postProcMode,
            'threshold': 0.600,
            'steepness': 6,
            'name': 'large-horz-edge-bright',
           },
         ],

        # Large horizontal edges (dark to bright)
        [
          # Phase-270 (detects dark to bright edges)
          {
            'scale': scale['large'],
            'orient': orient['horz'],
            'phase': phase['edge-dark-to-bright'],
            'mode': postProcMode,
            'threshold': 0.600,
            'steepness': 6,
            'name': 'large-horz-edge-dark',
           },
         ],

        #+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+
        # Large SE-NW lines and edges
        #+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+

        # Large SE-NW lines
        [
          # Phase-0 (detects bright line against dark background)
          {
            'scale': scale['large'],
            'orient': orient['nwse'],
            'phase': phase['line-bright'],
            'mode': postProcMode,
            'threshold': 0.600,
            'steepness': 6,
            'name': 'large-nwse-line-bright',
           },
        ],
        # [
        #   # Phase-180 (detects dark line against bright background)
        #   {
        #    'scale': scale['large'],
        #    'orient': orient['nwse'],
        #    'phase': phase['line-dark'],
        #    'mode': postProcMode,
        #    'threshold': 0.600,
        #    'steepness': 6,
        #    'name': 'large-nwse-line-dark',
        #   },
        # ],

        # Large SW-NE edges (bright to dark)
        [
          # Phase-90 (detects bright to dark edges)
          {
            'scale': scale['large'],
            'orient': orient['nwse'],
            'phase': phase['edge-bright-to-dark'],
            'mode': postProcMode,
            'threshold': 0.600,
            'steepness': 6,
            'name': 'large-nwse-edge-bright',
           },
         ],

        # Large SW-NE edges (dark to bright)
        [
          # Phase-270 (detects dark to bright edges)
          {
            'scale': scale['large'],
            'orient': orient['nwse'],
            'phase': phase['edge-dark-to-bright'],
            'mode': postProcMode,
            'threshold': 0.600,
            'steepness': 6,
            'name': 'large-nwse-edge-dark',
           },
         ],
    ]

  def createDebugImages(self, imageSet):
    """
    Organize and return Gabor response images and composite images, given the
    filter responses in imageSet.

    Returns a list of dictionaries, one per filter, where each
    dictionary item contains:
       'scale'    - scale
       'orient'   - orientation
       'image'    - raw response image for this scale and orientation
       'name'     - filter name
       'counter'  - internal image counter, incremented for each new source image

    Returns a list of dictionaries, one per scale, where each
    dictionary item contains:
        'scale'   - scale
        'image'   - composite image for this scale, over all orientations
        'counter' - internal image counter, incremented for each new source image
    """

    # Collect the raw responses
    rawResponses = []
    scaleDict = {}
    for responseIndex in range(len(imageSet)):
      r = {}
      r['image'] = imageSet[responseIndex].split()[0]
      r['name'] = self._gaborBankParams[responseIndex][0]['name']
      r['scale'] = self._gaborBankParams[responseIndex][0]['scale']
      r['orient'] = self._gaborBankParams[responseIndex][0]['orient']
      r['counter'] = self._imageCounter
      rawResponses.append(r)

      # Accounting for composite image creation
      if self._gaborBankParams[responseIndex][0]['scale'] not in scaleDict:
        scaleDict[self._gaborBankParams[responseIndex][0]['scale']] = []
      scaleDict[self._gaborBankParams[responseIndex][0]['scale']] += [(responseIndex,
                      self._gaborBankParams[responseIndex][0]['orient'])]

    # Create the composite responses
    compositeResponses = []
    for scale in scaleDict.keys():
      # Accumulate image indices for each orientation in this scale
      composites = {'red': [], 'green': [], 'blue': []}
      for responseIndex, orient in scaleDict[scale]:
        if orient == 90:
          composites['red'] += [responseIndex]
        elif orient == 0:
          composites['green'] += [responseIndex]
        else:
          composites['blue'] += [responseIndex]

      # Generate RGB composite images for this scale
      images = {'red': None, 'green': None, 'blue': None}
      bands = []
      for color in composites.keys():
        imgList = composites[color]
        im = Image.new('L',imageSet[0].size)
        for indx in imgList:
          im2 = rawResponses[indx]['image']
          im = ImageChops.add(im, im2)
        images[color] = im
        bands += [images[color]]

      # Make final composite for this scale
      m = Image.merge(mode='RGB', bands=bands)
      r = {}
      r['scale'] = scale
      r['counter'] = self._imageCounter
      r['image'] = m
      compositeResponses.append(r)

    return rawResponses, compositeResponses

  #+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+
  # Private (internal) methods
  #+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+

  def _createDir(self, directory):
    """
    Create a directory for writing debug images
    """
    print 'Wiping directory tree: [%s]...' % directory
    if os.path.exists(directory):
      shutil.rmtree(directory)
    print 'Creating directory: [%s]...' % directory
    os.makedirs(directory)

  def _saveDebugImages(self, imageSet, amplification=1.0):
    """
    Save the Gabor responses to disk.
    """
    # Create composites and organize response images
    rawResponses, compositeResponses = self.createDebugImages(imageSet)

    # Create output directory if needed
    if self._imageCounter == 0:
      self._createDir(self._debugOutputDir)
      self._createDir(self._debugCompositeDir)

    # Save the responses
    for r in rawResponses:
      dstPath = os.path.join(self._debugOutputDir,
                             'image-%04d.%s.png' % (r['counter'], r['name']))
      r['image'].save(dstPath)

    # Save each composite image
    for c in compositeResponses:
      dstPath = os.path.join(self._debugCompositeDir,
                             'image-%04d.composite-%02d.png' % (c['counter'], c['scale']))
      fp = open(dstPath, 'w')
      c['image'].save(fp)
      fp.close()


  def _makeGaborBank(self):
    """
    Generate a complete bank of Gabor filters to the specification
    contained in self._gaborBankParams.

    Return value: a list of sub-lists, each sub-list containing a
    2-tuple of the form (filter, filterParams), where filter is
    a numpy array instantiating the Gabor filter, and 'filterParams'
    is a reference to the relevant filter description from
    'self._gaborBankParams'.
    """
    gaborBank = []
    for responseSet in self._gaborBankParams:
      filterSet = []
      for filterParams in responseSet:
        filter = self._makeGaborFilter(filterParams)
        self._normalizeGain(filter, filterParams)
        filterSet += [(filter, filterParams)]
      gaborBank += [filterSet]
    if self._debugMode:
      print '_makeGaborBank: %d filters generated' % len(gaborBank)

    return gaborBank


  def _normalizeGain(self, filter, filterParams):
    """
    This method normalizes the gain of the filter by adding a DC offset. On the
    two orientations used by the convolve1D function, it noramlizes rows and
    columns individually to ensure that the filter remains separable.

    Previously, this method always shifted the mean of the entire 2D filter,
    which produced a filter that was no longer separable. This method performs
    differently on the 0 and 90 degree orientations. It normalizes rows when
    rows hold the sinusoidal component of the Gabor filter (orientation 0) and
    columns otherwise (orientation 90). It falls back to normalizing the entire
    filter together for the other two orientations, because they are not used
    by the convolve1D function and are not separable.

    It is only necessary on even-symmetric filters (phase 0 and 180), as the
    sinusoidal components of the odd-symmetric filters always have 0 mean.

    The responses created using this method and those created using the old
    version of _normalizeGain differ slightly when orientation is 0 or 90 and
    phase is 0 or 180, but not significantly.
    """

    # if filterParams['orient'] in (0, 90):
    if filterParams['phase'] in (0, 180):
      # Only necessary on the even-symmetric filters
      if filterParams['orient'] in (0, 45):
        for row in xrange(filter.shape[0]):
          filter[row,:] -= filter[row,:].mean()
      else:
        for col in xrange(filter.shape[1]):
          filter[:,col] -= filter[:,col].mean()
    # else:
    #   # Old normalizeGain method
    #   filter -= filter.mean()


  def _applyDefaults(self, filterParams):
    """
    Instantiate filter parameters 'aspect', 'width', and 'wavelength'
    (if they are not already explicitly specified in the filter
    description) by deriving values based upon existing values.
    """

    if 'aspect' not in filterParams:
      filterParams['aspect'] = self._defAspectRatio

    defWidth, defWavelength = self._chooseParams(filterParams['scale'])
    if 'width' not in filterParams:
      filterParams['width'] = defWidth
    if 'wavelength' not in filterParams:
      filterParams['wavelength'] = defWavelength


  def _makeGaborFilter(self, filterParams):
    """
    Generate a single Gabor filter.
    @param filterParams -- a dict containing a description of the
              desired filter.
    """
    self._applyDefaults(filterParams)

    if self._debugMode:
      print '_makeGaborFilter():'
      print '  maskRadius:', filterParams['scale']
      print '  effectiveWidth:', filterParams['width']
      print '  aspectRatio:', filterParams['aspect']
      print '  orientation:', filterParams['orient']
      print '  wavelength:', filterParams['wavelength']
      print '  phase:', filterParams['phase']

    # Deg-to-rad
    orientation = filterParams['orient']
    if self._convolutionMethod == '1D' and orientation in (45, 135):
      # Rotate the filter 45 degrees counterclockwise for separable convolution
      # Filter will be vertical or horizontal, and image will be rotated
      orientation -= 45
    orientation = self._deg2rad(orientation)
    phase = self._deg2rad(filterParams['phase'])

    # Create the mask lattice
    maskDim = filterParams['scale'] * 2 + 1
    # if filterParams['phase'] in [0, 180]:
    #     maskDim += 1

    g = numpy.zeros((maskDim, maskDim), dtype=dtype)
    x = numpy.zeros((maskDim, maskDim), dtype=dtype)
    y = numpy.zeros((maskDim, maskDim), dtype=dtype)

    # Create the latttice points
    halfWidth = 0.5 * float(maskDim - 1)
    for j in range(maskDim):
        for i in range(maskDim):
            x[j,i] = float(i-halfWidth)
            y[j,i] = float(j-halfWidth)

    # Generate gabor mask
    cosTheta = math.cos(orientation)
    sinTheta = math.sin(orientation)
    sinusoidalConstant = 2.0 * math.pi / filterParams['wavelength']
    exponentialConstant = -0.5 / (filterParams['width']* filterParams['width'])
    aspectConstant = filterParams['aspect'] * filterParams['aspect']
    for j in range(maskDim):
        for i in range(maskDim):
            x0 = x[j,i] * cosTheta + y[j,i] * sinTheta
            y0 = y[j,i] * cosTheta - x[j,i] * sinTheta
            sinusoidalTerm = math.cos((sinusoidalConstant * x0) + phase)
            exponentialTerm = math.exp(exponentialConstant * (x0*x0 + aspectConstant*y0*y0))
            g[j,i] = exponentialTerm * sinusoidalTerm

    # Suppress lobes (optional)
    if self._suppressLobes:
      g = self._doLobeSuppression(g, orientation, phase)

    return g


  def _deg2rad(self, degrees):
    """
    Utility macro for converting from degrees to radians.
    """
    return degrees * math.pi / 180.0


  def _rad2deg(self, radians):
    """
    Utility macro for converting from radians to degrees.
    """
    return radians * 180.0 / math.pi


  def _chooseParams(self, scale):
    """
    Returns a 2-tuple of (width, wavelength) containing reasonable
    default values for a particular 'scale'.
    """
    paramTable = [
        (1,   0.8,  1.3),
        (2,   1.7,  2.4),
        (3,   2.8,  3.5),
        (4,   3.6,  4.6),
        (5,   4.5,  5.6),
        (6,   5.4,  6.8),
        (7,   6.3,  7.9),
        (8,   7.3,  9.1),
        (9,   8.2, 10.3),
        (10,  9.2, 11.5),
        (11, 10.2, 12.7),
        (12, 11.3, 14.1),
        (13, 12.3, 15.4),
        (14, 13.4, 16.8),
        (15, 14.6, 18.2),
        (16, 15.8, 19.7),
        (17, 17.0, 21.2),
        (18, 18.2, 22.8),
            ]
    for paramSet in paramTable:
      if paramSet[0] == scale:
        return paramSet[1], paramSet[2]


  def _suppressionRules(self):
    """
    Return a table of algorithmic parameters used in lobe suppression.
    The table contains a list of 2-tuples of the form:
    (orientation, phase), (numPosZones, numNegZones, deltaX, deltaY)
    """
    return [
        # Verticals
        ((0,   0), (1, 1, 1, 0)),
        ((0,  90), (0, 1, 1, 0)),
        ((0, 180), (1, 1, 1, 0)),
        ((0, 270), (1, 0, 1, 0)),

        # SW-NE
        ((45,   0), (1, 1, 1, 1)),
        ((45,  90), (0, 1, 1, 1)),
        ((45, 180), (1, 1, 1, 1)),
        ((45, 270), (1, 0, 1, 1)),

        # Horizontals
        ((90,   0), (1, 1, 0, 1)),
        ((90,  90), (0, 1, 0, 1)),
        ((90, 180), (1, 1, 0, 1)),
        ((90, 270), (1, 0, 0, 1)),

        # SE-NW
        ((135,   0), (1, 1, -1, 1)),
        ((135,  90), (0, 1, -1, 1)),
        ((135, 180), (1, 1, -1, 1)),
        ((135, 270), (1, 0, -1, 1)),
           ]


  def _findSuppressionRules(self, orientation, phase):
    """
    Return a set of algorithmic parameters for performing lobe suppression
    given the specified values of 'orientation' and 'phase'.
    """
    # (orientation, phase) ==> (numPosZones, numNegZones, deltaX, deltaY)
    phase = int(self._rad2deg(phase))
    orientation = int(self._rad2deg(orientation))
    return [x[1] for x in self._suppressionRules() if x[0][0] == orientation and x[0][1] == phase][0]


  def _doLobeSuppression(self, filter, orientation, phase):
    """
    Suppress (set to zero) the filter values outside of the main (primary)
    lobes for particular filter.
    Returns a copy of the filter with lobes suppressed.
    """
    # Obtain rules for scanning
    (numPosZones, numNegZones, deltaX, deltaY) = self._findSuppressionRules(orientation, phase)

    # Do the actual lobe-scanning
    sideLen = filter.shape[0]
    lobeFound = self._scanForLobes(filter, numPosZones, numNegZones, deltaX, deltaY,
      sideLen, self._rad2deg(phase))

    # Zero out the lobes found
    if lobeFound is not None:
      filter = self._zeroLobes(filter, lobeFound, (deltaX, deltaY))

    return filter


  def _zeroLobes(self, filter, (lobeX, lobeY), (deltaX, deltaY)):
    """
    Perform the actual suppression of the secondary lobes of a
    filter, assuming these secondary lobes have already been located
    and are identified by the 2-tuple (lobeX, lobeY).

    @param (lobeX lobeY) -- The position (relative to the filter's own
              local coordinate frame) where the lobe suppression should
              commence.
    @param (deltaX, deltaY) -- The direction (from filter center) in
              which the search for secondary lobes proceeded.
    """
    wipes = []
    sideLen = filter.shape[0]

    # Vertical scan ==> horizontal wipes
    if deltaX == 0 and deltaY > 0:
      for y in range(lobeY, sideLen):
        wipes += [(lobeX, y, 1, 0)]
        wipes += [(lobeX, y, -1, 0)]
      for y in range(sideLen - lobeY - 1, -1, -1):
        wipes += [(lobeX, y, 1, 0)]
        wipes += [(lobeX, y, -1, 0)]

    # Horizontal scan ==> vertical wipes
    elif deltaX > 0 and deltaY == 0:
      for x in range(lobeX, sideLen):
        wipes += [(x, lobeY, 0, 1)]
        wipes += [(x, lobeY, 0, -1)]
      for x in range(sideLen - lobeX - 1, -1, -1):
        wipes += [(x, lobeY, 0, 1)]
        wipes += [(x, lobeY, 0, -1)]

    # SW-bound scan ==> NW-SE wipes
    elif deltaX < 0 and deltaY > 0:
      for k in range(lobeY, sideLen):
        wipes += [(sideLen-k-1, k, -1, -1)]
        wipes += [(sideLen-k-1, k, 1, 1)]
        wipes += [(max(0, sideLen-k-2), k, -1, -1)]
        wipes += [(max(0, sideLen-k-2), k, 1, 1)]
      for k in range(sideLen - lobeY - 1, -1, -1):
        wipes += [(sideLen-k-1, k, -1, -1)]
        wipes += [(sideLen-k-1, k, 1, 1)]
        wipes += [(min(sideLen-1, sideLen-k), k, -1, -1)]
        wipes += [(min(sideLen-1, sideLen-k), k, 1, 1)]

    # SE-bound scan ==> SW-NE wipes
    elif deltaX > 0 and deltaY > 0:
      for k in range(lobeY, sideLen):
        wipes += [(k, k, -1, 1)]
        wipes += [(k, k, 1, -1)]
        wipes += [(min(sideLen-1, k+1), k, -1, 1)]
        wipes += [(min(sideLen-1, k+1), k, 1, -1)]
      for k in range(sideLen - lobeY - 1, -1, -1):
        wipes += [(k, k, -1, 1)]
        wipes += [(k, k, 1, -1)]
        wipes += [(max(0, k-1), k, -1, 1)]
        wipes += [(max(0, k-1), k, 1, -1)]

    # Do the wipes
    for wipe in wipes:
      filter = self._wipeStripe(filter, wipe)

    return filter


  def _wipeStripe(self, filter, (x, y, dX, dY)):
    """
    Zero out a particular row, column, or diagonal within the filter.

    @param filter -- The filter to be modified (a numpy array).
    @param x, y -- The starting point (in filter coordinate frame)
              of a stripe that is to be zeroed out.
    @param dX, dY -- The direction to proceed (and zero) until
              the edge of the filter is encountered.
    """
    sideLen = filter.shape[0]
    while True:
      filter[y, x] = 0.0
      x += dX
      y += dY
      if min(x, y) < 0 or max(x, y) >= sideLen:
        return filter


  def _scanForLobes(self, filter, numPosZones, numNegZones, deltaX, deltaY,
      sideLen, phase):
    """
    Search a filter for the location of it's secondary lobes, starting from
    the center of the filter.

    @param filter -- The filter to be searched for secondary lobe(s).
    @param numPosZones -- The number of regions of positive filter values that
              are contained within the primary lobes.
    @param numNegZones -- The number of regions of negative filter values that
              are contained within the primary lobes.
    @param deltaX, deltaY -- The direction in which to proceed (either vertically,
              horizontally, or diagonally) during the search.
    @param sideLen -- The length of one side of the filter (in pixels.)
    @param phase -- Phase of the filter (in degrees).
    """
    # Choose starting point
    x = sideLen/2
    y = sideLen/2

    if sideLen % 2 == 0:
      x += min(deltaX, 0)
      y += min(deltaY, 0)

    elif phase in (90, 270):
      # Odd-symmetric filters of odd dimension are 0 in the center
      # Skip the center point in order to not over-count lobes
      x += deltaX
      y += deltaY

    posZonesFound = 0
    negZonesFound = 0
    curZone = None
    while True:
      # Determine zone type
      maskVal = filter[y,x]
      if maskVal > 0.0 and curZone is not 'positive':
        curZone = 'positive'
        posZonesFound += 1
      elif maskVal < 0.0 and curZone is not 'negative':
        curZone = 'negative'
        negZonesFound += 1

      # Are we done?
      if posZonesFound > numPosZones or negZonesFound > numNegZones:
        return (x, y)
      else:
        x += deltaX
        y += deltaY

      # Hit the edge?
      if x == sideLen or x < 0 or y == sideLen or y < 0:
        return None


  def _doProcessing(self, imageData, maskData=None):
    """
    Apply the bank of pre-computed Gabor filters against a submitted image.
    @param imageData -- The image to be filtered.
    @param maskData -- Used if _wipeOutsideMask is True.
    """
    if self._debugMode:
      print 'GaborFilter._doProcessing(): imageData:'
      print imageData

    imageData *= (1.0/255.0)

    filterIndex = 0
    responseSet = []

    if self._convolutionMethod == '1D':
      # Create rotated version of image if necessary
      createRotated = False
      imageDataRotated = None
      for filterSet in self._gaborBank:
        for (filter, filterSpecs) in filterSet:
          if filterSpecs['orient'] in (45, 135):
            createRotated = True
            break
        if createRotated:
          break
      if createRotated:
        y, x = imageData.shape
        z = int(round((x+y) * 1/(2**.5)))
        imageDataRotated = numpy.zeros((z,z), dtype)
        self._rotation.rotate(imageData, imageDataRotated, y, x, z)

    if self._wipeOutsideMask:
      # Get all the scales
      scales = []
      for filter in self._gaborBank:
        if filter[0][1]['scale'] not in scales:
          scales.append(filter[0][1]['scale'])
      # Create an eroded mask for each scale
      erodedMaskData = {}
      for scale in scales:
        # Create a uniform filter at the specified scale
        filter = numpy.ones((scale*2 + 1), dtype)
        # Convolve the filter with the mask
        convolution = Float32SeparableConvolution2D()
        convolution.init(maskData.shape[0], maskData.shape[1],
          filter.shape[0], filter.shape[0], filter, filter)
        erodedMaskData[scale] = numpy.zeros(maskData.shape, dtype)
        convolution.compute(maskData, erodedMaskData[scale])

    for filterSet in self._gaborBank:
      filterResponse = []
      for (filter, filterSpecs) in filterSet:

        if self._debugMode:
          print 'Applying filter:  phase=%f  scale=%f  orient=%d' % (filterSpecs['phase'],
                                                                     filterSpecs['scale'],
                                                                     filterSpecs['orient'])

        # Perform the convolution
        if self._convolutionMethod == '2D':
          response = convolve2D(imageData, filter)

        elif self._convolutionMethod == '1D':
          response = convolve1D(imageData, imageDataRotated, filter,
            phase=filterSpecs['phase'], orientation=filterSpecs['orient'],
            rotation=self._rotation)

        elif self._convolutionMethod == 'FFT':
          response = convolveFFT(imageData, filter)

        else:
          raise RuntimeError("Unknown convolution method: "
            + self._convolutionMethod)

        if self._wipeOutsideMask:
          # Zero the response outside the mask
          mask = erodedMaskData[filterSpecs['scale']]
          maskMax = 255 * (filterSpecs['scale'] * 2 + 1) ** 2
          response[mask < maskMax] = 0

        maxResponse = response.max()
        if maxResponse > 0.0:
          response *= (1.0 / maxResponse)

        postProcessingMode = filterSpecs['mode']
        threshold          = filterSpecs['threshold']
        steepness          = filterSpecs['steepness']

        # Perform post-processing
        if postProcessingMode != 'raw':
          response = self._postProcess(response, postProcessingMode, threshold, steepness)

        filterResponse += [response]

      # Combine sequential filters to compute energy
      if len(filterResponse) > 1:
        if self._debugMode:
          print 'Computing combined energy...'
        combinedResponse = self._combineResponses(filterResponse)
      else:
        combinedResponse = filterResponse[0]

      responseSet += [combinedResponse]

    return responseSet


  def _combineResponses(self, responseSet):
    """
    Combine a list of one or more individual Gabor response maps
    into a single combined response map.
    Uses L2 norm to combine the responses.
    """
    combinedResponse = numpy.zeros(responseSet[0].shape)
    for response in responseSet:
      combinedResponse += (response * response)
    combinedResponse = numpy.sqrt(combinedResponse.clip(min=0.0, max=1.0))
    return combinedResponse


  def _postProcess(self, preResponse, postProcessingMode, threshold=0.0, steepness=500.0):
      """performs post-processing on the raw Gabor responses.

      Modes are as follows:
      'raw'       Perform no post-processing. Outputs will be in (-1.0, +1.0)
      'abs'       Output is absolute value of raw response.  Output will
                  lie in range (0.0, +1.0)
      'positive'  Clip negative raw response values to zero.  Output will lie
                  in range (0.0, +1.0)
      'rectified' Output is square of input.  Output will lie in the
                  range (0.0, +1.0)
      'power'     Positive raw responses are squared; negative raw responses
                  are clipped to zero.  Output will be in (0.0, +1.0)
      'hardFull'  Apply a hard threshold to the input; if raw response is
                  >= 'threshold', then output is 1.0, else output is -1.0.
                  Output lies in (-1.0, +1.0)
      'hardHalf'  Apply a hard threshold to the input; if raw response is
                  >= 'threshold', then output is 1.0, else output is 0.0.
                  Output lies in (0.0, +1.0)
      'sigmoidFull': Apply a sigmoidal threshold function to the input using
                  'threshold' and 'steepness'.  Output lies in (-1.0, +1.0).
      'sigmoidHalf': Apply a sigmoidal threshold function to the input using
                  'threshold' and 'steepness'.  Output lies in (0.0, +1.0).
      """

      # No processing
      if postProcessingMode == 'raw':
         postResponse = preResponse

      # Compute absolute value
      elif postProcessingMode == 'abs':
        postResponse = abs(preResponse)

      # Negative values set to 0.0
      elif postProcessingMode == 'positive':
        postResponse = preResponse.clip(min=0.0, max=1.0)

      # Compute square of response
      elif postProcessingMode == 'rectified':
        postResponse = preResponse * preResponse

      # Compute square of response for positive values
      elif postProcessingMode == 'power':
        intResponse = preResponse.clip(min=0.0, max=1.0)
        postResponse = intResponse * intResponse

      # Compute polynomial response for positive values
      elif postProcessingMode == 'polynomial':
        #intResponse = preResponse.clip(min=0.0, max=1.0)
        #postResponse = intResponse ** int(steepness)
        gain = 1.0 / ((1.0 - threshold) ** steepness)
        intResponse = (preResponse - threshold).clip(min=0.0, max=1.0)
        postResponse = gain * (intResponse ** steepness)

      # If output is > threshold, set to 1.0, else -1.0
      elif postProcessingMode == 'hardFull':
        postResponse = (preResponse > threshold).astype(dtype) * 2.0 - 1.0

      # If output is > threshold, set to 1.0, else 0.0
      elif postProcessingMode == 'hardHalf':
        postResponse = (preResponse > threshold).astype(dtype)

      # Sigmoid
      elif postProcessingMode == 'sigmoidHalf':
        postResponse = 1.0 / (numpy.exp(numpy.clip(steepness * (threshold - preResponse), -40.0, 40.0)) + 1.0)

      # Sigmoid
      elif postProcessingMode == 'sigmoidFull':
        postResponse = (2.0 / (numpy.exp(numpy.clip(steepness * (threshold - preResponse), -40.0, 40.0)) + 1.0)) - 1.0

      return postResponse


  def _printGaborFilter(self, g):
    """
    Print a Gabor filter mask in reasonably pretty format.
    @param g -- numpy array embodying the filter to be printed.
    """
    for j in range(g.shape[0]):
      for i in range(g.shape[1]):
        print '%7.4f' % g[j,i],
      print


  def _convertToPIL(self, responseSet, amplification=1.0):
    """
    Convert a list of gabor responses (represented as numpy arrays)
    into a list of PIL images.
    """

    imageSet = []
    for responseIndex in range(len(responseSet)):
      response = responseSet[responseIndex]
      maxPixelVal = 255.0
      halfMaxPixelVal = 0.5 * maxPixelVal
      # im = Image.new('L', (response.shape[1], response.shape[0]))
      # im.putdata(response.flatten(), scale=maxPixelVal * amplification, offset=0.0)
      im = Image.fromarray(((response*maxPixelVal*amplification).clip(min=0, max=255.0)).astype(numpy.uint8))
      imageSet += [im]
    return imageSet


def convolve2D(image, filter):
  """
  Convolve 2D filter with 2D image.
  """

  filterDim = filter.shape[0]
  # filterRadius = (filterDim - 1) / 2
  filterRadius = filterDim / 2
  flatFilter = filter.flatten()

  numPosnX = image.shape[1] - filterDim + 1
  numPosnY = image.shape[0] - filterDim + 1

  # response = scipy.signal.convolve(image, filter, mode='same')
  response = numpy.zeros(image.shape, dtype=dtype)
  for j in range(numPosnY):
    for i in range(numPosnX):
      response[j+filterRadius, i+filterRadius] = numpy.inner(flatFilter,
                                                 image[j:j+filterDim,
                                                 i:i+filterDim].flatten())

  return response

def convolve1D(image, imageRotated, filter, phase, orientation, rotation):
  """
  Convolve 2D filter with 2D image by approximating the filter as the outer
  product of two 1D filters and performing two separate convolutions.

  Results nearly match convolve2D if the filter is separable, with an average
  pixel intensity difference of about 1 / 1,000,000.
  """

  # Separate the 2D filter into two 1D filters
  if orientation in (0, 45):
    filterX = filter[filter.shape[0]/2,:].copy()
    if phase in (0, 180):
      filterY = filter[:,filter.shape[1]/2].copy()
    else:
      # Sinusoid is zero in the center
      filterY = filter[:,filter.shape[1]/2-1].copy()
  elif orientation in (90, 135):
    filterY = filter[:,filter.shape[1]/2].copy()
    if phase in (0, 180):
      filterX = filter[filter.shape[1]/2,:].copy()
    else:
      # Sinusoid is zero in the center
      filterX = filter[filter.shape[1]/2-1,:].copy()
  else:
    raise RuntimeError("convolve1D cannot do orientation %d" % orientation)

  if phase in (180, 270):
    if orientation in (0, 45):
      filterY *= -1
    elif orientation in (90, 135):
      filterX *= -1

  if orientation in (45, 135):
    imageToUse = imageRotated
  else:
    imageToUse = image

  # Process the 2D convolution as two 1D convolutions
  convolution = Float32SeparableConvolution2D()
  convolution.init(imageToUse.shape[0], imageToUse.shape[1],
    filterX.shape[0], filterY.shape[0], filterX, filterY)
  response = numpy.zeros(imageToUse.shape, dtype)
  convolution.compute(imageToUse, response)

  if orientation in (45, 135):
    # The image has been rotated
    # Unrotate the image
    y, x = image.shape
    z = response.shape[0]
    unrotatedResponse = numpy.zeros((y, x), dtype)
    rotation.unrotate(unrotatedResponse, response, y, x, z)
    response = unrotatedResponse
    # Zero the edges
    response[:filter.shape[0]/2,:] = 0.0
    response[:,:filter.shape[1]/2] = 0.0
    response[-filter.shape[0]/2+1:,:] = 0.0
    response[:,-filter.shape[1]/2+1:] = 0.0

  return response

def convolveFFT(image, filter):
  """
  Convolve 2D filter with 2D image using FFT -> multiply -> IFFT.

  Results nearly match convolve2D, with an average pixel intensity difference
  of about 1 / 1,000,000.
  """

  size = [
    pow(2,int(math.ceil(math.log(image.shape[0] + filter.shape[0] - 1, 2)))),
    pow(2,int(math.ceil(math.log(image.shape[1] + filter.shape[1] - 1, 2))))
  ]

  image2 = numpy.zeros(size)
  image2[0:image.shape[0], 0:image.shape[1]] = image
  image2 = numpy.fft.fft2(image2)

  filter = numpy.fliplr(numpy.flipud(filter))
  filter2 = numpy.zeros(size)
  filter2[0:filter.shape[0], 0:filter.shape[1]] = filter
  filter2 = numpy.fft.fft2(filter2)

  response = numpy.fft.ifft2(image2 * filter2)

  x = (filter.shape[0] - 1) / 2
  y = (filter.shape[1] - 1) / 2
  response = response[x:x+image.shape[0], y:y+image.shape[1]]
  response[:filter.shape[0]/2,:] = 0.0
  response[:,:filter.shape[1]/2] = 0.0
  response[-filter.shape[0]/2+1:,:] = 0.0
  response[:,-filter.shape[1]/2+1:] = 0.0

  return response
