pythonimage-processingcomputer-vision# Implementing log Gabor filter bank

I was reading this paper "Self-Invertible 2D Log-Gabor Wavelets" it defines 2D log gabor filter as such:

The paper also states that the filter only covers one side of the frequency space and shows that in this image

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.

**Implementation:**

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.

`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`

To calculate the filter value at each index another transform is made where we go to the log-polar space

`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)`

**Problems:**

**The angle:**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:**The coverage:**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):**Middle Column (orientation $\pi / 2$):**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:

**Update:**
Adding the impulse response of the filter in spatial domain, as you can see there is an obvious distortion in -4 & 4 orientations:

Solution

After a lot of code analysis, I found that my implementation was correct but the `getPolar`

function was messed up, so the code above should work just fine. This is the a new code without the `getPolar`

function if any one was looking for it:

```
number_scales = 5 # scale resolution
number_orientations = 8 # orientation resolution
N = 128 # image dimensions
def getFilter(f_0, theta_0):
# filter configuration
scale_bandwidth = 0.996 * math.sqrt(2/3)
angle_bandwidth = 0.996 * (1/math.sqrt(2)) * (np.pi/number_orientations)
# x,y grid
extent = np.arange(-N/2, N/2 + N%2)
x, y = np.meshgrid(extent,extent)
mid = int(N/2)
## orientation component ##
theta = np.arctan2(y,x)
center_angle = ((np.pi/number_orientations) * theta_0) if (f_0 % 2) \
else ((np.pi/number_orientations) * (theta_0+0.5))
# 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 = np.cos(theta)
sintheta = np.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 = np.arctan2(ds,dc)
orientation_component = np.exp(-0.5 * (dtheta/angle_bandwidth)**2)
## frequency componenet ##
# go to polar space
raw = np.sqrt(x**2+y**2)
# set origin to 1 as in the log space zero is not defined
raw[mid,mid] = 1
# go to log space
raw = np.log2(raw)
center_scale = math.log2(N) - f_0
draw = raw-center_scale
frequency_component = np.exp(-0.5 * (draw/ scale_bandwidth)**2)
# reset origin to zero (not needed as it is already 0?)
frequency_component[mid,mid] = 0
return frequency_component * orientation_component
```

- Python Jinja2 LaTeX Table
- Getting attributes of a class
- How can I print many significant figures in Python?
- How to allow list append() method to return the new list
- Calculate Last Friday of Month in Pandas
- Python type hint for Iterable[str] that isn't str
- How to iterate over a list in chunks
- How to exit the entire application from a Python thread?
- Running shell command and capturing the output
- How do I pass a variable by reference?
- Convert range(r) to list of strings of length 2 in python
- How can I get the start and end dates for each week?
- how to use send_message() in python-telegram-bot
- Python conditional replacement based on element type
- How can I count the number of items in an arbitrary iterable (such as a generator)?
- Find longest consecutive range of numbers in list
- Insert text in braces with asyncpg
- How does one put a link / url to the web-site's home page in Django?
- How to determine if a path is a subdirectory of another?
- Custom Keybindings for Ipython terminal
- FastAPI asynchronous background tasks blocks other requests?
- How to make sure that information from one file is duplicated into several text documents, without specific lines
- Installing a Python environment with Anaconda
- sklearn pipeline model predicting same results for all input
- Brew command not found after installing Anaconda Python
- How to get an XPath from selenium webelement or from lxml?
- Pipe PuTTY console to Python script
- How to align the axes of a figure in matplotlib?
- Persist ParentDocumentRetriever of langchain
- How to reset index in a pandas dataframe?