fft_filter Subroutine

public subroutine fft_filter(tab, long, larg, cutoff, bf_tab, multi_fft, pad, ext, type_apo, shift)

Classical Gaussian filter

Arguments

Type IntentOptional Attributes Name
real(kind=R8), intent(in), dimension(1:long, 1:larg) :: tab

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

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

2D array height
real(kind=R8), intent(in) :: cutoff

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

2D array out
logical(kind=I4), intent(in) :: multi_fft

multiple fft at once ?
real(kind=R8), intent(in), optional :: pad

fft padding
character(len=*), intent(in), optional :: ext

extension
character(len=*), intent(in), optional :: type_apo

apodization type
real(kind=R8), intent(in), optional, dimension(2) :: shift

surface shift fraction along x and y

Calls

proc~~fft_filter~~CallsGraph proc~fft_filter fft_filter calc_fftw3 calc_fftw3 proc~fft_filter->calc_fftw3 extend extend proc~fft_filter->extend proc~gaussian_filter gaussian_filter proc~fft_filter->proc~gaussian_filter proc~shifting shifting proc~fft_filter->proc~shifting proc~top_hat_filter top_hat_filter proc~fft_filter->proc~top_hat_filter tab_calc_fftw3 tab_calc_fftw3 proc~fft_filter->tab_calc_fftw3

Called by

