top_hat_filter Subroutine

private subroutine top_hat_filter(long, larg, xc, top_filt)

Top-hat 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) :: top_filt

2D array out

Called by

proc~~top_hat_filter~~CalledByGraph proc~top_hat_filter top_hat_filter proc~fft_filter fft_filter proc~fft_filter->proc~top_hat_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 top_hat_filter(long, larg, xc, top_filt)
   !================================================================================================
   !! Top-hat 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) :: top_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 = top_hat_function(xi, xj, xc)
         top_filt(       +i,        +j) = tmp
         top_filt(long+2 -i,        +j) = tmp
         top_filt(       +i, larg+2 -j) = tmp
         top_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 = top_hat_function(xi, xj, xc)
         top_filt(i,         j) = tmp
         top_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 = top_hat_function(xi, xj, xc)
         top_filt(        i, j) = tmp
         top_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)
      top_filt(i, j) = top_hat_function(xi, xj, xc)

   contains
      !-----------------------------------------
      real(kind=R8) function top_hat_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)

         real(kind=R8) :: val

         val = (xi / xc) **2 + (xj / (xc * (long - 1) / (larg - 1))) **2

         top_hat_function = 0
         if ( val < 1. ) top_hat_function = 1

      return
      endfunction top_hat_function
      !-----------------------------------------
   endsubroutine top_hat_filter