program test_smooth
use data_arch, only : I4, R8
use miscellaneous, only : get_unit
use surfile, only : read_surf, write_surf, init_scal, SCALE_SURF, unit2IUf
use filter, only : median_smooth, median_filter, soften, morpho_filter, fft_filter, PAD_FFT_FILTER
use fftw3, only : fftw_plan_with_nthreads, init_fftw3, end_fftw3, extend
!$ use omp_lib
implicit none
real(kind=R8), allocatable, dimension(:,:) :: heights, heights_out, heights_copy, mask, ext_heights
real(kind=R8) :: dx, dy, dz, fft_cutoff, pad
integer(kind=I4) :: nx, ny, nx2, ny2, n_th
type(SCALE_SURF) :: scal_surf, scal_mask
!========================================================================= padding / windowing
call read_surf( nom_fic = "sur/AB-Zinv-ART-8-21-200x200.sur", & ! IN
tab_s = heights, & ! OUT
scal = scal_surf ) ! OUT
nx = scal_surf%xres
ny = scal_surf%yres
pad = 1.5!PAD_FFT_FILTER
nx2 = 2 * ( nint(pad * nx)/2 )
ny2 = 2 * ( nint(pad * ny)/2 )
allocate( ext_heights(1:nx2, 1:ny2) )
call extend( tab_in = heights(1:nx, 1:ny), & !
tab_out = ext_heights(1:nx2, 1:ny2), & !
nx = nx, & !
ny = ny, & !
nx2 = nx2, & !
ny2 = ny2, & !
ext = 'constant', & !
type_apo = 'hann__') !
scal_surf%xres = nx2
scal_surf%yres = ny2
call write_surf( nom_fic = "out/test_extend.sur", & !
tab_s = ext_heights(1:nx2, 1:ny2), & !
scal = scal_surf ) !
deallocate( ext_heights )
!========================================================================= fft_filter 5 µm
call read_surf( nom_fic = "sur/AB-Zinv-ART-8-21-200x200.sur", & ! IN
tab_s = heights, & ! OUT
scal = scal_surf ) ! OUT
nx = scal_surf%xres
ny = scal_surf%yres
dx = scal_surf%dx * unit2IUf(scal_surf%dx_unit)
dy = scal_surf%dy * unit2IUf(scal_surf%dy_unit)
dz = 0
write(*,*) 'nx, ny = ', nx, ny
write(*,*) 'dx, dy = ', dx, dy
allocate( heights_out(1:nx, 1:ny) )
n_th = omp_get_max_threads()
call fftw_plan_with_nthreads(nthreads = n_th)
nx2 = 2 * ( nint(PAD_FFT_FILTER * nx)/2 )
ny2 = 2 * ( nint(PAD_FFT_FILTER * ny)/2 )
fft_cutoff = dx / 5.e-6
call init_fftw3(long = nx2, & !
larg = ny2 ) !
call fft_filter(tab = heights(1:nx, 1:ny), & !
long = nx, & !
larg = ny, & !
cutoff = fft_cutoff, & !
bf_tab = heights_out(1:nx, 1:ny), & !
multi_fft = .FALSE.) !
call write_surf( nom_fic = "out/test_fft_filter_005µm.sur", & !
tab_s = heights_out(1:nx, 1:ny), & !
scal = scal_surf ) !
write(*,*) 'gaussian filter cutoff = 5 µm'
!========================================================================= fft_filter 80 µm
fft_cutoff = dx / 80.e-6
call fft_filter(tab = heights(1:nx, 1:ny), & !
long = nx, & !
larg = ny, & !
cutoff = fft_cutoff, & !
bf_tab = heights_out(1:nx, 1:ny), & !
multi_fft = .FALSE.) !
call write_surf( nom_fic = "out/test_fft_filter_080µm.sur", & !
tab_s = heights_out(1:nx, 1:ny), & !
scal = scal_surf ) !
write(*,*) 'gaussian filter cutoff = 80 µm'
!========================================================================= fft_filter 150 µm
fft_cutoff = dx / 150.e-6
call fft_filter(tab = heights(1:nx, 1:ny), & !
long = nx, & !
larg = ny, & !
cutoff = fft_cutoff, & !
bf_tab = heights_out(1:nx, 1:ny), & !
multi_fft = .FALSE.) !
call write_surf( nom_fic = "out/test_fft_filter_150µm.sur", & !
tab_s = heights_out(1:nx, 1:ny), & !
scal = scal_surf ) !
write(*,*) 'gaussian filter cutoff = 150 µm'
call end_fftw3()
!========================================================================= ROLL +, RAY = 10 µm
call morpho_filter( mtype = "dilation", & ! IN
tabin = heights(1:nx, 1:ny), & ! IN
tabou = heights_out(1:nx, 1:ny), & ! OUT
long = nx, & ! IN
larg = ny, & ! IN
scale_xyz = [dx, dy, dz], & ! IN
ray = 10.e-6_R8, & ! IN
omp = .true., & ! IN
nb_div = 50 ) ! IN
call write_surf( nom_fic = "out/test_roll_smooth_dilation_10µm.sur", & !
tab_s = heights_out(1:nx, 1:ny), & !
scal = scal_surf ) !
write(*,*) 'morpho filter: dilation disk 10 µm'
!------------------------------------------------------------------------- ROLL -, RAY = 10 µm
call morpho_filter( mtype = "erosion", & ! IN
tabin = heights(1:nx, 1:ny), & ! IN
tabou = heights_out(1:nx, 1:ny), & ! OUT
long = nx, & ! IN
larg = ny, & ! IN
scale_xyz = [dx, dy, dz], & ! IN
ray = 10.e-6_R8, & ! IN
omp = .true., & ! IN
nb_div = 50 ) ! IN
call write_surf( nom_fic = "out/test_roll_smooth_erosion_10µm.sur", & !
tab_s = heights_out(1:nx, 1:ny), & !
scal = scal_surf ) !
write(*,*) 'morpho filter: erosion disk 10 µm'
!------------------------------------------------------------------------- CLOSING, RAY = 10 µm
call morpho_filter( mtype = "closing", & ! IN
tabin = heights(1:nx, 1:ny), & ! IN
tabou = heights_out(1:nx, 1:ny), & ! OUT
long = nx, & ! IN
larg = ny, & ! IN
scale_xyz = [dx, dy, dz], & ! IN
ray = 10.e-6_R8, & ! IN
omp = .true., & ! IN
nb_div = 50 ) ! IN
call write_surf( nom_fic = "out/test_roll_smooth_closing_10µm.sur", & !
tab_s = heights_out(1:nx, 1:ny), & !
scal = scal_surf ) !
write(*,*) 'morpho filter: closing disk 10 µm'
!------------------------------------------------------------------------- OPENING, RAY = 10 µm
call morpho_filter( mtype = "opening", & ! IN
tabin = heights(1:nx, 1:ny), & ! IN
tabou = heights_out(1:nx, 1:ny), & ! OUT
long = nx, & ! IN
larg = ny, & ! IN
scale_xyz = [dx, dy, dz], & ! IN
ray = 10.e-6_R8, & ! IN
omp = .true., & ! IN
nb_div = 50 ) ! IN
call write_surf( nom_fic = "out/test_roll_smooth_opening_10µm.sur", & !
tab_s = heights_out(1:nx, 1:ny), & !
scal = scal_surf ) !
write(*,*) 'morpho filter: opening disk 10 µm'
!========================================================================= KERNEL SIZE = 3
allocate( heights_copy(1:nx, 1:ny) )
heights_copy(1:nx, 1:ny) = heights(1:nx, 1:ny)
call median_smooth( tab = heights(1:nx, 1:ny), & ! INOUT
long = nx, & ! IN
larg = ny, & ! IN
kernel = 3, & ! IN
omp = .true. ) ! IN
call write_surf( nom_fic = "out/test_med_smooth_3.sur", & !
tab_s = heights(1:nx, 1:ny), & !
scal = scal_surf ) !
write(*,*) 'median filter kernel: 3'
!========================================================================= KERNEL SIZE = 5
heights(1:nx, 1:ny) = heights_copy(1:nx, 1:ny)
call median_smooth( tab = heights(1:nx, 1:ny), & ! INOUT
long = nx, & ! IN
larg = ny, & ! IN
kernel = 5, & ! IN
omp = .true. ) ! IN
call write_surf( nom_fic = "out/test_med_smooth_5.sur", & !
tab_s = heights(1:nx, 1:ny), & !
scal = scal_surf ) !
write(*,*) 'median filter kernel: 5'
!========================================================================= KERNEL SIZE = 9
heights(1:nx, 1:ny) = heights_copy(1:nx, 1:ny)
call median_smooth( tab = heights(1:nx, 1:ny), & ! INOUT
long = nx, & ! IN
larg = ny, & ! IN
kernel = 9, & ! IN
omp = .true. ) ! IN
call write_surf( nom_fic = "out/test_med_smooth_9.sur", & !
tab_s = heights(1:nx, 1:ny), & !
scal = scal_surf ) !
write(*,*) 'median filter kernel: 9'
!========================================================================= SOFTEN NO MASK
heights(1:nx, 1:ny) = heights_copy(1:nx, 1:ny)
call soften( tabin = heights(1:nx, 1:ny), & ! IN
tabout = heights_out(1:nx, 1:ny), & ! OUT
long = nx, & ! IN
larg = ny ) ! IN
call write_surf( nom_fic = "out/test_soften_no_mask.sur", & !
tab_s = heights_out(1:nx, 1:ny), & !
scal = scal_surf ) !
write(*,*) 'simple smooth with no mask'
!========================================================================= SOFTEN WITH MASK
call read_surf( nom_fic = "sur/test_mask.sur", & ! IN
tab_s = mask, & ! OUT
scal = scal_mask ) ! OUT
call soften( tabin = heights(1:nx, 1:ny), & ! IN
tabout = heights_out(1:nx, 1:ny), & ! OUT
mask = nint(mask(1:nx, 1:ny)), & ! IN
long = nx, & ! IN
larg = ny ) ! IN
call write_surf( nom_fic = "out/test_soften_with_mask.sur", & !
tab_s = heights_out(1:nx, 1:ny), & !
scal = scal_surf ) !
write(*,*) 'simple smooth in upper mask'
!========================================================================= MEDIAN FILTER SNB = 10, KERN = 15
call median_filter( tab = heights(1:nx, 1:ny), & ! INOUT
long = nx, & ! IN
larg = ny, & ! IN
kernel = 15, & ! IN
snb = 10, & ! IN
sig = 3._R8, & ! IN
omp = .true. ) ! IN
call write_surf( nom_fic = "out/test_med_filt_k15_s10.sur", & !
tab_s = heights(1:nx, 1:ny), & !
scal = scal_surf ) !
deallocate( heights, heights_copy, heights_out, mask )
stop
endprogram test_smooth