89 integer :: object_index = -1
104 procedure, pass(this),
public :: get_object_index => &
106 procedure, pass(this),
public :: get_parent_index => &
108 procedure, pass(this),
public :: get_left_index => &
110 procedure, pass(this),
public :: get_right_index => &
121 generic ::
operator(.lt.) => less
122 generic ::
operator(.gt.) => greater
134 integer :: allocated_node_count = 0
136 integer :: node_capacity = 0
137 integer :: growth_size = 1
144 procedure, pass(this),
public :: insert_object => &
150 procedure, pass(this),
public :: get_root_index => &
152 procedure, pass(this),
public :: get_parent_index => &
154 procedure, pass(this),
public :: get_left_index => &
156 procedure, pass(this),
public :: get_right_index => &
160 procedure, pass(this),
public :: get_root_node => &
162 procedure, pass(this),
public :: get_parent_node => &
164 procedure, pass(this),
public :: get_left_node => &
166 procedure, pass(this),
public :: get_right_node => &
171 procedure, pass(this),
public :: query_overlaps => &
200 this%object_index = -1
221 integer :: object_index
223 object_index = this%object_index
229 integer :: parent_index
231 parent_index = this%parent_node_index
237 integer :: left_index
239 left_index = this%left_node_index
245 integer :: right_index
247 right_index = this%right_node_index
253 real(kind=
dp),
dimension(3),
intent(in) :: p
254 real(kind=
dp) :: distance
256 distance = 0.5_rp * this%aabb%get_diameter() &
257 - norm2(this%aabb%get_center() - p)
277 if (this%is_leaf())
then
281 & this%object_index .gt. 0
286 & this%object_index .eq. -1
300 res = this%aabb .lt. other%aabb
310 res = this%aabb .gt. other%aabb
321 integer,
intent(in) :: initial_capacity
326 this%allocated_node_count = 0
327 this%next_free_node_index = 1
328 this%node_capacity = initial_capacity
329 this%growth_size = initial_capacity
331 if (
allocated(this%nodes))
deallocate(this%nodes)
332 allocate(this%nodes(initial_capacity))
334 do i = 1, initial_capacity
335 this%nodes(i)%next_node_index = i + 1
343 class(*),
dimension(:),
intent(in) :: objects
345 integer :: i_obj, i_node, i
348 integer :: start_layer, end_layer
350 type(
aabb_t),
dimension(:),
allocatable :: box_list
351 integer,
dimension(:),
allocatable :: sorted_indices
353 call this%init(
size(objects) * 2)
359 allocate(box_list(
size(objects)))
361 do i_obj = 1,
size(objects)
362 box_list(i_obj) =
get_aabb(objects(i_obj))
364 call sort(box_list, sorted_indices)
366 do i = 1,
size(sorted_indices)
367 i_obj = sorted_indices(i)
368 i_node = this%allocate_node()
369 this%nodes(i_node)%aabb =
get_aabb(objects(i_obj))
370 this%nodes(i_node)%object_index = i_obj
375 end_layer =
size(objects)
377 do while (.not. done)
380 do i = start_layer, end_layer - 1, 2
381 i_node = this%allocate_node()
383 this%nodes(i_node)%aabb =
merge(this%nodes(i)%aabb, &
384 this%nodes(i + 1)%aabb)
386 this%nodes(i_node)%left_node_index = i
387 this%nodes(i_node)%right_node_index = i + 1
389 this%nodes(i)%parent_node_index = i_node
390 this%nodes(i + 1)%parent_node_index = i_node
395 if (mod(end_layer - start_layer, 2) .eq. 0)
then
396 i_node = this%allocate_node()
397 this%nodes(i_node)%aabb = this%nodes(end_layer)%aabb
398 this%nodes(i_node)%left_node_index = end_layer
401 this%nodes(end_layer)%parent_node_index = i_node
405 start_layer = end_layer + 1
406 end_layer = this%allocated_node_count
409 done = start_layer .eq. end_layer
413 this%root_node_index = this%allocated_node_count
415 if (this%get_size() .ne.
size(objects))
then
416 print *,
"this%get_size() = ", this%get_size()
417 print *,
"size(objects) = ",
size(objects)
424 subroutine sort(array, indices)
425 type(
aabb_t),
dimension(:),
intent(in) :: array
426 integer,
intent(inout),
dimension(:),
allocatable :: indices
427 logical,
dimension(:),
allocatable :: visited
432 allocate(indices(
size(array)))
433 allocate(visited(
size(array)))
437 do i = 1,
size(array)
439 do imin = 1,
size(array)
440 if (.not. visited(imin) .and. minidx .eq. -1) minidx = imin
441 if (minidx .gt. -1)
then
442 if (visited(imin) .and. array(imin) .lt. array(minidx)) minidx = imin
447 visited(minidx) = .true.
463 call simple_stack%init(this%allocated_node_count)
465 tmp = this%get_root_index()
467 call simple_stack%push(tmp)
470 do while (.not. simple_stack%is_empty())
471 idx = simple_stack%pop()
474 if (this%nodes(idx)%is_leaf())
then
477 tmp = this%get_left_index(idx)
478 call simple_stack%push(tmp)
479 tmp = this%get_right_index(idx)
480 call simple_stack%push(tmp)
492 integer :: root_index
494 root_index = this%root_node_index
501 integer,
intent(in) :: node_index
502 integer :: parent_index
504 parent_index = this%nodes(node_index)%parent_node_index
511 integer,
intent(in) :: node_index
512 integer :: left_index
514 left_index = this%nodes(node_index)%left_node_index
521 integer,
intent(in) :: node_index
522 integer :: right_index
524 right_index = this%nodes(node_index)%right_node_index
533 integer,
intent(in) :: node_index
536 node = this%nodes(node_index)
544 root_node = this%nodes(this%root_node_index)
551 integer,
intent(in) :: node_index
554 parent_node = this%nodes(this%nodes(node_index)%parent_node_index)
560 integer,
intent(in) :: node_index
563 left_node = this%nodes(this%nodes(node_index)%left_node_index)
570 integer,
intent(in) :: node_index
573 right_node = this%nodes(this%nodes(node_index)%right_node_index)
578 integer,
intent(in) :: node_index
581 out_box = this%nodes(node_index)%aabb
589 class(*),
intent(in) :: object
590 integer,
intent(in) :: object_index
592 integer :: node_index
594 node_index = this%allocate_node()
595 this%nodes(node_index)%aabb =
get_aabb(object)
596 this%nodes(node_index)%object_index = object_index
598 call this%insert_leaf(node_index)
604 class(*),
intent(in) :: object
605 integer,
intent(in) :: object_index
606 integer,
intent(inout),
allocatable :: overlaps(:)
609 type(
aabb_t) :: object_box
611 integer :: root_index, left_index, right_index
612 integer :: node_index
614 if (
allocated(overlaps))
deallocate(overlaps)
615 allocate(overlaps(0))
618 root_index = this%get_root_index()
620 call simple_stack%init()
621 call simple_stack%push(root_index)
623 do while (.not. simple_stack%is_empty())
624 node_index = simple_stack%pop()
628 if (this%nodes(node_index)%aabb%overlaps(object_box))
then
629 if (this%nodes(node_index)%is_leaf())
then
630 if (this%nodes(node_index)%object_index .ne. object_index)
then
631 overlaps = [this%nodes(node_index)%object_index, overlaps]
634 left_index = this%get_left_index(node_index)
635 call simple_stack%push(left_index)
636 right_index = this%get_right_index(node_index)
637 call simple_stack%push(right_index)
641 call simple_stack%free()
650 integer :: node_index
653 call this%resize_node_pool(this%node_capacity + this%growth_size)
656 node_index = this%next_free_node_index
658 associate(new_node => this%nodes(node_index))
659 this%next_free_node_index = new_node%next_node_index
665 this%next_free_node_index = new_node%next_node_index
666 this%allocated_node_count = this%allocated_node_count + 1
674 integer,
intent(in) :: node_index
676 this%nodes(node_index)%next_node_index = this%next_free_node_index
677 this%next_free_node_index = node_index
678 this%allocated_node_count = this%allocated_node_count - 1
684 integer,
intent(in) :: leaf_node_index
686 integer :: tree_node_index
688 real(kind=rp) :: cost_left
689 real(kind=rp) :: cost_right
696 type(aabb_t) :: combined_aabb
697 real(kind=rp) :: new_parent_node_cost
698 real(kind=rp) :: minimum_push_down_cost
699 type(aabb_t) :: new_left_aabb
700 type(aabb_t) :: new_right_aabb
701 integer :: leaf_sibling_index
703 integer :: old_parent_index
704 integer :: new_parent_index
709 leaf_node = this%nodes(leaf_node_index)
713 this%root_node_index = leaf_node_index
723 tree_node_index = this%root_node_index
724 tree_node = this%get_node(tree_node_index)
725 do while (.not. tree_node%is_leaf())
729 left_node = this%get_left_node(tree_node_index)
730 right_node = this%get_right_node(tree_node_index)
734 combined_aabb = merge(tree_node%aabb, leaf_node%get_aabb())
736 new_parent_node_cost = 2.0_rp * combined_aabb%get_surface_area()
737 minimum_push_down_cost = 2.0_rp * ( &
738 & combined_aabb%get_surface_area() &
739 & - tree_node%aabb%get_surface_area()&
744 if (left_node%is_leaf())
then
745 new_left_aabb = merge(leaf_node%aabb, left_node%get_aabb())
746 cost_left = new_left_aabb%get_surface_area() + minimum_push_down_cost
748 new_left_aabb = merge(leaf_node%aabb, left_node%get_aabb())
750 & new_left_aabb%get_surface_area() &
751 & - left_node%aabb%get_surface_area()&
752 & ) + minimum_push_down_cost
755 if (right_node%is_leaf())
then
756 new_right_aabb = merge(leaf_node%aabb, right_node%aabb)
757 cost_right = new_right_aabb%get_surface_area() + &
758 minimum_push_down_cost
760 new_right_aabb = merge(leaf_node%aabb, right_node%aabb)
762 & new_right_aabb%get_surface_area() &
763 & - right_node%aabb%get_surface_area() &
764 & ) + minimum_push_down_cost
770 if (new_parent_node_cost < cost_left .and. &
771 new_parent_node_cost < cost_right)
then
776 if (cost_left .lt. cost_right)
then
777 tree_node_index = tree_node%get_left_index()
779 tree_node_index = tree_node%get_right_index()
784 tree_node = this%get_node(tree_node_index)
789 leaf_sibling_index = tree_node_index
790 leaf_sibling = this%nodes(leaf_sibling_index)
791 old_parent_index = this%get_parent_index(leaf_sibling_index)
792 new_parent_index = this%allocate_node()
793 new_parent = this%nodes(new_parent_index)
794 new_parent%parent_node_index = old_parent_index
795 new_parent%aabb = merge(leaf_node%aabb, leaf_sibling%aabb)
797 if (leaf_node .lt. leaf_sibling)
then
798 new_parent%left_node_index = leaf_node_index
799 new_parent%right_node_index = leaf_sibling_index
801 new_parent%left_node_index = leaf_sibling_index
802 new_parent%right_node_index = leaf_node_index
805 leaf_node%parent_node_index = new_parent_index
806 leaf_sibling%parent_node_index = new_parent_index
810 this%root_node_index = new_parent_index
814 old_parent = this%nodes(old_parent_index)
815 if (old_parent%left_node_index .eq. leaf_sibling_index)
then
816 old_parent%left_node_index = new_parent_index
818 old_parent%right_node_index = new_parent_index
820 this%nodes(old_parent_index) = old_parent
823 this%nodes(leaf_node_index) = leaf_node
824 this%nodes(leaf_sibling_index) = leaf_sibling
825 this%nodes(new_parent_index) = new_parent
828 tree_node_index = leaf_node%parent_node_index
830 call this%fix_upwards_tree(tree_node_index)
839 type(stack_i4_t) :: simple_stack
840 integer :: current_index
841 integer :: root_index, left_index, right_index
848 root_index = this%get_root_index()
850 call simple_stack%init(this%node_capacity)
851 call simple_stack%push(root_index)
853 do while (.not. simple_stack%is_empty())
854 current_index = simple_stack%pop()
857 valid = valid .and. this%nodes(current_index)%is_valid()
859 if (.not. this%nodes(current_index)%is_leaf())
then
860 left_index = this%get_left_index(current_index)
861 right_index = this%get_right_index(current_index)
863 call simple_stack%push(left_index)
864 call simple_stack%push(right_index)
874 integer,
intent(in) :: tree_start_index
878 integer :: tree_node_index
880 tree_node_index = tree_start_index
882 left_node = this%get_left_node(tree_node_index)
883 right_node = this%get_right_node(tree_node_index)
885 this%nodes(tree_node_index)%aabb = merge(left_node%aabb, right_node%aabb)
887 tree_node_index = this%get_parent_index(tree_node_index)
894 type(stack_i4_t) :: simple_stack
896 integer :: current_index
897 integer :: root_index, left_index, right_index
899 root_index = this%get_root_index()
900 call simple_stack%init(this%node_capacity)
901 call simple_stack%push(root_index)
903 do while (.not. simple_stack%is_empty())
904 current_index = simple_stack%pop()
907 left_index = this%get_left_index(current_index)
908 right_index = this%get_right_index(current_index)
910 call simple_stack%push(this%nodes(current_index)%left_node_index)
911 call simple_stack%push(this%nodes(current_index)%right_node_index)
913 write(*, *)
"i = ", current_index
914 write(*, *)
" Parent : ", this%get_parent_index(current_index)
915 write(*, *)
" Children: ", this%get_left_index(current_index), &
916 this%get_right_index(current_index)
918 write(*, *)
" object_index = ", this%nodes(current_index)%object_index
926 integer,
intent(in) :: new_capacity
928 type(
aabb_node_t),
dimension(:),
allocatable :: temp
931 allocate(temp(new_capacity))
932 temp(:this%node_capacity) = this%nodes(:this%node_capacity)
934 do i = this%allocated_node_count, new_capacity
935 temp(i)%next_node_index = i + 1
939 call move_alloc(temp, this%nodes)
941 this%node_capacity = new_capacity
942 this%next_free_node_index = this%allocated_node_count + 1
Axis Aligned Bounding Box (aabb) Tree data structure.
pure type(aabb_t) function aabb_tree_get_aabb(this, node_index)
subroutine aabb_tree_fix_upwards_tree(this, tree_start_index)
Fixes the tree upwards.
pure logical function aabb_node_is_valid(this)
Returns true if the node is a valid node.
subroutine aabb_node_init(this)
Initializes the AABB node.
pure type(aabb_node_t) function aabb_tree_get_right_node(this, node_index)
Returns the right node of the node at the given index.
pure type(aabb_node_t) function aabb_tree_get_parent_node(this, node_index)
Returns the parent node of the node at the given index.
pure type(aabb_node_t) function aabb_tree_get_root_node(this)
Returns the root node of the tree.
subroutine aabb_tree_insert_object(this, object, object_index)
Inserts an object into the tree.
pure type(aabb_node_t) function aabb_tree_get_node(this, node_index)
Returns the node at the given index.
subroutine aabb_tree_print(this)
Prints the tree.
pure logical function aabb_node_greater(this, other)
Returns true if the node is greater than the other node.
integer function aabb_tree_allocate_node(this)
Allocates a new node in the tree.
pure type(aabb_t) function aabb_node_get_aabb(this)
Returns the Axis Aligned Bounding Box (aabb) of the node.
pure integer function aabb_tree_get_root_index(this)
Returns the index of the root node.
pure integer function aabb_node_get_right_index(this)
Returns the right index of the node.
subroutine aabb_tree_deallocate_node(this, node_index)
Deallocates a node in the tree.
pure integer function aabb_node_get_left_index(this)
Returns the left index of the node.
pure logical function aabb_node_is_leaf(this)
Returns true if the node is a leaf node.
real(kind=dp) function aabb_node_min_distance(this, p)
Get the minimum possible distance from the aabb to a point.
subroutine aabb_tree_init(this, initial_capacity)
Initializes the AABB tree.
subroutine aabb_tree_insert_leaf(this, leaf_node_index)
Inserts a leaf into the tree.
subroutine aabb_tree_query_overlaps(this, object, object_index, overlaps)
Queries the tree for overlapping objects.
subroutine aabb_tree_build_tree(this, objects)
Builds the tree.
subroutine sort(array, indices)
Return a list of sorted indices of the aabb nodes.
pure logical function aabb_node_less(this, other)
Returns true if the node is less than the other node.
pure integer function aabb_tree_get_parent_index(this, node_index)
Returns the index of the parent node of the node at the given index.
pure integer function aabb_node_get_parent_index(this)
Returns the parent index of the node.
pure integer function aabb_tree_get_left_index(this, node_index)
Returns the index of the left node of the node at the given index.
integer, parameter, public aabb_null_node
pure integer function aabb_node_get_object_index(this)
Returns the object index of the node.
logical function aabb_tree_valid_tree(this)
Validates the tree.
pure type(aabb_node_t) function aabb_tree_get_left_node(this, node_index)
Returns the left node of the node at the given index.
integer function aabb_tree_get_size(this)
Returns the size of the tree, in number of leaves.
subroutine aabb_tree_resize_node_pool(this, new_capacity)
Resizes the node pool.
pure integer function aabb_tree_get_right_index(this, node_index)
Returns the index of the right node of the node at the given index.
Axis Aligned Bounding Box (aabb) implementation in Fortran.
type(aabb_t) function, public get_aabb(object, padding)
Construct the aabb of a predefined object.
integer, parameter, public dp
integer, parameter, public rp
Global precision used in computations.
Implements a dynamic stack ADT.
Defines a triangular element.
Axis Aligned Bounding Box (aabb) data structure.
Node type for the Axis Aligned Bounding Box (aabb) Tree.
Axis Aligned Bounding Box (aabb) Tree.