From 95c8ac3e464cf4cef5e1b7a2c50a5184b61a2ed0 Mon Sep 17 00:00:00 2001 From: Danny Hermes Date: Mon, 20 Nov 2017 23:02:19 -0800 Subject: [PATCH] Adding `interior_combine()` to Fortran implementation. --- src/bezier/surface_intersection.f90 | 128 +++++++++++++++++++- tests/fortran/test_surface_intersection.f90 | 126 ++++++++++++++++++- 2 files changed, 247 insertions(+), 7 deletions(-) diff --git a/src/bezier/surface_intersection.f90 b/src/bezier/surface_intersection.f90 index 71507ae4..b61fbf84 100644 --- a/src/bezier/surface_intersection.f90 +++ b/src/bezier/surface_intersection.f90 @@ -25,11 +25,11 @@ module surface_intersection evaluate_barycentric, jacobian_both, subdivide_nodes, compute_edge_nodes implicit none private & - LocateCandidate, MAX_LOCATE_SUBDIVISIONS, LOCATE_EPS, & + LocateCandidate, MAX_LOCATE_SUBDIVISIONS, LOCATE_EPS, MAX_EDGES, & newton_refine_solve, split_candidate, allocate_candidates, & update_candidates, ignored_edge_corner, ignored_double_corner, & ignored_corner, classify_tangent_intersection, no_intersections, & - remove_node + remove_node, finalize_segment public & Intersection, CurvedPolygonSegment, & IntersectionClassification_FIRST, IntersectionClassification_SECOND, & @@ -40,7 +40,7 @@ module surface_intersection SurfaceContained_FIRST, SurfaceContained_SECOND, newton_refine, & locate_point, classify_intersection, add_st_vals, & surfaces_intersection_points, get_next, to_front, add_segment, & - surfaces_intersect + interior_combine, surfaces_intersect ! NOTE: This (for now) is not meant to be C-interoperable. type :: Intersection @@ -68,6 +68,7 @@ module surface_intersection ! NOTE: These values are also defined in equivalent Python source. integer(c_int), parameter :: MAX_LOCATE_SUBDIVISIONS = 20 real(c_double), parameter :: LOCATE_EPS = 0.5_dp**47 + integer(c_int), parameter :: MAX_EDGES = 10 ! Values of IntersectionClassification enum: integer(c_int), parameter :: IntersectionClassification_FIRST = 0 integer(c_int), parameter :: IntersectionClassification_SECOND = 1 @@ -1140,6 +1141,127 @@ subroutine add_segment(curr_node, next_node, index, segments) end subroutine add_segment + subroutine finalize_segment(segment_ends, index, count) + + integer(c_int), allocatable, intent(inout) :: segment_ends(:) + integer(c_int), intent(inout) :: index + integer(c_int), intent(in) :: count + ! Variables outside of signature. + integer(c_int) :: curr_size + integer(c_int), allocatable :: segment_ends_swap(:) + + if (allocated(segment_ends)) then + curr_size = size(segment_ends) + if (curr_size < index) then + allocate(segment_ends_swap(index)) + segment_ends_swap(:curr_size) = segment_ends(:curr_size) + call move_alloc(segment_ends_swap, segment_ends) + end if + else + allocate(segment_ends(index)) + end if + + segment_ends(index) = count + index = index + 1 + + end subroutine finalize_segment + + subroutine interior_combine( & + num_intersections, intersections, & + segment_ends, segments, status) + + ! NOTE: This subroutine is not meant to be part of the interface for this + ! module, but it is (for now) public, so that it can be tested. + ! NOTE: This assumes, but does not check, that ``num_intersections > 0``. + ! NOTE: If ``segments`` is already allocated / with data, the data will be + ! overwritten in place. + + ! Possible error states: + ! * Status_SUCCESS: On success. + ! * Status_UNKNOWN: If a curved polygon requires more than ``MAX_EDGES`` + ! sides. (This could be due to either a particular + ! complex intersection or a programming error which + ! causes ``at_start`` to never be true.) + + integer(c_int), intent(in) :: num_intersections + type(Intersection), intent(in) :: intersections(num_intersections) + integer(c_int), allocatable, intent(out) :: segment_ends(:) + type(CurvedPolygonSegment), allocatable, intent(inout) :: segments(:) + integer(c_int), intent(out) :: status + ! Variables outside of signature. + integer(c_int) :: unused(num_intersections) + integer(c_int) :: segment_index, segment_ends_index + integer(c_int) :: remaining, i, start + type(Intersection) :: curr_node, next_node + logical(c_bool) :: at_start + + ! In this case, we know we'll have at least ``num_intersections`` segments + ! (and may have more due to corners). + if (allocated(segments)) then + if (size(segments) < num_intersections) then + deallocate(segments) + allocate(segments(num_intersections)) + end if + else + allocate(segments(num_intersections)) + end if + + status = Status_SUCCESS + + ! Set all of the unused indices. + unused = [ (i, i = 1, num_intersections) ] + remaining = num_intersections + segment_index = 1 + segment_ends_index = 1 + do while (remaining > 0) + ! "Pop" off intersection from the end to start a curved polygon. + start = unused(remaining) + remaining = remaining - 1 + + curr_node = intersections(start) + at_start = .FALSE. + ! Now move around the edge segments until the curved polygon's + ! entire boundary has been traversed. + edge_loop: do i = 1, MAX_EDGES + ! Follow along the current edge to the next node. + call get_next( & + num_intersections, intersections, unused, remaining, & + start, curr_node, next_node, at_start) + call add_segment(curr_node, next_node, segment_index, segments) + ! After adding the segment, check if we are back where we started. + if (at_start) then + exit edge_loop + end if + + ! Now the ``next_node`` becomes the current node, but we may + ! need to "rotate" it to the next edge if on a corner. + call to_front( & + num_intersections, intersections, unused, remaining, & + start, next_node, curr_node, at_start) + ! After "rotating" to the next edge, check if we are back where + ! we started. + if (at_start) then + exit edge_loop + end if + end do edge_loop + + if (at_start) then + ! NOTE: We pass ``segment_index - 1`` because it has been incremented + ! one past the last written index. + call finalize_segment( & + segment_ends, segment_ends_index, segment_index - 1) + else + ! If the loop terminated without reaching the start node, then + ! we have encountered an error. + ! LCOV_EXCL_START + status = Status_UNKNOWN + return + ! LCOV_EXCL_STOP + end if + end do + + end subroutine interior_combine + subroutine surfaces_intersect( & num_nodes1, nodes1, degree1, & num_nodes2, nodes2, degree2, & diff --git a/tests/fortran/test_surface_intersection.f90 b/tests/fortran/test_surface_intersection.f90 index 05e1cab1..7e65573f 100644 --- a/tests/fortran/test_surface_intersection.f90 +++ b/tests/fortran/test_surface_intersection.f90 @@ -27,7 +27,7 @@ module test_surface_intersection SurfaceContained_FIRST, SurfaceContained_SECOND, newton_refine, & locate_point, classify_intersection, add_st_vals, & surfaces_intersection_points, get_next, to_front, add_segment, & - surfaces_intersect + interior_combine, surfaces_intersect use types, only: dp use unit_test_helpers, only: print_status implicit none @@ -35,7 +35,8 @@ module test_surface_intersection test_newton_refine, test_locate_point, test_classify_intersection, & test_add_st_vals, test_surfaces_intersection_points, & intersection_check, intersection_equal, test_get_next, test_to_front, & - test_add_segment, test_surfaces_intersect + test_add_segment, test_interior_combine, segment_check, & + test_surfaces_intersect public surface_intersection_all_tests contains @@ -51,6 +52,7 @@ subroutine surface_intersection_all_tests(success) call test_get_next(success) call test_to_front(success) call test_add_segment(success) + call test_interior_combine(success) call test_surfaces_intersect(success) end subroutine surface_intersection_all_tests @@ -839,7 +841,7 @@ function intersection_check( & intersection_, s, t, index_first, index_second, & interior_curve) result(same) - type(intersection), intent(in) :: intersection_ + type(Intersection), intent(in) :: intersection_ real(c_double), intent(in) :: s, t integer(c_int), intent(in) :: index_first, index_second integer(c_int), intent(in) :: interior_curve @@ -856,7 +858,7 @@ end function intersection_check function intersection_equal(intersection1, intersection2) result(same) - type(intersection), intent(in) :: intersection1, intersection2 + type(Intersection), intent(in) :: intersection1, intersection2 logical(c_bool) :: same same = intersection_check( & @@ -1179,6 +1181,122 @@ subroutine test_add_segment(success) end subroutine test_add_segment + subroutine test_interior_combine(success) + logical(c_bool), intent(inout) :: success + ! Variables outside of signature. + logical :: case_success + integer :: case_id + character(16) :: name + integer(c_int) :: first, second + type(Intersection) :: intersections(4) + integer(c_int), allocatable :: segment_ends(:) + type(CurvedPolygonSegment), allocatable :: segments(:) + integer(c_int) :: status + + case_id = 1 + name = "interior_combine" + + first = IntersectionClassification_FIRST + second = IntersectionClassification_SECOND + + ! CASE 1: ``segments`` and ``segment_ends`` are not allocated, + ! from 1L-6L (ID: 23). + intersections(1) = Intersection(0.25_dp, 0.75_dp, 1, 3, second) + intersections(2) = Intersection(0.75_dp, 0.25_dp, 3, 1, first) + + case_success = ( & + .NOT. allocated(segments) .AND. & + .NOT. allocated(segment_ends)) + call interior_combine( & + 2, intersections(:2), segment_ends, segments, status) + case_success = ( & + case_success .AND. & + allocated(segment_ends) .AND. & + all(segment_ends == [4]) .AND. & + allocated(segments) .AND. & + size(segments) == 4 .AND. & + segment_check(segments(1), 0.75_dp, 1.0_dp, 3) .AND. & + segment_check(segments(2), 0.0_dp, 0.25_dp, 1) .AND. & + segment_check(segments(3), 0.75_dp, 1.0_dp, 6) .AND. & + segment_check(segments(4), 0.0_dp, 0.25_dp, 4) .AND. & + status == Status_SUCCESS) + call print_status(name, case_id, case_success, success) + + ! CASE 2: ``segments`` and ``segment_ends`` must be re-allocated; + ! two disjoint intersections, from 13L-38Q (ID: 39). + intersections(1) = Intersection(0.625_dp, 0.25_dp, 1, 1, first) + intersections(2) = Intersection(0.375_dp, 0.75_dp, 1, 1, second) + intersections(3) = Intersection(0.3125_dp, 0.25_dp, 1, 2, first) + intersections(4) = Intersection(0.6875_dp, 0.75_dp, 1, 3, second) + + case_success = ( & + allocated(segment_ends) .AND. & + size(segment_ends) == 1 .AND. & + allocated(segments) .AND. & + size(segments) == 4) + call interior_combine( & + 4, intersections, segment_ends, segments, status) + case_success = ( & + case_success .AND. & + allocated(segment_ends) .AND. & + all(segment_ends == [3, 6]) .AND. & + allocated(segments) .AND. & + size(segments) == 6 .AND. & + segment_check(segments(1), 0.75_dp, 1.0_dp, 6) .AND. & + segment_check(segments(2), 0.0_dp, 0.25_dp, 4) .AND. & + segment_check(segments(3), 0.625_dp, 0.6875_dp, 1) .AND. & + segment_check(segments(4), 0.3125_dp, 0.375_dp, 1) .AND. & + segment_check(segments(5), 0.75_dp, 1.0_dp, 4) .AND. & + segment_check(segments(6), 0.0_dp, 0.25_dp, 5) .AND. & + status == Status_SUCCESS) + call print_status(name, case_id, case_success, success) + + ! CASE 3: ``segments`` has been allocated but is not large enough; + ! intersection terminates after node is "rotated" via + ! ``to_front()``, from 1L-28Q (ID: 25) with rounded parameters + ! from ``s = 0.40587727318528016`` and + ! ``t == 0.59412272681471978``. + intersections(1) = Intersection(0.375_dp, 0.5625_dp, 2, 2, second) + intersections(2) = Intersection(0.1875_dp, 0.0_dp, 1, 3, first) + ! Note that the corner is the last intersection, which will be used + ! as ``start``. + deallocate(segments) + allocate(segments(1)) ! Too small. + case_success = ( & + allocated(segment_ends) .AND. & + size(segment_ends) == 2 .AND. & + allocated(segments) .AND. & + size(segments) == 1) + call interior_combine( & + 2, intersections(:2), segment_ends, segments, status) + case_success = ( & + case_success .AND. & + allocated(segment_ends) .AND. & + all(segment_ends == [3]) .AND. & + allocated(segments) .AND. & + size(segments) == 3 .AND. & + segment_check(segments(1), 0.1875_dp, 1.0_dp, 1) .AND. & + segment_check(segments(2), 0.0_dp, 0.375_dp, 2) .AND. & + segment_check(segments(3), 0.5625_dp, 1.0_dp, 5) .AND. & + status == Status_SUCCESS) + call print_status(name, case_id, case_success, success) + + end subroutine test_interior_combine + + function segment_check(segment, start, end_, edge_index) result(same) + + type(CurvedPolygonSegment), intent(in) :: segment + real(c_double), intent(in) :: start, end_ + integer(c_int), intent(in) :: edge_index + logical(c_bool) :: same + + same = ( & + segment%start == start .AND. & + segment%end_ == end_ .AND. & + segment%edge_index == edge_index) + + end function segment_check + subroutine test_surfaces_intersect(success) logical(c_bool), intent(inout) :: success ! Variables outside of signature.