kdtree.f90

kd tree java implementation 2

Module kd_tree

  ! K-D tree routines in Fortran 90 by Matt Kennel.
  ! Original program was written in Sather by Steve Omohundro and
  ! Matt Kennel.  Only the Euclidean metric works so far.


  ! Institute for Nonlinear science, UC San Diego, 1997
  ! [e-mail]
  ! .. Parameters ..
  ! you choose this.
  Integer, Parameter :: bucket_size = 10
  ! ..
  ! .. Derived Type Declarations ..
  ! Global information about the tree
  ! pointer to the actual
  ! data array
  ! dimensionality and total # of points
  ! permuted index into the
  ! data, so that
  ! indexes[l..u]
  ! of some bucket represent the indexes of the actual
  ! points in that bucket.
  ! root pointer of the tree
  ! an internal tree node
  ! the dimension to cut
  ! where to cut the dimension
  ! indices of points included in this node,
  ! referring back to indices
  ! child pointers
  ! One of these is created for each search.
  ! best squared distance found so far
  ! best index found so far
  ! Query vector
  ! indexes of best found so far
  ! squared distances found so far
  ! how many best distances we are searching
  ! for
  ! how many have been found so far, i.e. il(1:nfound)
  ! are the only valid indexes.
  ! exclude points within
  ! 'correltime'
  ! of 'centeridx'

  Type :: tree_node
      Integer :: dnum
      Real :: val
      Integer :: l, u
      Type (tree_node), Pointer :: left, right
  End Type tree_node

  Type :: tree_master_record
      Real, Dimension (:,:), Pointer :: the_data
      Integer :: dim, n
      Integer, Dimension (:), Pointer :: indexes
      Type (tree_node), Pointer :: root
  End Type tree_master_record

  Type :: tree_search_record
      Private
      Real :: bsd
      Real, Dimension (:), Pointer :: qv
      Integer, Dimension (:), Pointer :: il
      Real, Dimension (:), Pointer :: dsl
      Integer :: nbst, nfound
      Integer :: centeridx, correltime
      Logical :: linfinity_metric
  End Type tree_search_record
  ! ..
  ! set to true if we're doing linfinity metric
