program apply_gabor_filter

implicit none

integer, parameter :: pr = selected_real_kind(15,3)

integer :: i,j,k,l,m,n

real(pr) :: pi = dacos(-1.0d0)

real(pr) :: theta,f,sigma_x,sigma_y
integer :: x1,x2
real(pr) :: G

real(pr), dimension(512,512) :: IMG,CIMG
real(pr), dimension(:,:), allocatable :: GM

!----------------------------------------------------------------------------------------!

open(1,file='lena_gray.txt')
do i = 1, 512
end do
close(1)

!----------------------------------------------------------------------------------------!
! User Inputs

theta = 90.0 ! degrees

f = 0.1

sigma_x = 6
sigma_y = 6

!----------------------------------------------------------------------------------------!
! Create Gabor filter

theta = ( pi * theta  ) / 180.0 ! converting to radian

k = 1
l = 1
call Gabor_Filter(x1,x2,theta,f,sigma_x,sigma_y,G)
GM(k,l) = G
l = l + 1
end do
k = k + 1
end do

open(1,file='gabor_filter_01.txt')
end do
close(1)

!----------------------------------------------------------------------------------------!
! Apply Gabor Filter

CIMG = 0

do i = 1, 512
do j = 1, 512
k = 1
l = 1
m = i - x1
n = j - x2
if( 0 < m .AND. m <= 512)then
if( 0 < n .AND. n <= 512)then
CIMG(i,j) = CIMG(i,j) + GM(k,l) * IMG(m,n)
end if
end if
l = l + 1
end do
k = k + 1
end do
end do
end do

!----------------------------------------------------------------------------------------!

open(1,file='image_gabor_filter_01.txt')
do i = 1, 512
write(1,*) (CIMG(i,j),j=1,512)
end do
close(1)

!----------------------------------------------------------------------------------------!

end program apply_gabor_filter

!----------------------------------------------------------------------------------------!
! Gabor Filter Rotation

subroutine Change_Base(x1,x2,theta,x_theta,y_theta)

implicit none

integer, parameter :: pr = selected_real_kind(15,3)

integer :: x1
integer :: x2
real(pr) :: theta
real(pr) :: x_theta
real(pr) :: y_theta

x_theta = x1 * dcos(theta) + x2 * dsin(theta)
y_theta = x2 * dcos(theta) - x1 * dsin(theta)

end subroutine Change_Base

!----------------------------------------------------------------------------------------!
! Gabor Filter

subroutine Gabor_Filter(x1,x2,theta,f,sigma_x,sigma_y,G)

implicit none

integer, parameter :: pr = selected_real_kind(15,3)

integer :: x1
integer :: x2
real(pr) :: theta
real(pr) :: f
real(pr) :: sigma_x
real(pr) :: sigma_y

real(pr) :: x_theta
real(pr) :: y_theta

real(pr) :: r1,r2,arg

real(pr) :: pi = dacos(-1.0d0)
real(pr) :: G

call Change_Base(x1,x2,theta,x_theta,y_theta)

r1 = x_theta / sigma_x
r2 = y_theta / sigma_y

arg = - 0.5 * ( r1**2 + r2**2 )

G = dexp(arg) * dcos(2.0*pi*f*x_theta)

end subroutine Gabor_Filter

!----------------------------------------------------------------------------------------!