gaussian_filter Subroutine

private subroutine gaussian_filter(long, larg, xc, gauss_filt)

Gaussian kernel

Arguments

Type IntentOptional Attributes Name
integer(kind=I4), intent(in) :: long

2D array length
integer(kind=I4), intent(in) :: larg

2D array width
real(kind=R8), intent(in) :: xc

the cut-off wavelength
real(kind=R8), intent(out), dimension(1:long, 1:larg) :: gauss_filt

2D array out

Called by

proc~~gaussian_filter~~CalledByGraph proc~gaussian_filter gaussian_filter proc~fft_filter fft_filter proc~fft_filter->proc~gaussian_filter proc~multiple_anisotropy multiple_anisotropy proc~multiple_anisotropy->proc~fft_filter proc~test_peaks_and_pits_curvatures test_peaks_and_pits_curvatures proc~test_peaks_and_pits_curvatures->proc~fft_filter program~test_smooth test_smooth program~test_smooth->proc~fft_filter program~test_anisotropy test_anisotropy program~test_anisotropy->proc~multiple_anisotropy program~test_grad_curv test_grad_curv program~test_grad_curv->proc~test_peaks_and_pits_curvatures

Source Code

   subroutine gaussian_filter(long, larg, xc, gauss_filt)
   !================================================================================================
   !! Gaussian kernel
   !------------------------------------------------------------------------------------------------
   implicit none
   integer(kind=I4), intent(in )                            :: long        !! *2D array length*
   integer(kind=I4), intent(in )                            :: larg        !! *2D array width*
   real   (kind=R8), intent(in )                            :: xc          !! *the cut-off wavelength*
   real   (kind=R8), intent(out), dimension(1:long, 1:larg) :: gauss_filt  !! *2D array out*

      integer(kind=I4) :: i, j
      real   (kind=R8) :: tmp, xi, xj

      real   (kind=R8), parameter :: const = sqrt( log(2._R8)/PI_R8 )

      do j = 2, larg/2 +1
      do i = 2, long/2 +1
         xi = (i-1) ; xj = (j-1)
         xi = xi/(long -1) ; xj = xj/(larg -1)
         tmp = gaussian_function(xi, xj, xc)
         gauss_filt(       +i,        +j) = tmp
         gauss_filt(long+2 -i,        +j) = tmp
         gauss_filt(       +i, larg+2 -j) = tmp
         gauss_filt(long+2 -i, larg+2 -j) = tmp
      enddo
      enddo
      do j = 2, larg/2 +1
         i = 1
         xi = (i-1) ; xj = (j-1)
         xi = xi/(long -1) ; xj = xj/(larg -1)
         tmp = gaussian_function(xi, xj, xc)
         gauss_filt(i,         j) = tmp
         gauss_filt(i, larg+2 -j) = tmp
      enddo
      do i = 2, long/2 +1
         j = 1
         xi = (i-1) ; xj = (j-1)
         xi = xi/(long -1) ; xj = xj/(larg -1)
         tmp = gaussian_function(xi, xj, xc)
         gauss_filt(        i, j) = tmp
         gauss_filt(long+2 -i, j) = tmp
      enddo
      i = 1
      j = 1
      xi = (i-1) ; xj = (j-1)
      xi = xi/(long -1) ; xj = xj/(larg -1)
      gauss_filt(i, j) = gaussian_function(xi, xj, xc)

   contains
      !-----------------------------------------
      real(kind=R8) function gaussian_function(xi, xj, xc)
      implicit none
      real(kind=R8), intent(in) :: xi
      real(kind=R8), intent(in) :: xj
      real(kind=R8), intent(in) :: xc ! fréquence de coupure, plus exactement proportion : (freq coup) / (nb points)

         gaussian_function = exp( -PI_R8 * const**2 * ( xi**2 + xj**2 ) / (xc**2) ) ! sqrt(ln(2)/pi)=0.47

      return
      endfunction gaussian_function
      !-----------------------------------------
   endsubroutine gaussian_filter