label_surf_summits Subroutine

public subroutine label_surf_summits(tab, nx, ny, valleys, peaks, saddles, nb_summits)

Function to output the extrema of a 2D array, as peaks, valleys or saddles.

Arguments

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

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

number of pixels along x
integer(kind=I4), intent(in) :: ny

number of pixels along y
integer(kind=I4), intent(out), dimension(:,:), allocatable :: valleys

list of valley coordinates
integer(kind=I4), intent(out), dimension(:,:), allocatable :: peaks

list of peaks coordinates
integer(kind=I4), intent(out), dimension(:,:), allocatable :: saddles

list of saddles coordinates
integer(kind=I4), intent(out), dimension(1:3) :: nb_summits

number of extrema of each kind

Calls

proc~~label_surf_summits~~CallsGraph proc~label_surf_summits label_surf_summits proc~gradient_corner gradient_corner proc~label_surf_summits->proc~gradient_corner proc~labelize_point labelize_point proc~label_surf_summits->proc~labelize_point selectcase selectcase proc~label_surf_summits->selectcase proc~deriv_n deriv_N proc~labelize_point->proc~deriv_n

Called by

proc~~label_surf_summits~~CalledByGraph proc~label_surf_summits label_surf_summits proc~peaks_and_pits_curvatures peaks_and_pits_curvatures proc~peaks_and_pits_curvatures->proc~label_surf_summits proc~test_label_surf_summits test_label_surf_summits proc~test_label_surf_summits->proc~label_surf_summits proc~test_peaks_and_pits_curvatures test_peaks_and_pits_curvatures proc~test_peaks_and_pits_curvatures->proc~peaks_and_pits_curvatures program~test_grad_curv test_grad_curv program~test_grad_curv->proc~test_label_surf_summits program~test_grad_curv->proc~test_peaks_and_pits_curvatures

Source Code

   subroutine label_surf_summits(tab, nx, ny, valleys, peaks, saddles, nb_summits)
   !================================================================================================
   !! Function to output the extrema of a 2D array, as peaks, valleys or saddles.
   implicit none
   integer(kind=I4), intent(in ) :: nx                                     !! *number of pixels along x*
   integer(kind=I4), intent(in ) :: ny                                     !! *number of pixels along y*
   real   (kind=R8), intent(in ), dimension(1:nx, 1:ny) :: tab             !! *input 2D array*
   integer(kind=I4), intent(out), dimension(1:3)        :: nb_summits      !! *number of extrema of each kind*
   integer(kind=I4), intent(out), dimension(:,:), allocatable :: valleys   !! *list of valley coordinates*
   integer(kind=I4), intent(out), dimension(:,:), allocatable :: peaks     !! *list of peaks coordinates*
   integer(kind=I4), intent(out), dimension(:,:), allocatable :: saddles   !! *list of saddles coordinates*

      integer(kind=I4) :: i, j
      integer(kind=I4) :: ip ! peak counter
      integer(kind=I4) :: iv ! valley counter
      integer(kind=I4) :: is ! saddle counter

      character(len=1) :: label

      integer(kind=I4), dimension(1:3, 1:3) :: tp

      integer(kind=I4), dimension(:,:), allocatable :: topo ! point kind 2D array

      real(kind=R8), dimension(1:3, 1:3) :: ht

      real(kind=R8), dimension(1:3) :: gdx, gdy ! gradient vectors

      allocate( topo(1:nx, 1:ny) )

      topo(1:nx, 1:ny) = 0

      ! Loop through each point in the surface
      ip = 1 ; iv = 1 ; is = 1
      do i = 1 + 1, nx - 1

         do j = 1 + 1, ny - 1

            tp(1:3, 1:3) = topo(i-1:i+1, j-1:j+1)
            ht(1:3, 1:3) =  tab(i-1:i+1, j-1:j+1)

            ! if the gradients along x (resp. y) have the same sign, there is no extremum in the middle node
            call gradient_corner( hgt = ht(1:3, 1:3),   &  ! in
                                  gdx = gdx(1:3),       &  ! out
                                  gdy = gdy(1:3) )         ! out

!~             if ( all( [gdx(6) > 0, gdx(8) > 0, gdx(9) > 0] ) .or. all( [gdy(5) > 0, gdy(7) > 0, gdy(9) > 0] ) ) cycle
!~             if ( all( [gdx(6) > 0, gdx(8) > 0, gdx(9) > 0] ) .or. all( [gdy(5) < 0, gdy(7) < 0, gdy(9) < 0] ) ) cycle
!~             if ( all( [gdx(6) < 0, gdx(8) < 0, gdx(9) < 0] ) .or. all( [gdy(5) < 0, gdy(7) < 0, gdy(9) < 0] ) ) cycle
!~             if ( all( [gdx(6) < 0, gdx(8) < 0, gdx(9) < 0] ) .or. all( [gdy(5) > 0, gdy(7) > 0, gdy(9) > 0] ) ) cycle

            if ( all( gdx(1:3) > 0 ) .or. all( gdy(1:3) > 0 ) ) cycle
            if ( all( gdx(1:3) > 0 ) .or. all( gdy(1:3) < 0 ) ) cycle
            if ( all( gdx(1:3) < 0 ) .or. all( gdy(1:3) < 0 ) ) cycle
            if ( all( gdx(1:3) < 0 ) .or. all( gdy(1:3) > 0 ) ) cycle

            ! condition to avoid summits glued to each other: if a summit has been found in the neighborhood, cycle.
            if ( any( [tp(1, 1) > 0, &  !
                       tp(1, 2) > 0, &  !
                       tp(1, 3) > 0, &  !
                       tp(2, 1) > 0] ) ) cycle

            call labelize_point( height = ht(1:3, 1:3), &  ! in
                                 label  = label )          ! out

            selectcase(label)
               case('V') ; topo(i, j) = 1 ; iv = iv + 1 ! one valley more detected
               case('S') ; topo(i, j) = 2 ; is = is + 1 ! one saddle more detected
               case('P') ; topo(i, j) = 3 ; ip = ip + 1 ! one peak   more detected
            endselect

         enddo

      enddo

      nb_summits = [iv - 1, is - 1, ip - 1]

      ! now the number of extrema is known
      allocate( valleys(1:nb_summits(1), 1:2) )
      allocate( saddles(1:nb_summits(2), 1:2) )
      allocate( peaks  (1:nb_summits(3), 1:2) )

      ip = 1 ; iv = 1 ; is = 1
      do i = 1 + 1, nx - 1

         do j = 1 + 1, ny - 1

            selectcase( topo(i, j) )
               case(1) ; valleys(iv, 1) = i ;  valleys(iv, 2) = j ; iv = iv + 1
               case(2) ; saddles(is, 1) = i ;  saddles(is, 2) = j ; is = is + 1
               case(3) ; peaks  (ip, 1) = i ;  peaks  (ip, 2) = j ; ip = ip + 1
            endselect

         enddo

      enddo

      deallocate( topo )

   return
   endsubroutine label_surf_summits