median_filter Subroutine

public subroutine median_filter(tab, long, larg, snb, kernel, sig, omp)

A bit more complex filter: the overall height standard deviation is taken into account

Arguments

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

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

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

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

patch number along a direction
integer(kind=I4), intent(in) :: kernel

kernel size
real(kind=R8), intent(in) :: sig

error std
logical(kind=I4), intent(in) :: omp

if multithreading

Calls

proc~~median_filter~~CallsGraph proc~median_filter median_filter proc~calc_median calc_median proc~median_filter->proc~calc_median proc~calc_moments calc_moments proc~median_filter->proc~calc_moments proc~median_smooth median_smooth proc~median_filter->proc~median_smooth sort_array2 sort_array2 proc~calc_median->sort_array2 proc~calc_moments_1d calc_moments_1D proc~calc_moments->proc~calc_moments_1d proc~median_smooth->proc~calc_median omp_get_num_procs omp_get_num_procs proc~median_smooth->omp_get_num_procs proc~median_smooth->sort_array2

Called by

proc~~median_filter~~CalledByGraph proc~median_filter median_filter program~test_smooth test_smooth program~test_smooth->proc~median_filter

Source Code

   subroutine median_filter(tab, long, larg, snb, kernel, sig, omp)
   !================================================================================================
   !! A bit more complex filter: the overall height standard deviation is taken into account
   !------------------------------------------------------------------------------------------------
   implicit none
   integer(kind=I4), intent(in   )                            :: long   !! *2D array length*
   integer(kind=I4), intent(in   )                            :: larg   !! *2D array width*
   integer(kind=I4), intent(in   )                            :: snb    !! *patch number along a direction*
   integer(kind=I4), intent(in   )                            :: kernel !! *kernel size*
   logical(kind=I4), intent(in   )                            :: omp    !! *if multithreading*
   real   (kind=R8), intent(in   )                            :: sig    !! *error std*
   real   (kind=R8), intent(inout), dimension(1:long, 1:larg) :: tab    !! *2D array*

      integer(kind=I4) :: i, j, k
      real(kind=R8)    :: md

      integer(kind=I4), dimension(1:snb + 1) :: li, lj

      real(kind=R8), dimension(1:snb * snb)  :: ect

      real(kind=R8), allocatable, dimension(:,:) :: tab_tmp1, tab_tmp2

      type( moment_stat ) :: mx_smooth

      allocate( tab_tmp1(1:long, 1:larg) )
      allocate( tab_tmp2(1:long, 1:larg) )

      ! first determine the difference between the input surface and a median filtered one
      tab_tmp1(1:long, 1:larg) = tab(1:long, 1:larg)
      call median_smooth(  tab    = tab_tmp1(1:long, 1:larg),  &  !
                           kernel = kernel,                    &  !
                           long   = long,                      &  !
                           larg   = larg,                      &  !
                           omp    = omp )                         !
      tab_tmp2(1:long, 1:larg) = tab(1:long, 1:larg) - tab_tmp1(1:long, 1:larg)

      ! bounds when patching domain
      li(1) = 1 ; li(snb + 1) = long
      lj(1) = 1 ; lj(snb + 1) = larg

      do i = 2, snb
         li(i) = li(i - 1) + int( real(long, kind=R8)/snb, kind=I4 )
         lj(i) = lj(i - 1) + int( real(larg, kind=R8)/snb, kind=I4 )
      enddo

      k  = 0
      do j = 1, snb
      do i = 1, snb

         call calc_moments(   tab = reshape( tab_tmp2( li(i):li(i + 1) - 1, lj(j):lj(j + 1) - 1 ),    &  !
                                             [( li(i) - li(i + 1) ) * ( lj(j) - lj(j + 1) )] ),       &  !
                               mx = mx_smooth,                                                        &  !
                           nb_mom = 2 )                                                                  !
         k = k +1
         ect(k) = mx_smooth%si

      enddo
      enddo

      call calc_median( tab = ect( 1:snb*snb ), &  !
                         md = md )                 !

      call calc_moments(   tab = reshape( tab_tmp2(1:long, 1:larg),  &  !
                                          [long * larg] ),           &  !
                            mx = mx_smooth,                          &  !
                        nb_mom = 2 )                                    !

      where( abs(tab_tmp2(1:long, 1:larg) -mx_smooth%mu) > sig*md ) tab(1:long, 1:larg) = tab_tmp1(1:long, 1:larg)

      deallocate( tab_tmp1, tab_tmp2 )

   return
   endsubroutine median_filter