Contains

  Subroutine destroy_tree(tp)
    ! Deallocates all memory for the tree, except input data matrix
    ! .. Structure Arguments ..
    Type (tree_master_record), Pointer :: tp
    ! ..
    Call destroy_node(tp%root)
    Deallocate (tp%indexes)
    tp%indexes => NULL() !#Nullify (tp%indexes)
    Return
  end Subroutine destroy_tree
  Recursive Subroutine destroy_node(np)
    ! .. Structure Arguments ..
    Type (tree_node), Pointer :: np
    ! ..
    ! ..
    If (ASSOCIATED(np%left)) Then
       Call destroy_node(np%left)
       Deallocate (np%left)
       np%left => NULL() 
    End If
    If (ASSOCIATED(np%right)) Then
       Call destroy_node(np%right)
       Deallocate (np%right)
       np%right => NULL() 
    End If
    Return
  End Subroutine destroy_node


  Function create_tree(input_data) Result (master_record)
    ! create the actual tree structure, given an input array of data.
    ! Arguments
    ! .. Function Return Value ..
    Type (tree_master_record), Pointer :: master_record
    ! ..
    ! .. Array Arguments ..
    Real, Target :: input_data(:,:)
    ! ..
    ! ..
    Allocate (master_record)
    master_record%the_data => input_data
    ! pointer assignment
    master_record%n = SIZE(input_data,1)
    master_record%dim = SIZE(input_data,2)
    Print *, 'Creating KD tree with N = ', master_record%n, &
     ' and dim = ', master_record%dim
    Call build_tree(master_record)

  Contains

    Subroutine build_tree(tp)
      ! .. Structure Arguments ..
      Type (tree_master_record), Pointer :: tp
      ! ..
      ! .. Local Scalars ..
      Integer :: j
      ! ..
      Allocate (tp%indexes(tp%n))
      Do j = 1, tp%n
         tp%indexes(j) = j
      End Do
      tp%root => build_tree_for_range(tp,1,tp%n)
    End Subroutine build_tree

    Recursive Function build_tree_for_range(tp,l,u) Result (res)
      ! .. Function Return Value ..
      Type (tree_node), Pointer :: res
      ! ..
      ! .. Structure Arguments ..
      Type (tree_master_record), Pointer :: tp
      ! ..
      ! .. Scalar Arguments ..
      Integer, Intent (In) :: l, u
      ! ..
      ! .. Local Scalars ..
      Integer :: c, m
      ! ..
      If ((u-l)<=bucket_size) Then
         Allocate (res)
         res%dnum = 0
         res%val = 0.0
         res%l = l
         res%u = u
         Nullify (res%left,res%right)
      Else
         Allocate (res)
         c = most_spread_coordinate(tp,l,u)
         m = (l+u)/2
         Call select_on_coordinate(tp,c,m,l,u)
         ! moves indexes around
         res%dnum = c
         res%val = tp%the_data(tp%indexes(m),c)
         res%l = l
         res%u = u
         res%left => build_tree_for_range(tp,l,m)
         res%right => build_tree_for_range(tp,m+1,u)
      End If
    End Function build_tree_for_range

    Subroutine select_on_coordinate(tp,c,k,li,ui)
      ! Move elts of ind around between l and u, so that the kth
      ! element
      ! is >= those below, <= those above, in the coordinate c.
      ! .. Structure Arguments ..
      Type (tree_master_record), Pointer :: tp
      ! ..
      ! .. Scalar Arguments ..
      Integer, Intent (In) :: c, k, li, ui
      ! ..
      ! .. Local Scalars ..
      Integer :: i, l, m, s, t, u
      ! ..
      ! .. Local Arrays ..
      Real, Pointer :: v(:,:)
      Integer, Pointer :: ind(:)
      ! ..
      v => tp%the_data
      ind => tp%indexes
      l = li
      u = ui
      Do While (l<u)
         t = ind(l)
         m = l
         Do i = l + 1, u
            If (v(ind(i),c)<v(t,c)) Then
               m = m + 1
               s = ind(m)
               ind(m) = ind(i)
               ind(i) = s
            End If
         End Do
         s = ind(l)
         ind(l) = ind(m)
         ind(m) = s
         If (m<=k) l = m + 1
         If (m>=k) u = m - 1
      End Do
    End Subroutine select_on_coordinate

    Function most_spread_coordinate(tp,l,u) Result (res)
      ! Of indices in l..u find the axis which has the largest spread,
      ! and
      ! return its index
      ! .. Function Return Value ..
      Integer :: res
      ! ..
      ! .. Structure Arguments ..
      Type (tree_master_record), Pointer :: tp
      ! ..
      ! .. Scalar Arguments ..
      Integer, Intent (In) :: l, u
      ! ..
      ! .. Local Scalars ..
      Real :: bsp, sp
      Integer :: i = 0
      ! ..
      res = 0
      bsp = -1.0
      Do i = 1, tp%dim
         sp = spread_in_coordinate(tp,i,l,u)
         If (sp>bsp) Then
            res = i
            bsp = sp
         End If
      End Do
    End Function most_spread_coordinate

    Function spread_in_coordinate(tp,c,l,u) Result (res)
      ! the spread in coordinate 'c', between l and u
      ! for easier local access
      ! ibid
      ! .. Function Return Value ..
      Real :: res
      ! ..
      ! .. Structure Arguments ..
      Type (tree_master_record), Pointer :: tp
      ! ..
      ! .. Scalar Arguments ..
      Integer, Intent (In) :: c, l, u
      ! ..
      ! .. Local Scalars ..
      Real :: last, lmax, lmin, max, min, t
      Integer :: i
      ! ..
      ! .. Local Arrays ..
      Real, Pointer :: v(:,:)
      Integer, Pointer :: ind(:)
      ! ..
      v => tp%the_data
      ind => tp%indexes
      min = v(ind(l),c)
      max = min
      Do i = l + 2, u, 2
         lmin = v(ind(i-1),c)
         lmax = v(ind(i),c)
         If (lmin>lmax) Then
            t = lmin
            lmin = lmax
            lmax = t
         End If
         If (min>lmin) min = lmin
         If (max<lmax) max = lmax
      End Do
      If (i==u+1) Then
         last = v(ind(u),c)
         If (min>last) min = last