title: Implementing log Gabor filter bank
body: <p>I was reading this paper <a href="http://staff.utia.cas.cz/sroubekf/papers/gabor_07.pdf">"Self-Invertible 2D Log-Gabor Wavelets"</a> it defines 2D log gabor filter as such:</p>
  
  <p><a href="https://i.stack.imgur.com/dAUHE.png"><img src="https://i.stack.imgur.com/dAUHE.png" alt="enter image description here"></a>
  <a href="https://i.stack.imgur.com/tuL08.png"><img src="https://i.stack.imgur.com/tuL08.png" alt="enter image description here"></a></p>
  
  <p>The paper also states that the filter only covers one side of the frequency space and shows that in this image</p>
  
  <p><a href="https://i.stack.imgur.com/4cIS3.png"><img src="https://i.stack.imgur.com/4cIS3.png" alt="enter image description here"></a></p>
  
  <p>On my attempt to implement the filter I get results that do not match with what is said in the paper. Let me start with my implementation then I will state the problems.</p>
  
  <p><strong>Implementation:</strong></p>
  
  <ol>
  <li><p>I created a 2d array that  contains the filter and transformed each index so that the origin of the frequency domain is at the center of the array with positive x-axis going right and positive y-axis going up.</p>
  
  <pre class="lang-py prettyprint-override"><code>number_scales = 5         # scale resolution
  number_orientations = 9   # orientation resolution
  N = constantDim           # image dimensions
  
  def getLogGaborKernal(scale, angle, logfun=math.log2, norm = True):
      # setup up filter configuration
      center_scale = logfun(N) - scale          
      center_angle = ((np.pi/number_orientations) * angle) if (scale % 2) \
                  else ((np.pi/number_orientations) * (angle+0.5))
      scale_bandwidth =  0.996 * math.sqrt(2/3)
      angle_bandwidth =  0.996 * (1/math.sqrt(2)) * (np.pi/number_orientations)
  
      # 2d array that will hold the filter
      kernel = np.zeros((N, N))
      # get the center of the 2d array so we can shift origin
      middle = math.ceil((N/2)+0.1)-1
  
      # calculate the filter
      for x in range(0,constantDim):
          for y in range(0,constantDim):
              # get the transformed x and y where origin is at center
              # and positive x-axis goes right while positive y-axis goes up
              x_t, y_t = (x-middle),-(y-middle)
              # calculate the filter value at given index
              kernel[y,x] = logGaborValue(x_t,y_t,center_scale,center_angle,
          scale_bandwidth, angle_bandwidth,logfun)
  
      # normalize the filter energy
      if norm:
          Kernel = kernel / np.sum(kernel**2)
      return kernel
  </code></pre></li>
  <li><p>To calculate the filter value at each index another transform is made where we go to the log-polar space</p>
  
  <pre class="lang-py prettyprint-override"><code>def logGaborValue(x,y,center_scale,center_angle,scale_bandwidth,
                angle_bandwidth, logfun):
      # transform to polar coordinates
      raw, theta = getPolar(x,y)
      # if we are at the center, return 0 as in the log space
      # zero is not defined
      if raw == 0:
          return 0
  
      # go to log polar coordinates
      raw = logfun(raw)
  
      # calculate (theta-center_theta), we calculate cos(theta-center_theta) 
      # and sin(theta-center_theta) then use atan to get the required value,
      # this way we can eliminate the angular distance wrap around problem
      costheta, sintheta = math.cos(theta), math.sin(theta)
      ds = sintheta * math.cos(center_angle) - costheta * math.sin(center_angle)    
      dc = costheta * math.cos(center_angle) + sintheta * math.sin(center_angle)  
      dtheta = math.atan2(ds,dc)
  
      # final value, multiply the radial component by the angular one
      return math.exp(-0.5 * ((raw-center_scale) / scale_bandwidth)**2) * \
              math.exp(-0.5 * (dtheta/angle_bandwidth)**2)
  </code></pre></li>
  </ol>
  
  <p><strong>Problems:</strong></p>
  
  <ol>
  <li><p><strong>The angle:</strong> the paper stated that indexing the angles from 1->8 would produce good coverage of the orientation, but in my implementation angles from 1->n don't cover except for half orientations. Even the vertical orientation is not correctly covered. This can be shown in this figure which contains sets of filters of scale 3 and orientations ranging from 1->8:</p>
  
  <p><a href="https://i.stack.imgur.com/fJz4o.png"><img src="https://i.stack.imgur.com/fJz4o.png" alt="enter image description here"></a></p></li>
  <li><p><strong>The coverage:</strong> from filters above it is clear the filter covers both sides of the space which is not what the paper says. This can be made more explicit by using 9 orientations ranging from -4 -> 4. The following image contains all the filters in one image to show how it covers both sides of the spectrum (this image is created by taking the maximum at each location from all filters):</p>
  
  <p><a href="https://i.stack.imgur.com/gNAT3.png"><img src="https://i.stack.imgur.com/gNAT3.png" alt="enter image description here"></a></p></li>
  <li><p><strong>Middle Column (orientation $\pi / 2$):</strong> in the first figure in orientation from 3 -> 8 it can be seen that the filter vanishes at orientation $ \pi / 2$. Is this normal? This can be seen too when I combine all the filters(of all 5 scales and 9 orientations) in one image:</p>
  
  <p><a href="https://i.stack.imgur.com/BVk49.png"><img src="https://i.stack.imgur.com/BVk49.png" alt="enter image description here"></a></p></li>
  </ol>
  
  <p><strong>Update:</strong>
  Adding the impulse response of the filter in spatial domain, as you can see there is an obvious distortion in -4 &amp; 4 orientations:</p>
  
  <p><a href="https://i.stack.imgur.com/gl1nY.png"><img src="https://i.stack.imgur.com/gl1nY.png" alt="enter image description here"></a></p>