proc~~fft_filter~~CalledByGraph proc~fft_filter fft_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 fft_filter(tab, long, larg, cutoff, bf_tab, multi_fft, pad, ext, type_apo, shift)
   !================================================================================================
   !! Classical Gaussian filter
   !------------------------------------------------------------------------------------------------
   implicit none
   integer(kind=I4), intent(in ) :: long                                !! *2D array width*
   integer(kind=I4), intent(in ) :: larg                                !! *2D array height*
   real   (kind=R8), intent(in ) :: cutoff                              !! *cut-off wavelength*
   logical(kind=I4), intent(in ) :: multi_fft                           !! *multiple fft at once ?*
   real   (kind=R8), intent(in ), optional :: pad                       !! *fft padding*
   character(len=*), intent(in ), optional :: ext                       !! *extension*
   character(len=*), intent(in ), optional :: type_apo                  !! *apodization type*
   real   (kind=R8), intent(in ), dimension(2), optional    :: shift    !! *surface shift fraction along x and y*
   real   (kind=R8), intent(in ), dimension(1:long, 1:larg) :: tab      !! *2D array in*
   real   (kind=R8), intent(out), dimension(1:long, 1:larg) :: bf_tab   !! *2D array out*

      integer(kind=I4) :: i, j, nx2, ny2, iex, iey, ibx, iby
      real   (kind=R8) :: o_pad
      logical(kind=I4) :: with_pad

      character(len=:), allocatable :: o_ext
      character(len=:), allocatable :: o_type_apo

      complex(kind=R8) :: cdx, cdy

      real   (kind=R8), dimension(:,:), allocatable :: tab_ext, gauss_tab, top_tab

      complex(kind=R8), dimension(:,:), allocatable :: cmpl1, cmpl2, shift_tab

      with_pad = .true.

      if ( .not.present( ext ) ) then

         o_ext = 'symmetry'

      else

         o_ext = ext

      endif

      if ( .not.present( type_apo ) ) then

         o_type_apo = 'no_apo'

      else

         o_type_apo = type_apo

      endif

      if ( .not.present( pad ) ) then

         ! default padding
         o_pad = PAD_FFT_FILTER

      else

         ! no padding
         if ( pad < 0. ) then

            o_pad = 1

            nx2 = long
            ny2 = larg

            with_pad = .false.

         else

            ! padding with argument
            o_pad = pad

         endif

      endif

      if ( with_pad ) then

         ! make even surface dimensions

         nx2 = 2 * ( nint(o_pad * long)/2 )
         ny2 = 2 * ( nint(o_pad * larg)/2 )

      endif

      allocate( shift_tab(1:nx2, 1:ny2) )

      shift_tab(1:nx2, 1:ny2) = 1

      if ( present( shift ) ) then

         ! shift surface with frequencies rotation
         call shifting(nx2, ny2, shift, shift_tab)

      endif

      allocate( tab_ext(1:nx2, 1:ny2)  )     !

      allocate( cmpl1(1:nx2, 1:ny2),      &  !
                cmpl2(1:nx2, 1:ny2) )        !

      if ( nx2 > long ) then

         ibx = ceiling( (nx2 - long)/2. ) ; iex = ibx + long - 1
         iby = ceiling( (ny2 - larg)/2. ) ; iey = iby + larg - 1

         call extend(   tab_in = tab(1:long, 1:larg),       &  !
                       tab_out = tab_ext(1:nx2, 1:ny2),     &  !
                            nx = long,                      &  !
                            ny = larg,                      &  !
                           nx2 = nx2,                       &  !
                           ny2 = ny2,                       &  !
                           ext = o_ext,                     &  !
                      type_apo = o_type_apo )                  !

      else

         tab_ext(1:nx2, 1:ny2) = tab(1:long, 1:larg)

      endif

      cmpl1(1:nx2, 1:ny2) = cmplx( tab_ext(1:nx2, 1:ny2), 0, kind = R8 )

      if (multi_fft) then

         call tab_calc_fftw3(   sens = FORWARD,                  &  !
                              tab_in = cmpl1(1:nx2, 1:ny2),      &  !
                              tab_ou = cmpl2(1:nx2, 1:ny2),      &  !
                                long = nx2,                      &  !
                                larg = ny2)                         !

      else

         call calc_fftw3(   sens = FORWARD,                    &  !
                          tab_in = cmpl1(1:nx2, 1:ny2),        &  !
                          tab_ou = cmpl2(1:nx2, 1:ny2),        &  !
                            long = nx2,                        &  !
                            larg = ny2)                           !

      endif

      if ( cutoff < 0. ) then

         ! TOP HAT FILTER
         allocate( top_tab(1:nx2, 1:ny2) )

         call top_hat_filter(  long     = nx2,                     &  !
                               larg     = ny2,                     &  !
                               xc       = -cutoff,                 &  !
                               top_filt = top_tab(1:nx2, 1:ny2) )     !

         cmpl1(1:nx2, 1:ny2) = cmpl2(1:nx2, 1:ny2) * top_tab(1:nx2, 1:ny2) * shift_tab(1:nx2, 1:ny2)

         deallocate( top_tab )

      else

         ! GAUSSIAN FILTER
         allocate( gauss_tab(1:nx2, 1:ny2) )

         call gaussian_filter( long       = nx2,                     &  !
                               larg       = ny2,                     &  !
                               xc         = cutoff,                  &  !
                               gauss_filt = gauss_tab(1:nx2, 1:ny2) )   !

         cmpl1(1:nx2, 1:ny2) = cmpl2(1:nx2, 1:ny2) * gauss_tab(1:nx2, 1:ny2)

         deallocate( gauss_tab )

      endif

      if (multi_fft) then

         call tab_calc_fftw3(   sens = BACKWARD,               &  !
                              tab_in = cmpl1(1:nx2, 1:ny2),    &  !
                              tab_ou = cmpl2(1:nx2, 1:ny2),    &  !
                                long = nx2,                    &  !
                                larg = ny2)                       !

      else

         call calc_fftw3(   sens = BACKWARD,                   &  !
                          tab_in = cmpl1(1:nx2, 1:ny2),        &  !
                          tab_ou = cmpl2(1:nx2, 1:ny2),        &  !
                            long = nx2,                        &  !
                            larg = ny2)                           !

      endif

      tab_ext(1:nx2, 1:ny2) = real(cmpl2(1:nx2, 1:ny2), kind=R8)

      if ( nx2 > long ) then

         bf_tab(1:long, 1:larg) = tab_ext(ibx:iex, iby:iey)

      else

         bf_tab(1:long, 1:larg) = tab_ext(1:nx2, 1:ny2)

      endif

      deallocate(cmpl1, cmpl2, shift_tab)
      deallocate(tab_ext)

   return
   endsubroutine fft_filter