subroutine bin_search_one_2d ( bin, nset, pset, nbin, bin_start, bin_next, & ptest, found_a_neighbor, i_min, dist_min_sq, compares ) !*****************************************************************************80 ! !! BIN_SEARCH_ONE_2D searches one cell in a 2D array of bins. ! ! Discussion: ! ! The bins are presumed to have been set up by successive calls to: ! ! R82VEC_BIN_EVEN2, ! R82VEC_BINNED_REORDER, and ! R82VEC_BINNED_SORT_A. ! ! Modified: ! ! 26 March 2001 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer BIN(2), the indices of the cell to be examined. ! ! Input, integer NSET, the number of points in the set. ! ! Input, real ( kind = 8 ) PSET(2,NSET), the coordinates of the points ! in the set. ! ! Input, integer NBIN(2), the number of cells in the horizontal and ! vertical directions. ! ! Input, integer BIN_START(NBIN(1),NBIN(2)), BIN_LAST(NBIN(1),NBIN(2)), ! indicates the index of the first and last element in the bin, or -1 ! if none. ! ! Input, integer BIN_NEXT(NSET), the index of the next element of the bin ! containing this element. ! ! Input, real ( kind = 8 ) PTEST(2), the coordinates of the test point. ! ! Input/output, logical FOUND_A_NEIGHBOR, is set to TRUE if at least ! one point of PSET is found in the current bin. Otherwise, it retains its ! input value. ! ! Input/output, integer I_MIN, the index of the nearest neighbor in ! PSET to PTEST, if at least one neighbor has been found. ! ! Input/output, real ( kind = 8 ) DIST_MIN_SQ, the square of the distance ! from the nearest neighbor in PSET to PTEST, if at least one neighbor ! has been found. ! ! Input/output, integer COMPARES, the number of elements of PSET whose ! distance to PTEST has been computed. ! implicit none integer, parameter :: ndim = 2 integer nbin(ndim) integer nset integer bin(ndim) integer bin_next(nset) integer bin_start(nbin(1),nbin(2)) integer compares real ( kind = 8 ) dist_min_sq real ( kind = 8 ) dist_sq logical found_a_neighbor integer i_min integer node real ( kind = 8 ) pset(ndim,nset) real ( kind = 8 ) ptest(ndim) node = bin_start(bin(1),bin(2)) do while ( 0 < node ) found_a_neighbor = .true. dist_sq = sum ( ( ptest(1:ndim) - pset(1:ndim,node) )**2 ) compares = compares + 1 if ( dist_sq < dist_min_sq ) then dist_min_sq = dist_sq i_min = node end if node = bin_next(node) end do return end subroutine bin_to_r8_even ( nbin, bin, a, b, cmin, cmax ) !*****************************************************************************80 ! !! BIN_TO_R8_EVEN returns the limits for a given "bin" in [A,B]. ! ! Discussion: ! ! The interval from A to B is divided into NBIN-2 equal subintervals or bins. ! An initial bin takes everything less than A, and a final bin takes ! everything greater than B. ! ! Example: ! ! NBIN = 7, A = 10, B = 20 ! ! BIN CMIN CMAX ! ! 1 -HUGE 10 ! 2 10 12 ! 3 12 14 ! 4 14 16 ! 5 16 18 ! 6 18 20 ! 7 20 HUGE ! ! Modified: ! ! 08 December 2004 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer NBIN, the number of bins. NBIN is normally at least 3. ! If NBIN is 1 or 2, then everything is assigned to bin 1. ! ! Input, integer BIN, the index of the bin to be considered. ! If BIN is less than 1, or greater than NBIN, the user will get what ! the user deserves. ! ! Input, real ( kind = 8 ) A, B, the lower and upper limits of the ! bin interval. While A is expected to be less than B, the code ! should return useful results if A is actually greater than B. ! ! Output, real ( kind = 8 ) CMIN, CMAX, the minimum and maximum ! limits on the bin. ! implicit none real ( kind = 8 ) a real ( kind = 8 ) b integer bin real ( kind = 8 ) cmax real ( kind = 8 ) cmin integer nbin ! ! Take care of special cases. ! if ( nbin <= 2 ) then cmin = -huge ( cmin ) cmax = huge ( cmax ) return end if if ( b == a ) then cmin = -huge ( cmin ) cmax = huge ( cmax ) return end if ! ! Compute the bin limits. ! if ( bin == 1 ) then cmin = -huge ( cmin ) cmax = a else if ( bin < nbin ) then cmin = ( real ( nbin - bin, kind = 8 ) * a & + real ( bin - 2, kind = 8 ) * b ) & / real ( nbin - 2, kind = 8 ) cmax = ( real ( nbin - bin - 1, kind = 8 ) * a & + real ( bin - 1, kind = 8 ) * b ) & / real ( nbin - 2, kind = 8 ) else if ( bin == nbin ) then cmin = b cmax = huge ( cmax ) else cmin = -huge ( cmin ) cmax = huge ( cmax ) end if return end subroutine bin_to_r8_even2 ( nbin, bin, a, b, cmin, cmax ) !*****************************************************************************80 ! !! BIN_TO_R8_EVEN2 returns the limits for a given "bin" in [A,B]. ! ! Discussion: ! ! The interval from A to B is divided into NBIN equal subintervals or bins. ! ! Example: ! ! NBIN = 5, A = 10, B = 20 ! ! BIN CMIN CMAX ! ! 1 10 12 ! 2 12 14 ! 3 14 16 ! 4 16 18 ! 5 18 20 ! ! Modified: ! ! 05 February 2001 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer NBIN, the number of bins. ! ! Input, integer BIN, the index of the bin to be considered. ! If BIN is less than 1, or greater than NBIN, the user will get what ! the user deserves. ! ! Input, real ( kind = 8 ) A, B, the lower and upper limits of the bin ! interval. While A is expected to be less than B, the code should ! return useful results if A is actually greater than B. ! ! Output, real ( kind = 8 ) CMIN, CMAX, the minimum and maximum limits ! on the bin. ! implicit none real ( kind = 8 ) a real ( kind = 8 ) b integer bin real ( kind = 8 ) cmax real ( kind = 8 ) cmin integer nbin ! ! Compute the bin limits. ! if ( bin < 1 ) then cmin = - huge ( cmin ) cmax = a else if ( bin <= nbin ) then cmin = ( real ( nbin - bin + 1, kind = 8 ) * a & + real ( bin - 1, kind = 8 ) * b ) & / real ( nbin, kind = 8 ) cmax = ( real ( nbin - bin, kind = 8 ) * a & + real ( bin, kind = 8 ) * b ) & / real ( nbin, kind = 8 ) else if ( nbin < bin ) then cmin = b cmax = huge ( cmax ) end if return end subroutine bin_to_r82_even ( nbin, bin, a, b, cmin, cmax ) !*****************************************************************************80 ! !! BIN_TO_R82_EVEN returns the limits for a given R82 "bin" in [A,B]. ! ! Discussion: ! ! The interval from A to B is divided into NBIN-2 equal subintervals or bins. ! An initial bin takes everything less than A, and a final bin takes ! everything greater than B. ! ! Example: ! ! NBIN = 7, A(1) = 5, B(1) = 15 ! A(2) = 0, B(2) = 20 ! ! BIN CMIN CMAX ! ------ ----------- -------- ! 1, 1 -HUGE -HUGE 5 0 ! 2, 2 5 0 7 4 ! 3, 3 7 4 9 8 ! 4, 4 9 8 11 12 ! 5, 5 11 12 13 16 ! 6, 6 13 16 15 20 ! 7, 7 15 20 HUGE HUGE ! ! Modified: ! ! 08 December 2004 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer NBIN, the number of bins. NBIN is normally at least 3. ! If NBIN is 1 or 2, then everything is assigned to bin 1. ! ! Input, integer BIN(2), the index of the bin to be considered. ! If BIN(I) is less than 1, or greater than NBIN, the user will get what ! the user deserves. ! ! Input, real ( kind = 8 ) A(2), B(2), the lower and upper limits of the ! bin interval. While A(I) is expected to be less than B(I), the code ! should return useful results if A(I) is actually greater than B(I). ! ! Output, real ( kind = 8 ) CMIN(2), CMAX(2), the minimum and maximum ! limits on the bin. ! implicit none integer, parameter :: ndim = 2 real ( kind = 8 ) a(ndim) real ( kind = 8 ) b(ndim) integer bin(ndim) real ( kind = 8 ) cmax(ndim) real ( kind = 8 ) cmin(ndim) integer i integer nbin do i = 1, ndim call bin_to_r8_even ( nbin, bin(i), a(i), b(i), cmin(i), cmax(i) ) end do return end subroutine bin_to_r82_even2 ( nbin, bin, a, b, cmin, cmax ) !*****************************************************************************80 ! !! BIN_TO_R82_EVEN2 returns the limits for a given D82 "bin" in [A,B]. ! ! Discussion: ! ! The interval from A to B is divided into NBIN equal subintervals or bins. ! ! Example: ! ! NBIN = 5, A(1) = 5, B(1) = 15 ! A(2) = 0, B(2) = 20 ! ! BIN CMIN CMAX ! ------ ----------- -------- ! 1, 1 5 0 7 4 ! 2, 2 7 4 9 8 ! 3, 3 9 8 11 12 ! 4, 4 11 12 13 16 ! 5, 5 13 16 15 20 ! ! Modified: ! ! 07 February 2001 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer NBIN, the number of bins. ! ! Input, integer BIN(2), the index of the bin to be considered. ! ! Input, real ( kind = 8 ) A(2), B(2), the lower and upper limits of the ! bin interval. While A(I) is expected to be less than B(I), the code ! should return useful results if A(I) is actually greater than B(I). ! ! Output, real ( kind = 8 ) CMIN(2), CMAX(2), the minimum and maximum ! limits on the bin. ! implicit none integer, parameter :: ndim = 2 real ( kind = 8 ) a(ndim) real ( kind = 8 ) b(ndim) integer bin(ndim) real ( kind = 8 ) cmax(ndim) real ( kind = 8 ) cmin(ndim) integer i integer nbin do i = 1, ndim call bin_to_r8_even2 ( nbin, bin(i), a(i), b(i), cmin(i), cmax(i) ) end do return end subroutine bin_to_r82_even3 ( nbin, bin, a, b, cmin, cmax ) !*****************************************************************************80 ! !! BIN_TO_R82_EVEN3 returns the limits for a given R82 "bin" in [A,B]. ! ! Discussion: ! ! The interval from A(I) to B(I) is divided into NBIN(I) equal ! subintervals or bins. ! ! Example: ! ! NBIN = (/ 4, 5, /) ! ! A(1) = 5, B(1) = 15 ! A(2) = 0, B(2) = 20 ! ! BIN CMIN CMAX ! ------ ----------- -------- ! 1, 1 5 0 7 4 ! 2, 2 7 4 9 8 ! 3, 3 9 8 11 12 ! 4, 4 11 12 13 16 ! 5, 5 13 16 15 20 ! ! Modified: ! ! 18 March 2001 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer NBIN(2), the number of bins in each dimension. ! ! Input, integer BIN(2), the index of the bin to be considered. ! ! Input, real ( kind = 8 ) A(2), B(2), the lower and upper limits of the ! bin interval. While A(I) is expected to be less than B(I), the code ! should return useful results if A(I) is actually greater than B(I). ! ! Output, real ( kind = 8 ) CMIN(2), CMAX(2), the minimum and maximum ! limits on the bin. ! implicit none integer, parameter :: ndim = 2 real ( kind = 8 ) a(ndim) real ( kind = 8 ) b(ndim) integer bin(ndim) real ( kind = 8 ) cmax(ndim) real ( kind = 8 ) cmin(ndim) integer i integer nbin(ndim) do i = 1, ndim call bin_to_r8_even2 ( nbin(i), bin(i), a(i), b(i), cmin(i), cmax(i) ) end do return end subroutine bin_to_r83_even2 ( nbin, bin, a, b, cmin, cmax ) !*****************************************************************************80 ! !! BIN_TO_R83_EVEN2 returns the limits for a given R83 "bin" in [A,B]. ! ! Discussion: ! ! The interval from A to B is divided into NBIN equal subintervals or bins. ! ! Modified: ! ! 09 February 2001 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer NBIN, the number of bins. ! ! Input, integer BIN(), the index of the bin to be considered. ! ! Input, real ( kind = 8 ) A(3), B(3), the lower and upper limits of the ! bin interval. While A(I) is expected to be less than B(I), the code ! should return useful results if A(I) is actually greater than B(I). ! ! Output, real ( kind = 8 ) CMIN(3), CMAX(3), the minimum and maximum ! limits on the bin. ! implicit none integer, parameter :: ndim = 3 real ( kind = 8 ) a(ndim) real ( kind = 8 ) b(ndim) integer bin(ndim) real ( kind = 8 ) cmax(ndim) real ( kind = 8 ) cmin(ndim) integer i integer nbin do i = 1, ndim call bin_to_r8_even2 ( nbin, bin(i), a(i), b(i), cmin(i), cmax(i) ) end do return end subroutine bin_to_r83_even3 ( nbin, bin, a, b, cmin, cmax ) !*****************************************************************************80 ! !! BIN_TO_R83_EVEN3 returns the limits for a given R83 "bin" in [A,B]. ! ! Discussion: ! ! The interval from A(I) to B(I) is divided into NBIN(I) equal ! subintervals or bins. ! ! Modified: ! ! 30 March 2001 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer NBIN(3), the number of bins in each dimension. ! ! Input, integer BIN(3), the index of the bin to be considered. ! ! Input, real ( kind = 8 ) A(3), B(3), the lower and upper limits of the ! bin interval. While A(I) is expected to be less than B(I), the code ! should return useful results if A(I) is actually greater than B(I). ! ! Output, real ( kind = 8 ) CMIN(3), CMAX(3), the minimum and maximum ! limits on the bin. ! implicit none integer, parameter :: ndim = 3 real ( kind = 8 ) a(ndim) real ( kind = 8 ) b(ndim) integer bin(ndim) real ( kind = 8 ) cmax(ndim) real ( kind = 8 ) cmin(ndim) integer i integer nbin(ndim) do i = 1, ndim call bin_to_r8_even2 ( nbin(i), bin(i), a(i), b(i), cmin(i), cmax(i) ) end do return end subroutine bin_to_r8vec_even3 ( ndim, nbin, bin, a, b, cmin, cmax ) !*****************************************************************************80 ! !! BIN_TO_R8VEC_EVEN3 returns the limits for a given R8VEC "bin" in [A,B]. ! ! Discussion: ! ! The interval from A(I) to B(I) is divided into NBIN(I) equal ! subintervals or bins. ! ! Modified: ! ! 08 December 2004 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer NDIM, the dimension of the space. ! ! Input, integer NBIN(NDIM), the number of bins in each dimension. ! ! Input, integer BIN(NDIM), the index of the bin to be considered. ! ! Input, real ( kind = 8 ) A(NDIM), B(NDIM), the lower and upper limits ! of the bin interval. While A(I) is expected to be less than B(I), ! the code should return useful results if A(I) is actually greater ! than B(I). ! ! Output, real ( kind = 8 ) CMIN(NDIM), CMAX(NDIM), the minimum and ! maximum limits on the bin. ! implicit none integer ndim real ( kind = 8 ) a(ndim) real ( kind = 8 ) b(ndim) integer bin(ndim) real ( kind = 8 ) cmax(ndim) real ( kind = 8 ) cmin(ndim) integer i integer nbin(ndim) do i = 1, ndim call bin_to_r8_even2 ( nbin(i), bin(i), a(i), b(i), cmin(i), cmax(i) ) end do return end function diaedg ( x0, y0, x1, y1, x2, y2, x3, y3 ) !*****************************************************************************80 ! !! DIAEDG chooses a diagonal edge. ! ! Discussion: ! ! The routine determines whether 0--2 or 1--3 is the diagonal edge ! that should be chosen, based on the circumcircle criterion, where ! (X0,Y0), (X1,Y1), (X2,Y2), (X3,Y3) are the vertices of a simple ! quadrilateral in counterclockwise order. ! ! Modified: ! ! 19 February 2001 ! ! Author: ! ! Barry Joe, ! Department of Computing Science, ! University of Alberta, ! Edmonton, Alberta, Canada T6G 2H1 ! ! Reference: ! ! Barry Joe, ! GEOMPACK - a software package for the generation of meshes ! using geometric algorithms, ! Advances in Engineering Software, ! Volume 13, pages 325-331, 1991. ! ! Parameters: ! ! Input, real ( kind = 8 ) X0, Y0, X1, Y1, X2, Y2, X3, Y3, the ! coordinates of the vertices of a quadrilateral, given in ! counter clockwise order. ! ! Output, integer DIAEDG, chooses a diagonal: ! +1, if diagonal edge 02 is chosen; ! -1, if diagonal edge 13 is chosen; ! 0, if the four vertices are cocircular. ! implicit none real ( kind = 8 ) ca real ( kind = 8 ) cb integer diaedg real ( kind = 8 ) dx10 real ( kind = 8 ) dx12 real ( kind = 8 ) dx30 real ( kind = 8 ) dx32 real ( kind = 8 ) dy10 real ( kind = 8 ) dy12 real ( kind = 8 ) dy30 real ( kind = 8 ) dy32 real ( kind = 8 ) s real ( kind = 8 ) tol real ( kind = 8 ) tola real ( kind = 8 ) tolb real ( kind = 8 ) x0 real ( kind = 8 ) x1 real ( kind = 8 ) x2 real ( kind = 8 ) x3 real ( kind = 8 ) y0 real ( kind = 8 ) y1 real ( kind = 8 ) y2 real ( kind = 8 ) y3 tol = 100.0D+00 * epsilon ( tol ) dx10 = x1 - x0 dy10 = y1 - y0 dx12 = x1 - x2 dy12 = y1 - y2 dx30 = x3 - x0 dy30 = y3 - y0 dx32 = x3 - x2 dy32 = y3 - y2 tola = tol * max ( abs ( dx10 ), abs ( dy10 ), abs ( dx30 ), abs ( dy30 ) ) tolb = tol * max ( abs ( dx12 ), abs ( dy12 ), abs ( dx32 ), abs ( dy32 ) ) ca = dx10 * dx30 + dy10 * dy30 cb = dx12 * dx32 + dy12 * dy32 if ( tola < ca .and. tolb < cb ) then diaedg = -1 else if ( ca < -tola .and. cb < -tolb ) then diaedg = 1 else tola = max ( tola, tolb ) s = ( dx10 * dy30 - dx30 * dy10 ) * cb + ( dx32 * dy12 - dx12 * dy32 ) * ca if ( tola < s ) then diaedg = -1 else if ( s < -tola ) then diaedg = 1 else diaedg = 0 end if end if return end subroutine dtris2 ( point_num, point_xy, tri_num, tri_vert, tri_nabe ) !*****************************************************************************80 ! !! DTRIS2 constructs a Delaunay triangulation of 2D vertices. ! ! Discussion: ! ! The routine constructs the Delaunay triangulation of a set of 2D vertices ! using an incremental approach and diagonal edge swaps. Vertices are ! first sorted in lexicographically increasing (X,Y) order, and ! then are inserted one at a time from outside the convex hull. ! ! Modified: ! ! 25 August 2001 ! ! Author: ! ! Barry Joe, ! Department of Computing Science, ! University of Alberta, ! Edmonton, Alberta, Canada T6G 2H1 ! ! Reference: ! ! Barry Joe, ! GEOMPACK - a software package for the generation of meshes ! using geometric algorithms, ! Advances in Engineering Software, ! Volume 13, pages 325-331, 1991. ! ! Parameters: ! ! Input, integer POINT_NUM, the number of vertices. ! ! Input/output, real ( kind = 8 ) POINT_XY(2,POINT_NUM), the coordinates ! of the vertices. On output, the vertices have been sorted into ! dictionary order. ! ! Output, integer TRI_NUM, the number of triangles in the triangulation; ! TRI_NUM is equal to 2*POINT_NUM - NB - 2, where NB is the number ! of boundary vertices. ! ! Output, integer TRI_VERT(3,TRI_NUM), the nodes that make up each triangle. ! The elements are indices of POINT_XY. The vertices of the triangles are ! in counter clockwise order. ! ! Output, integer TRI_NABE(3,TRI_NUM), the triangle neighbor list. ! Positive elements are indices of TIL; negative elements are used for links ! of a counter clockwise linked list of boundary edges; LINK = -(3*I + J-1) ! where I, J = triangle, edge index; TRI_NABE(J,I) refers to ! the neighbor along edge from vertex J to J+1 (mod 3). ! implicit none integer point_num real ( kind = 8 ) cmax integer e integer i integer ierr integer indx(point_num) integer j integer k integer l integer ledg integer lr integer lrline integer ltri integer m integer m1 integer m2 integer n real ( kind = 8 ) point_xy(2,point_num) integer redg integer rtri integer stack(point_num) integer t real ( kind = 8 ) tol integer top integer tri_nabe(3,point_num*2) integer tri_num integer tri_vert(3,point_num*2) tol = 100.0D+00 * epsilon ( tol ) ierr = 0 ! ! Sort the vertices by increasing (x,y). ! call r82vec_sort_heap_index_a ( point_num, point_xy, indx ) call r82vec_permute ( point_num, point_xy, indx ) ! ! Make sure that the data points are "reasonably" distinct. ! m1 = 1 do i = 2, point_num m = m1 m1 = i k = 0 do j = 1, 2 cmax = max ( abs ( point_xy(j,m) ), abs ( point_xy(j,m1) ) ) if ( tol * ( cmax + 1.0D+00 ) & < abs ( point_xy(j,m) - point_xy(j,m1) ) ) then k = j exit end if end do if ( k == 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'DTRIS2 - Fatal error!' write ( *, '(a,i6)' ) ' Fails for point number I = ', i write ( *, '(a,i6)' ) ' M = ', m write ( *, '(a,i6)' ) ' M1 = ', m1 write ( *, '(a,2g14.6)' ) ' X,Y(M) = ', point_xy(1,m), point_xy(2,m) write ( *, '(a,2g14.6)' ) ' X,Y(M1) = ', point_xy(1,m1), point_xy(2,m1) ierr = 224 return end if end do ! ! Starting from points M1 and M2, search for a third point M that ! makes a "healthy" triangle (M1,M2,M) ! m1 = 1 m2 = 2 j = 3 do if ( point_num < j ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'DTRIS2 - Fatal error!' ierr = 225 return end if m = j lr = lrline ( point_xy(1,m), point_xy(2,m), point_xy(1,m1), & point_xy(2,m1), point_xy(1,m2), point_xy(2,m2), 0.0D+00 ) if ( lr /= 0 ) then exit end if j = j + 1 end do ! ! Set up the triangle information for (M1,M2,M), and for any other ! triangles you created because points were collinear with M1, M2. ! tri_num = j - 2 if ( lr == -1 ) then tri_vert(1,1) = m1 tri_vert(2,1) = m2 tri_vert(3,1) = m tri_nabe(3,1) = -3 do i = 2, tri_num m1 = m2 m2 = i+1 tri_vert(1,i) = m1 tri_vert(2,i) = m2 tri_vert(3,i) = m tri_nabe(1,i-1) = -3 * i tri_nabe(2,i-1) = i tri_nabe(3,i) = i - 1 end do tri_nabe(1,tri_num) = -3 * tri_num - 1 tri_nabe(2,tri_num) = -5 ledg = 2 ltri = tri_num else tri_vert(1,1) = m2 tri_vert(2,1) = m1 tri_vert(3,1) = m tri_nabe(1,1) = -4 do i = 2, tri_num m1 = m2 m2 = i+1 tri_vert(1,i) = m2 tri_vert(2,i) = m1 tri_vert(3,i) = m tri_nabe(3,i-1) = i tri_nabe(1,i) = -3 * i - 3 tri_nabe(2,i) = i - 1 end do tri_nabe(3,tri_num) = -3 * tri_num tri_nabe(2,1) = -3 * tri_num - 2 ledg = 2 ltri = 1 end if ! ! Insert the vertices one at a time from outside the convex hull, ! determine visible boundary edges, and apply diagonal edge swaps until ! Delaunay triangulation of vertices (so far) is obtained. ! top = 0 do i = j+1, point_num m = i m1 = tri_vert(ledg,ltri) if ( ledg <= 2 ) then m2 = tri_vert(ledg+1,ltri) else m2 = tri_vert(1,ltri) end if lr = lrline ( point_xy(1,m), point_xy(2,m), point_xy(1,m1), & point_xy(2,m1), point_xy(1,m2), point_xy(2,m2), 0.0D+00 ) if ( 0 < lr ) then rtri = ltri redg = ledg ltri = 0 else l = -tri_nabe(ledg,ltri) rtri = l / 3 redg = mod(l,3) + 1 end if call vbedg ( point_xy(1,m), point_xy(2,m), point_num, point_xy, tri_num, & tri_vert, tri_nabe, ltri, ledg, rtri, redg ) n = tri_num + 1 l = -tri_nabe(ledg,ltri) do t = l / 3 e = mod ( l, 3 ) + 1 l = -tri_nabe(e,t) m2 = tri_vert(e,t) if ( e <= 2 ) then m1 = tri_vert(e+1,t) else m1 = tri_vert(1,t) end if tri_num = tri_num + 1 tri_nabe(e,t) = tri_num tri_vert(1,tri_num) = m1 tri_vert(2,tri_num) = m2 tri_vert(3,tri_num) = m tri_nabe(1,tri_num) = t tri_nabe(2,tri_num) = tri_num - 1 tri_nabe(3,tri_num) = tri_num + 1 top = top + 1 if ( point_num < top ) then ierr = 8 write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'DTRIS2 - Fatal error!' write ( *, '(a)' ) ' Stack overflow.' return end if stack(top) = tri_num if ( t == rtri .and. e == redg ) then exit end if end do tri_nabe(ledg,ltri) = -3 * n - 1 tri_nabe(2,n) = -3 * tri_num - 2 tri_nabe(3,tri_num) = -l ltri = n ledg = 2 call swapec ( m, top, ltri, ledg, point_num, point_xy, tri_num, & tri_vert, tri_nabe, stack, ierr ) if ( ierr /= 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'DTRIS2 - Fatal error!' write ( *, '(a)' ) ' Error return from SWAPEC.' return end if end do ! ! Now account for the sorting that we did. ! do i = 1, 3 do j = 1, tri_num tri_vert(i,j) = indx ( tri_vert(i,j) ) end do end do call perm_inv ( point_num, indx ) call r82vec_permute ( point_num, point_xy, indx ) return end subroutine get_seed ( seed ) !*****************************************************************************80 ! !! GET_SEED returns a seed for the random number generator. ! ! Discussion: ! ! The seed depends on the current time, and ought to be (slightly) ! different every millisecond. Once the seed is obtained, a random ! number generator should be called a few times to further process ! the seed. ! ! Modified: ! ! 17 November 2004 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Output, integer SEED, a pseudorandom seed value. ! implicit none integer seed real ( kind = 8 ) temp character ( len = 10 ) time character ( len = 8 ) today integer values(8) character ( len = 5 ) zone call date_and_time ( today, time, zone, values ) temp = 0.0D+00 temp = temp + real ( values(2) - 1, kind = 8 ) / 11.0D+00 temp = temp + real ( values(3) - 1, kind = 8 ) / 30.0D+00 temp = temp + real ( values(5), kind = 8 ) / 23.0D+00 temp = temp + real ( values(6), kind = 8 ) / 59.0D+00 temp = temp + real ( values(7), kind = 8 ) / 59.0D+00 temp = temp + real ( values(8), kind = 8 ) / 999.0D+00 temp = temp / 6.0D+00 ! ! Force 0 < TEMP <= 1. ! do while ( temp <= 0.0D+00 ) temp = temp + 1.0D+00 end do do while ( 1.0D+00 < temp ) temp = temp - 1.0D+00 end do seed = int ( real ( huge ( 1 ), kind = 8 ) * temp ) ! ! Never use a seed of 0 or maximum integer. ! if ( seed == 0 ) then seed = 1 end if if ( seed == huge ( 1 ) ) then seed = seed - 1 end if return end function i4_modp ( i, j ) !*****************************************************************************80 ! !! I4_MODP returns the nonnegative remainder of integer division. ! ! Discussion: ! ! If ! NREM = I4_MODP ( I, J ) ! NMULT = ( I - NREM ) / J ! then ! I = J * NMULT + NREM ! where NREM is always nonnegative. ! ! The MOD function computes a result with the same sign as the ! quantity being divided. Thus, suppose you had an angle A, ! and you wanted to ensure that it was between 0 and 360. ! Then mod(A,360) would do, if A was positive, but if A ! was negative, your result would be between -360 and 0. ! ! On the other hand, I4_MODP(A,360) is between 0 and 360, always. ! ! Examples: ! ! I J MOD I4_MODP Factorization ! ! 107 50 7 7 107 = 2 * 50 + 7 ! 107 -50 7 7 107 = -2 * -50 + 7 ! -107 50 -7 43 -107 = -3 * 50 + 43 ! -107 -50 -7 43 -107 = 3 * -50 + 43 ! ! Modified: ! ! 02 March 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer I, the number to be divided. ! ! Input, integer J, the number that divides I. ! ! Output, integer I4_MODP, the nonnegative remainder when I is ! divided by J. ! implicit none integer i integer i4_modp integer j if ( j == 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'I4_MODP - Fatal error!' write ( *, '(a,i6)' ) ' I4_MODP ( I, J ) called with J = ', j stop end if i4_modp = mod ( i, j ) if ( i4_modp < 0 ) then i4_modp = i4_modp + abs ( j ) end if return end subroutine i4_swap ( i, j ) !*****************************************************************************80 ! !! I4_SWAP switches two integer values. ! ! Modified: ! ! 30 November 1998 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input/output, integer I, J. On output, the values of I and ! J have been interchanged. ! implicit none integer i integer j integer k k = i i = j j = k return end function i4_wrap ( ival, ilo, ihi ) !*****************************************************************************80 ! !! I4_WRAP forces an integer to lie between given limits by wrapping. ! ! Example: ! ! ILO = 4, IHI = 8 ! ! I I4_WRAP ! ! -2 8 ! -1 4 ! 0 5 ! 1 6 ! 2 7 ! 3 8 ! 4 4 ! 5 5 ! 6 6 ! 7 7 ! 8 8 ! 9 4 ! 10 5 ! 11 6 ! 12 7 ! 13 8 ! 14 4 ! ! Modified: ! ! 19 August 2003 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer IVAL, an integer value. ! ! Input, integer ILO, IHI, the desired bounds for the integer value. ! ! Output, integer I4_WRAP, a "wrapped" version of IVAL. ! implicit none integer i4_modp integer i4_wrap integer ihi integer ilo integer ival integer jhi integer jlo integer wide jlo = min ( ilo, ihi ) jhi = max ( ilo, ihi ) wide = jhi - jlo + 1 if ( wide == 1 ) then i4_wrap = jlo else i4_wrap = jlo + i4_modp ( ival - jlo, wide ) end if return end subroutine i4mat_print ( m, n, a, title ) !*****************************************************************************80 ! !! I4MAT_PRINT prints an integer matrix. ! ! Modified: ! ! 30 June 2003 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer M, the number of rows in A. ! ! Input, integer N, the number of columns in A. ! ! Input, integer A(M,N), the matrix to be printed. ! ! Input, character ( len = * ) TITLE, a title to be printed first. ! TITLE may be blank. ! implicit none integer m integer n integer a(m,n) integer ihi integer ilo integer jhi integer jlo character ( len = * ) title ilo = 1 ihi = m jlo = 1 jhi = n call i4mat_print_some ( m, n, a, ilo, jlo, ihi, jhi, title ) return end subroutine i4mat_print_some ( m, n, a, ilo, jlo, ihi, jhi, title ) !*****************************************************************************80 ! !! I4MAT_PRINT_SOME prints some of an integer matrix. ! ! Modified: ! ! 04 November 2003 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer M, N, the number of rows and columns. ! ! Input, integer A(M,N), an M by N matrix to be printed. ! ! Input, integer ILO, JLO, the first row and column to print. ! ! Input, integer IHI, JHI, the last row and column to print. ! ! Input, character ( len = * ) TITLE, an optional title. ! implicit none integer, parameter :: incx = 10 integer m integer n integer a(m,n) character ( len = 7 ) ctemp(incx) integer i integer i2hi integer i2lo integer ihi integer ilo integer inc integer j integer j2 integer j2hi integer j2lo integer jhi integer jlo character ( len = * ) title if ( 0 < len_trim ( title ) ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) trim ( title ) end if do j2lo = max ( jlo, 1 ), min ( jhi, n ), incx j2hi = j2lo + incx - 1 j2hi = min ( j2hi, n ) j2hi = min ( j2hi, jhi ) inc = j2hi + 1 - j2lo write ( *, '(a)' ) ' ' do j = j2lo, j2hi j2 = j + 1 - j2lo write ( ctemp(j2), '(i7)') j end do write ( *, '('' Col '',10a7)' ) ctemp(1:inc) write ( *, '(a)' ) ' Row' write ( *, '(a)' ) ' ' i2lo = max ( ilo, 1 ) i2hi = min ( ihi, m ) do i = i2lo, i2hi do j2 = 1, inc j = j2lo - 1 + j2 write ( ctemp(j2), '(i7)' ) a(i,j) end do write ( *, '(i5,1x,10a7)' ) i, ( ctemp(j), j = 1, inc ) end do end do write ( *, '(a)' ) ' ' return end subroutine i4mat_transpose_print ( m, n, a, title ) !*****************************************************************************80 ! !! I4MAT_TRANSPOSE_PRINT prints an IMAT, transposed. ! ! Modified: ! ! 28 December 2004 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer M, N, the number of rows and columns. ! ! Input, integer A(M,N), an M by N matrix to be printed. ! ! Input, character ( len = * ) TITLE, an optional title. ! implicit none integer m integer n integer a(m,n) character ( len = * ) title call i4mat_transpose_print_some ( m, n, a, 1, 1, m, n, title ) return end subroutine i4mat_transpose_print_some ( m, n, a, ilo, jlo, ihi, jhi, title ) !*****************************************************************************80 ! !! I4MAT_TRANSPOSE_PRINT_SOME prints some of the transpose of an integer matrix. ! ! Modified: ! ! 28 December 2004 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer M, N, the number of rows and columns. ! ! Input, integer A(M,N), an M by N matrix to be printed. ! ! Input, integer ILO, JLO, the first row and column to print. ! ! Input, integer IHI, JHI, the last row and column to print. ! ! Input, character ( len = * ) TITLE, an optional title. ! implicit none integer, parameter :: incx = 10 integer m integer n integer a(m,n) character ( len = 7 ) ctemp(incx) integer i integer i2 integer i2hi integer i2lo integer ihi integer ilo integer inc integer j integer j2hi integer j2lo integer jhi integer jlo character ( len = * ) title if ( 0 < len_trim ( title ) ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) trim ( title ) end if do i2lo = max ( ilo, 1 ), min ( ihi, m ), incx i2hi = i2lo + incx - 1 i2hi = min ( i2hi, m ) i2hi = min ( i2hi, ihi ) inc = i2hi + 1 - i2lo write ( *, '(a)' ) ' ' do i = i2lo, i2hi i2 = i + 1 - i2lo write ( ctemp(i2), '(i7)') i end do write ( *, '('' Row '',10a7)' ) ctemp(1:inc) write ( *, '(a)' ) ' Col' write ( *, '(a)' ) ' ' j2lo = max ( jlo, 1 ) j2hi = min ( jhi, m ) do j = j2lo, j2hi do i2 = 1, inc i = i2lo - 1 + i2 write ( ctemp(i2), '(i7)' ) a(i,j) end do write ( *, '(i5,1x,10a7)' ) j, ( ctemp(i), i = 1, inc ) end do end do write ( *, '(a)' ) ' ' return end subroutine index_box2_next_2d ( n1, n2, ic, jc, i, j, more ) !*****************************************************************************80 ! !! INDEX_BOX2_NEXT_2D produces indices on the surface of a box in 2D. ! ! Discussion: ! ! The box has center at (IC,JC), and has half-widths N1 and N2. ! The indices are exactly those which are between (IC-N1,JC-N2) and ! (IC+N1,JC+N2) with the property that at least one of I and J ! is an "extreme" value. ! ! Modified: ! ! 02 May 2001 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer N1, N2, the half-widths of the box, that is, the ! maximum distance allowed between (IC,JC) and (I,J). ! ! Input, integer IC, JC, the central cell of the box. ! ! Input/output, integer I, J. On input, the previous index set. ! On output, the next index set. On the first call, MORE should ! be set to FALSE, and the input values of I and J are ignored. ! ! Input/output, logical MORE. ! On the first call for a given box, the user should set MORE to FALSE. ! On return, the routine sets MORE to TRUE. ! When there are no more indices, the routine sets MORE to FALSE. ! implicit none integer i integer ic integer j integer jc logical more integer n1 integer n2 if ( .not. more ) then more = .true. i = ic - n1 j = jc - n2 return end if if ( i == ic + n1 .and. & j == jc + n2 ) then more = .false. return end if ! ! Increment J. ! j = j + 1 ! ! Check J. ! if ( jc + n2 < j ) then j = jc - n2 i = i + 1 else if ( j < jc + n2 .and. & ( i == ic - n1 .or. i == ic + n1 ) ) then return else j = jc + n2 return end if return end subroutine index_box2_next_3d ( n1, n2, n3, ic, jc, kc, i, j, k, more ) !*****************************************************************************80 ! !! INDEX_BOX2_NEXT_3D produces indices on the surface of a box in 3D. ! ! Discussion: ! ! The box has a central cell of (IC,JC,KC), with a half widths of ! (N1,N2,N3). The indices are exactly those between (IC-N1,JC-N2,KC-N3) and ! (IC+N1,JC+N2,KC+N3) with the property that at least one of I, J, and K ! is an "extreme" value. ! ! Modified: ! ! 02 May 2001 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer N1, N2, N3, the "half widths" of the box, that is, the ! maximum distances from the central cell allowed for I, J and K. ! ! Input, integer IC, JC, KC, the central cell of the box. ! ! Input/output, integer I, J, K. On input, the previous index set. ! On output, the next index set. On the first call, MORE should ! be set to FALSE, and the input values of I, J, and K are ignored. ! ! Input/output, logical MORE. ! On the first call for a given box, the user should set MORE to FALSE. ! On return, the routine sets MORE to TRUE. ! When there are no more indices, the routine sets MORE to FALSE. ! implicit none integer i integer ic integer j integer jc integer k integer kc logical more integer n1 integer n2 integer n3 if ( .not. more ) then more = .true. i = ic - n1 j = jc - n2 k = kc - n3 return end if if ( i == ic + n1 .and. & j == jc + n2 .and. & k == kc + n3 ) then more = .false. return end if ! ! Increment K. ! k = k + 1 ! ! Check K. ! if ( kc + n3 < k ) then k = kc - n3 j = j + 1 else if ( k < kc + n3 .and. & ( i == ic - n1 .or. i == ic + n1 .or. & j == jc - n2 .or. j == jc + n2 ) ) then return else k = kc + n3 return end if ! ! Check J. ! if ( jc + n2 < j ) then j = jc - n2 i = i + 1 else if ( j < jc + n2 .and. & ( i == ic - n1 .or. i == ic + n1 .or. & k == kc - n3 .or. k == kc + n3 ) ) then return else j = jc + n2 return end if return end subroutine i4vec_heap_d ( n, a ) !*****************************************************************************80 ! !! I4VEC_HEAP_D reorders an I4VEC into a descending heap. ! ! Discussion: ! ! A descending heap is an array A with the property that, for every index J, ! A(2*J) <= A(J) and A(2*J+1) <= A(J), (as long as the indices ! 2*J and 2*J+1 are legal). ! ! Diagram: ! ! A(1) ! / \ ! A(2) A(3) ! / \ / \ ! A(4) A(5) A(6) A(7) ! / \ / \ ! A(8) A(9) A(10) A(11) ! ! Modified: ! ! 15 April 1999 ! ! Author: ! ! John Burkardt ! ! Reference: ! ! Albert Nijenhuis, Herbert Wilf, ! Combinatorial Algorithms, ! Academic Press, 1978, second edition, ! ISBN 0-12-519260-6. ! ! Parameters: ! ! Input, integer N, the size of the input array. ! ! Input/output, integer A(N). ! On input, an unsorted array. ! On output, the array has been reordered into a heap. ! implicit none integer n integer a(n) integer i integer ifree integer key integer m ! ! Only nodes N/2 down to 1 can be "parent" nodes. ! do i = n/2, 1, -1 ! ! Copy the value out of the parent node. ! Position IFREE is now "open". ! key = a(i) ifree = i do ! ! Positions 2*IFREE and 2*IFREE + 1 are the descendants of position ! IFREE. (One or both may not exist because they exceed N.) ! m = 2 * ifree ! ! Does the first position exist? ! if ( n < m ) then exit end if ! ! Does the second position exist? ! if ( m + 1 <= n ) then ! ! If both positions exist, take the larger of the two values, ! and update M if necessary. ! if ( a(m) < a(m+1) ) then m = m + 1 end if end if ! ! If the large descendant is larger than KEY, move it up, ! and update IFREE, the location of the free position, and ! consider the descendants of THIS position. ! if ( a(m) <= key ) then exit end if a(ifree) = a(m) ifree = m end do ! ! Once there is no more shifting to do, KEY moves into the free spot IFREE. ! a(ifree) = key end do return end subroutine i4vec_indicator ( n, a ) !*****************************************************************************80 ! !! I4VEC_INDICATOR sets an I4VEC to the indicator vector. ! ! Modified: ! ! 09 November 2000 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer N, the number of elements of A. ! ! Output, integer A(N), the array to be initialized. ! implicit none integer n integer a(n) integer i do i = 1, n a(i) = i end do return end subroutine i4vec_print ( n, a, title ) !*****************************************************************************80 ! !! I4VEC_PRINT prints an I4VEC. ! ! Modified: ! ! 28 November 2000 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer N, the number of components of the vector. ! ! Input, integer A(N), the vector to be printed. ! ! Input, character ( len = * ) TITLE, a title to be printed first. ! TITLE may be blank. ! implicit none integer n integer a(n) integer big integer i character ( len = * ) title if ( 0 < len_trim ( title ) ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) trim ( title ) end if big = maxval ( abs ( a(1:n) ) ) write ( *, '(a)' ) ' ' if ( big < 1000 ) then do i = 1, n write ( *, '(2x,i6,2x,i4)' ) i, a(i) end do else if ( big < 1000000 ) then do i = 1, n write ( *, '(2x,i6,2x,i7)' ) i, a(i) end do else do i = 1, n write ( *, '(2x,i6,2x,i12)' ) i, a(i) end do end if return end subroutine i4vec_sort_heap_a ( n, a ) !*****************************************************************************80 ! !! I4VEC_SORT_HEAP_A ascending sorts an I4VEC using heap sort. ! ! Modified: ! ! 15 April 1999 ! ! Author: ! ! John Burkardt ! ! Reference: ! ! Albert Nijenhuis, Herbert Wilf, ! Combinatorial Algorithms, ! Academic Press, 1978, second edition, ! ISBN 0-12-519260-6. ! ! Parameters: ! ! Input, integer N, the number of entries in the array. ! ! Input/output, integer A(N). ! On input, the array to be sorted; ! On output, the array has been sorted. ! implicit none integer n integer a(n) integer n1 if ( n <= 1 ) then return end if ! ! 1: Put A into descending heap form. ! call i4vec_heap_d ( n, a ) ! ! 2: Sort A. ! ! The largest object in the heap is in A(1). ! Move it to position A(N). ! call i4_swap ( a(1), a(n) ) ! ! Consider the diminished heap of size N1. ! do n1 = n-1, 2, -1 ! ! Restore the heap structure of A(1) through A(N1). ! call i4vec_heap_d ( n1, a ) ! ! Take the largest object from A(1) and move it to A(N1). ! call i4_swap ( a(1), a(n1) ) end do return end subroutine i4vec_unique ( n, a, nuniq ) !*****************************************************************************80 ! !! I4VEC_UNIQUE finds the number of unique elements in a sorted I4VEC. ! ! Modified: ! ! 09 July 2000 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer N, the number of elements in A. ! ! Input/output, integer A(N). On input, the sorted ! integer array. On output, the unique elements in A. ! ! Output, integer NUNIQ, the number of unique elements in A. ! implicit none integer n integer a(n) integer itest integer nuniq nuniq = 0 if ( n <= 0 ) then return end if nuniq = 1 do itest = 2, n if ( a(itest) /= a(nuniq) ) then nuniq = nuniq + 1 a(nuniq) = a(itest) end if end do return end subroutine i4vec2_compare ( n, a1, a2, i, j, isgn ) !*****************************************************************************80 ! !! I4VEC2_COMPARE compares pairs of integers stored in two vectors. ! ! Modified: ! ! 22 October 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer N, the number of data items. ! ! Input, integer A1(N), A2(N), contain the two components of each item. ! ! Input, integer I, J, the items to be compared. ! ! Output, integer ISGN, the results of the comparison: ! -1, item I < item J, ! 0, item I = item J, ! +1, item J < item I. ! implicit none integer n integer a1(n) integer a2(n) integer i integer isgn integer j isgn = 0 if ( a1(i) < a1(j) ) then isgn = -1 else if ( a1(i) == a1(j) ) then if ( a2(i) < a2(j) ) then isgn = -1 else if ( a2(i) < a2(j) ) then isgn = 0 else if ( a2(j) < a2(i) ) then isgn = +1 end if else if ( a1(j) < a1(i) ) then isgn = +1 end if return end subroutine i4vec2_sort_a ( n, a1, a2 ) !*****************************************************************************80 ! !! I4VEC2_SORT_A ascending sorts a vector of pairs of integers. ! ! Discussion: ! ! Each item to be sorted is a pair of integers (I,J), with the I ! and J values stored in separate vectors A1 and A2. ! ! Modified: ! ! 27 June 2000 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer N, the number of items of data. ! ! Input/output, integer A1(N), A2(N), the data to be sorted.. ! implicit none integer n integer a1(n) integer a2(n) integer i integer indx integer isgn integer j ! ! Initialize. ! i = 0 indx = 0 isgn = 0 j = 0 ! ! Call the external heap sorter. ! do call sort_heap_external ( n, indx, i, j, isgn ) ! ! Interchange the I and J objects. ! if ( 0 < indx ) then call i4_swap ( a1(i), a1(j) ) call i4_swap ( a2(i), a2(j) ) ! ! Compare the I and J objects. ! else if ( indx < 0 ) then call i4vec2_compare ( n, a1, a2, i, j, isgn ) else if ( indx == 0 ) then exit end if end do return end subroutine i4vec2_unique ( n, a1, a2, nuniq ) !*****************************************************************************80 ! !! I4VEC2_UNIQUE keeps the unique elements in a array of pairs of integers. ! ! Discussion: ! ! Item I is stored as the pair A1(I), A2(I). ! ! The items must have been sorted, or at least it must be the ! case that equal items are stored in adjacent vector locations. ! ! If the items were not sorted, then this routine will only ! replace a string of equal values by a single representative. ! ! Modified: ! ! 09 July 2000 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer N, the number of items. ! ! Input/output, integer A1(N), A2(N). ! On input, the array of N items. ! On output, an array of NUNIQ unique items. ! ! Output, integer NUNIQ, the number of unique items. ! implicit none integer n integer a1(n) integer a2(n) integer itest integer nuniq nuniq = 0 if ( n <= 0 ) then return end if nuniq = 1 do itest = 2, n if ( a1(itest) /= a1(nuniq) .or. a2(itest) /= a2(nuniq) ) then nuniq = nuniq + 1 a1(nuniq) = a1(itest) a2(nuniq) = a2(itest) end if end do return end function lrline ( xu, yu, xv1, yv1, xv2, yv2, dv ) !*****************************************************************************80 ! !! LRLINE determines where a point lies in relation to a directed line. ! ! Discussion: ! ! LRLINE determines whether a point is to the left of, right of, ! or on a directed line parallel to a line through given points. ! ! Modified: ! ! 18 June 2001 ! ! Author: ! ! Barry Joe, ! Department of Computing Science, ! University of Alberta, ! Edmonton, Alberta, Canada T6G 2H1 ! ! Reference: ! ! Barry Joe, ! GEOMPACK - a software package for the generation of meshes ! using geometric algorithms, ! Advances in Engineering Software, ! Volume 13, pages 325-331, 1991. ! ! Parameters: ! ! Input, real ( kind = 8 ) XU, YU, XV1, YV1, XV2, YV2, are vertex ! coordinates; the directed line is parallel to and at signed distance ! DV to the left of the directed line from (XV1,YV1) to (XV2,YV2); (XU,YU) ! is the vertex for which the position relative to the directed line is ! to be determined. ! ! Input, real ( kind = 8 ) DV, the signed distance, positive for left. ! ! Output, integer LRLINE, is +1, 0, or -1 depending on whether (XU,YU) is ! to the right of, on, or left of the directed line. LRLINE is 0 if ! the line degenerates to a point. ! implicit none real ( kind = 8 ) dv real ( kind = 8 ) dx real ( kind = 8 ) dxu real ( kind = 8 ) dy real ( kind = 8 ) dyu integer lrline real ( kind = 8 ) t real ( kind = 8 ), parameter :: tol = 0.0000001D+00 real ( kind = 8 ) tolabs real ( kind = 8 ) xu real ( kind = 8 ) xv1 real ( kind = 8 ) xv2 real ( kind = 8 ) yu real ( kind = 8 ) yv1 real ( kind = 8 ) yv2 dx = xv2 - xv1 dy = yv2 - yv1 dxu = xu - xv1 dyu = yu - yv1 tolabs = tol * max ( abs ( dx ), abs ( dy ), abs ( dxu ), abs ( dyu ), & abs ( dv ) ) t = dy * dxu - dx * dyu + dv * sqrt ( dx * dx + dy * dy ) if ( tolabs < t ) then lrline = 1 else if ( -tolabs <= t ) then lrline = 0 else if ( t < -tolabs ) then lrline = -1 end if return end subroutine perm_inv ( n, p ) !*****************************************************************************80 ! !! PERM_INV inverts a permutation "in place". ! ! Modified: ! ! 25 July 2000 ! ! Parameters: ! ! Input, integer N, the number of objects being permuted. ! ! Input/output, integer P(N), the permutation, in standard index form. ! On output, P describes the inverse permutation ! implicit none integer n integer i integer i0 integer i1 integer i2 integer is integer p(n) if ( n <= 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'PERM_INV - Fatal error!' write ( *, '(a,i6)' ) ' Input value of N = ', n stop end if is = 1 do i = 1, n i1 = p(i) do while ( i < i1 ) i2 = p(i1) p(i1) = -i2 i1 = i2 end do is = -sign ( 1, p(i) ) p(i) = sign ( p(i), is ) end do do i = 1, n i1 = -p(i) if ( 0 <= i1 ) then i0 = i do i2 = p(i1) p(i1) = i0 if ( i2 < 0 ) then exit end if i0 = i1 i1 = i2 end do end if end do return end subroutine points_nearest_point_bins_2d ( nset, pset, nbin, bin_min, bin_max, & bin_start, bin_last, bin_next, p, i_min, dist_min, compares ) !*****************************************************************************80 ! !! POINTS_NEAREST_POINT_BINS_2D finds the nearest point to a given point in 2D. ! ! Discussion: ! ! A set of NSET points with coordinates PSET is assumed to lie within a ! rectangle. The limits of this rectangle are given in BIN_MIN and BIN_MAX. ! The rectangle is divided up into an NBIN by NBIN regular grid of cells. ! ! The cells are ordered lexicographically, as suggested by the following ! diagram for NBIN = 5: ! ! *----*----*----*----*----* ! | | | | | | ! | 21 | 22 | 23 | 24 | 25 | ! | | | | | | ! *----*----*----*----*----* ! | | | | | | ! | 16 | 17 | 18 | 19 | 20 | ! | | | | | | ! *----*----*----*----*----* ! | | | | | | ! | 11 | 12 | 13 | 14 | 15 | ! | | | | | | ! *----*----*----*----*----* ! | | | | | | ! | 6 | 7 | 8 | 9 | 10 | ! | | | | | | ! *----*----*----*----*----* ! | | | | | | ! | 1 | 2 | 3 | 4 | 5 | ! | | | | | | ! *----*----*----*----*----* ! ! The points in PSET are ordered by cell, and within a cell they ! are ordered lexicographically. Thus, for points P1 and P2, ! ! P1 < P2 implies that ! * P1 is in a lower ordered cell than P2, or ! * P1 is in the same cell as P2, but P1.X < P2.X, or ! * P1 is in the same cell as P2, P1.X = P2.X, but P1.Y < P2.Y. ! ! The arrays BIN_START(I,J) and BIN_LAST(I,J) are given the coordinates ! I, J of a cell, and return the lowest and highest index in PSET of a ! point in the I, J cell. All indices in between also belong to ! this cell. If the cell has no points, then both arrays are -1. ! ! ! After all this preprocessing, the algorithm for finding the nearest ! point goes as follows: ! ! 1) for a point P, determine the cell it belongs to. ! Note that this algorithm will NOT be valid if P lies outside ! the bin limits. ! ! 2) Search for a cell that has at least one member of PSET in it. ! We start at the cell containing P, but if there are no members ! there, we spiral out from the cell, one layer at a time. ! ! 3) Within this cell, find the point nearest to P. We now know that ! we don't need to search any cell whose points will all be further ! from P than this distance. ! ! 4) Now, search in all other cells that could have a point closer to ! P than what we have found so far. ! ! Modified: ! ! 15 February 2001 ! ! Author: ! ! John Burkardt ! ! Reference: ! ! Jon Bentley, Bruce Weide, Andrew Yao, ! Optimal Expected Time Algorithms for Closest Point Problems, ! ACM Transactions on Mathematical Software, ! Volume 6, Number 4, December 1980, pages 563-580. ! ! Parameters: ! ! Input, integer NSET, the number of points in the set. ! ! Input, real ( kind = 8 ) PSET(2,NSET), the points in the set. ! ! Input, integer NBIN, the number of cells. ! ! Input, real ( kind = 8 ) BIN_MIN(2), BIN_MAX(2), the minimum and ! maximum bin values. ! ! Input, integer BIN_START(NBIN,NBIN), BIN_LAST(NBIN,NBIN), indicates ! the index of the first and last element in the bin, or -1 if none. ! ! Input, integer BIN_NEXT(NSET), the index of the next element of the bin ! containing this element. ! ! Input, real ( kind = 8 ) P(2), the point to be tested. ! ! Output, integer I_MIN, the index of the nearest point in PSET to P. ! ! Output, real ( kind = 8 ) DIST_MIN, the distance between P and PSET(*,I_MIN). ! ! Output, integer COMPARES, the number of point-to-point comparisons. ! implicit none integer nbin integer, parameter :: ndim = 2 integer nset integer bin(ndim) integer bin_last(nbin,nbin) real ( kind = 8 ) bin_max(ndim) real ( kind = 8 ) bin_min(ndim) integer bin_start(nbin,nbin) integer bin_next(nset) real ( kind = 8 ) c_max(ndim) real ( kind = 8 ) c_min(ndim) integer compares real ( kind = 8 ) dist_min real ( kind = 8 ) dist_min_sq real ( kind = 8 ) dist_sq integer i integer i_min integer ic integer il integer j integer jc integer jl integer layer real ( kind = 8 ) layer_width logical more_bins integer node real ( kind = 8 ) p(ndim) real ( kind = 8 ) pset(ndim,nset) real ( kind = 8 ) search_radius compares = 0 ! ! Special cases. ! if ( nset <= 0 ) then dist_min = huge ( dist_min ) i_min = 0 return end if if ( nset == 1 ) then dist_min = sqrt ( sum ( ( p(1:ndim) - pset(1:ndim,1) )**2 ) ) compares = 1 i_min = 1 return end if ! ! Initialize. ! dist_min_sq = huge ( dist_min_sq ) i_min = 0 search_radius = 0.0D+00 layer_width = minval ( bin_max(1:ndim) - bin_min(1:ndim) ) & / real ( nbin, kind = 8 ) ! ! Determine the bin coordinates of the point P. ! call r82_to_bin_even2 ( nbin, bin_min, bin_max, p, bin ) ! ! Determine the radius of the ball of space that will be completely ! searched after this first bin, "layer 0", is done. ! call bin_to_r82_even2 ( nbin, bin, bin_min, bin_max, c_min, c_max ) ! ! Set the central bin of the layers. ! ic = bin(1) jc = bin(2) ! ! Set ! * the current layer, ! * the starting bin of the current layer, ! * the current bin ! layer = 0 il = ic jl = jc i = il j = jl do ! ! Search all legal bins in layer LAYER. ! do ! ! Search BIN I, J. ! if ( 1 <= i .and. i <= nbin .and. 1 <= j .and. j <= nbin ) then node = bin_start(i,j) do while ( 0 < node ) dist_sq = sum ( ( p(1:ndim) - pset(1:ndim,node) )**2 ) compares = compares + 1 if ( dist_sq < dist_min_sq ) then dist_min_sq = dist_sq i_min = node end if node = bin_next(node) end do end if ! ! We have now searched all points in bin I, J. ! ! Determine the next bin in this layer. ! ! But if it lies outside the region, discard it, and go on to the next one. ! Once you get past the last bin, exit because you're done the layer. ! more_bins = .true. do if ( i < ic + layer .and. j == jc - layer ) then i = i + 1 else if ( i == ic + layer .and. j < jc + layer ) then j = j + 1 else if ( ic - layer < i .and. j == jc + layer ) then i = i - 1 else if ( i == ic - layer .and. jc - layer + 1 < j ) then j = j - 1 else more_bins = .false. exit end if if ( 1 <= i .and. i <= nbin .and. & 1 <= j .and. j <= nbin ) then exit end if end do if ( .not. more_bins ) then exit end if end do ! ! We've completed this layer. ! Update the radius of the searched area. ! if ( layer == 0 ) then search_radius = min ( & minval ( abs ( p(1:ndim) - c_min(1:ndim) ) ), & minval ( abs ( p(1:ndim) - c_max(1:ndim) ) ) ) else search_radius = search_radius + layer_width end if ! ! We can stop if: ! ! * We've found at least one neighbor; ! AND ! * We've searched all legal boxes that contain the circle ! with P at the center and the nearest N on the circumference. ! if ( i_min /= 0 ) then dist_min = sqrt ( dist_min_sq ) if ( dist_min <= search_radius ) then exit end if end if ! ! Prepare to search the next layer. ! layer = layer + 1 il = ic - layer jl = jc - layer i = il j = jl end do return end subroutine points_nearest_point_bins_2d_2 ( nset, pset, nbin, bin_min, bin_max, & bin_start, bin_last, bin_next, ptest, i_min, dist_min, compares ) !*****************************************************************************80 ! !! POINTS_NEAREST_POINT_BINS_2D_2 finds the nearest point to a given point in 2D. ! ! Discussion: ! ! This code differs from POINTS_NEAREST_POINT_BINS_2D by calling ! a subroutine to compute the next bin index. ! ! A set of NSET points with coordinates PSET is assumed to lie within a ! rectangle. The limits of this rectangle are given in BIN_MIN and BIN_MAX. ! The rectangle is divided up into an NBIN by NBIN regular grid of cells. ! ! The cells are ordered lexicographically, as suggested by the following ! diagram for NBIN = 5: ! ! *----*----*----*----*----* ! | | | | | | ! | 21 | 22 | 23 | 24 | 25 | ! | | | | | | ! *----*----*----*----*----* ! | | | | | | ! | 16 | 17 | 18 | 19 | 20 | ! | | | | | | ! *----*----*----*----*----* ! | | | | | | ! | 11 | 12 | 13 | 14 | 15 | ! | | | | | | ! *----*----*----*----*----* ! | | | | | | ! | 6 | 7 | 8 | 9 | 10 | ! | | | | | | ! *----*----*----*----*----* ! | | | | | | ! | 1 | 2 | 3 | 4 | 5 | ! | | | | | | ! *----*----*----*----*----* ! ! The points in PSET are ordered by cell, and within a cell they ! are ordered lexicographically. Thus, for points P1 and P2, ! ! P1 < P2 implies that ! * P1 is in a lower ordered cell than P2, or ! * P1 is in the same cell as P2, but P1.X < P2.X, or ! * P1 is in the same cell as P2, P1.X = P2.X, but P1.Y < P2.Y. ! ! The arrays BIN_START(I,J) and BIN_LAST(I,J) are given the coordinates ! I, J of a cell, and return the lowest and highest index in PSET of a ! point in the I, J cell. All indices in between also belong to ! this cell. If the cell has no points, then both arrays are -1. ! ! ! After all this preprocessing, the algorithm for finding the nearest ! point goes as follows: ! ! 1) for a point PTEST, determine the cell it belongs to. ! Note that this algorithm will NOT be valid if PTEST lies outside ! the bin limits. ! ! 2) Search for a cell that has at least one member of PSET in it. ! We start at the cell containing PTEST, but if there are no members ! there, we spiral out from the cell, one layer at a time. ! ! 3) Within this cell, find the point nearest to PTEST. We now know that ! we don't need to search any cell whose points will all be further ! from PTEST than this distance. ! ! 4) Now, search in all other cells that could have a point closer to ! PTEST than what we have found so far. ! ! Modified: ! ! 08 February 2001 ! ! Author: ! ! John Burkardt ! ! Reference: ! ! Jon Bentley, Bruce Weide, Andrew Yao, ! Optimal Expected Time Algorithms for Closest Point Problems, ! ACM Transactions on Mathematical Software, ! Volume 6, Number 4, December 1980, pages 563-580. ! ! Parameters: ! ! Input, integer NSET, the number of points in the set. ! ! Input, real ( kind = 8 ) PSET(2,NSET), the points in the set. ! ! Input, integer NBIN, the number of cells. ! ! Input, real ( kind = 8 ) BIN_MIN(2), BIN_MAX(2), the minimum and ! maximum bin values. ! ! Input, integer BIN_START(NBIN,NBIN), BIN_LAST(NBIN,NBIN), indicates ! the index of the first and last element in the bin, or -1 if none. ! ! Input, integer BIN_NEXT(NSET), the index of the next element of the bin ! containing this element. ! ! Input, real ( kind = 8 ) PTEST(2), the coordinates of the test points. ! ! Output, integer I_MIN, the index of the nearest point in PSET to PTEST. ! ! Output, real ( kind = 8 ) DIST_MIN, the distance between PTEST and ! PSET(*,I_MIN). ! ! Output, integer COMPARES, the number of point-to-point comparisons. ! implicit none integer nbin integer, parameter :: ndim = 2 integer nset integer bin(ndim) integer bin_last(nbin,nbin) real ( kind = 8 ) bin_max(ndim) real ( kind = 8 ) bin_min(ndim) integer bin_start(nbin,nbin) integer bin_next(nset) real ( kind = 8 ) c_max(ndim) real ( kind = 8 ) c_min(ndim) integer compares real ( kind = 8 ) dist_min real ( kind = 8 ) dist_min_sq real ( kind = 8 ) dist_sq integer i integer i_min integer ic integer j integer jc integer layer real ( kind = 8 ) layer_width logical more_bins integer node real ( kind = 8 ) pset(ndim,nset) real ( kind = 8 ) ptest(ndim) real ( kind = 8 ) search_radius compares = 0 ! ! Special cases. ! if ( nset <= 0 ) then dist_min = huge ( dist_min ) i_min = 0 return end if if ( nset == 1 ) then dist_min = sqrt ( sum ( ( ptest(1:ndim) - pset(1:ndim,1) )**2 ) ) compares = 1 i_min = 1 return end if layer_width = minval ( bin_max(1:ndim) - bin_min(1:ndim) ) & / real ( nbin, kind = 8 ) ! ! Search for each of the NTEST points. ! dist_min_sq = huge ( dist_min_sq ) i_min = 0 search_radius = 0.0D+00 ! ! Determine the bin coordinates of the point P. ! call r82_to_bin_even2 ( nbin, bin_min, bin_max, ptest, bin ) ! ! Determine the radius of the ball of space that will be completely ! searched after this first bin, "layer 0", is done. ! call bin_to_r82_even2 ( nbin, bin, bin_min, bin_max, c_min, c_max ) ! ! Set the central bin of the layers. ! ic = bin(1) jc = bin(2) layer = 0 ! ! Search all legal bins in layer LAYER. ! do more_bins = .false. call index_box2_next_2d ( layer, layer, ic, jc, i, j, more_bins ) ! ! In layer LAYER, search each BIN I, J. ! do if ( 1 <= i .and. i <= nbin .and. 1 <= j .and. j <= nbin ) then node = bin_start(i,j) do while ( 0 < node ) dist_sq = sum ( ( ptest(1:ndim) - pset(1:ndim,node) )**2 ) compares = compares + 1 if ( dist_sq < dist_min_sq ) then dist_min_sq = dist_sq i_min = node end if node = bin_next(node) end do end if ! ! We have now searched all points in bin I, J. ! ! Determine the next bin in this layer. ! ! But if it lies outside the region, discard it, and go on to the next one. ! Once you get past the last bin, exit because you're done the layer. ! do call index_box2_next_2d ( layer, layer, ic, jc, i, j, more_bins ) if ( .not. more_bins ) then exit end if if ( 1 <= i .and. i <= nbin .and. & 1 <= j .and. j <= nbin ) then exit end if end do if ( .not. more_bins ) then exit end if end do ! ! We've completed layer LAYER. ! Update the radius of the searched area. ! if ( layer == 0 ) then search_radius = min ( & minval ( abs ( ptest(1:ndim) - c_min(1:ndim) ) ), & minval ( abs ( ptest(1:ndim) - c_max(1:ndim) ) ) ) else search_radius = search_radius + layer_width end if ! ! We are done with PTEST if: ! ! * We've found at least one neighbor; ! AND ! * We've searched all legal boxes that contain the circle ! with PTEST at the center and the nearest neighbor on the circumference. ! if ( i_min /= 0 ) then dist_min = sqrt ( dist_min_sq ) if ( dist_min <= search_radius ) then exit end if end if layer = layer + 1 end do ! ! We are now done with all the layers. ! return end subroutine points_nearest_point_bins_2d_3 ( nset, pset, nbin, bin_min, & bin_max, bin_start, bin_last, bin_next, ptest, i_min, dist_min, compares ) !*****************************************************************************80 ! !! POINTS_NEAREST_POINT_BINS_2D_3 finds the nearest point to a given point in 2D. ! ! Discussion: ! ! This code differs from POINTS_NEAREST_POINTS_BINS_2D by allowing the ! user to specify the number of bins in each dimension. ! ! A set of NSET points with coordinates PSET is assumed to lie within a ! rectangle. The limits of this rectangle are given in BIN_MIN and BIN_MAX. ! The rectangle is divided up into an NBIN(1) by NBIN(2) regular grid of ! cells. ! ! The cells are ordered lexicographically, as suggested by the following ! diagram for NBIN = (/ 5, 4 /) ! ! *----*----*----*----*----* ! | | | | | | ! | 16 | 17 | 18 | 19 | 20 | ! | | | | | | ! *----*----*----*----*----* ! | | | | | | ! | 11 | 12 | 13 | 14 | 15 | ! | | | | | | ! *----*----*----*----*----* ! | | | | | | ! | 6 | 7 | 8 | 9 | 10 | ! | | | | | | ! *----*----*----*----*----* ! | | | | | | ! | 1 | 2 | 3 | 4 | 5 | ! | | | | | | ! *----*----*----*----*----* ! ! The points in PSET are ordered by cell, and within a cell they ! are ordered lexicographically. Thus, for points P1 and P2, ! ! P1 < P2 implies that ! * P1 is in a lower ordered cell than P2, or ! * P1 is in the same cell as P2, but P1.X < P2.X, or ! * P1 is in the same cell as P2, P1.X = P2.X, but P1.Y < P2.Y. ! ! The arrays BIN_START(I,J) and BIN_LAST(I,J) are given the coordinates ! I, J of a cell, and return the lowest and highest index in PSET of a ! point in the I, J cell. All indices in between also belong to ! this cell. If the cell has no points, then both arrays are -1. ! ! ! After all this preprocessing, the algorithm for finding the nearest ! point goes as follows: ! ! 1) for a point PTEST, determine the cell it belongs to. ! Note that this algorithm will NOT be valid if PTEST lies outside ! the bin limits. ! ! 2) Search for a cell that has at least one member of PSET in it. ! We start at the cell containing PTEST, but if there are no members ! there, we spiral out from the cell, one layer at a time. ! ! 3) Within this cell, find the point nearest to PTEST. We now know that ! we don't need to search any cell whose points will all be further ! from PTEST than this distance. ! ! 4) Now, search in all other cells that could have a point closer to ! PTEST than what we have found so far. ! ! Modified: ! ! 18 March 2001 ! ! Author: ! ! John Burkardt ! ! Reference: ! ! Jon Bentley, Bruce Weide, Andrew Yao, ! Optimal Expected Time Algorithms for Closest Point Problems, ! ACM Transactions on Mathematical Software, ! Volume 6, Number 4, December 1980, pages 563-580. ! ! Parameters: ! ! Input, integer NSET, the number of points in the set. ! ! Input, real ( kind = 8 ) PSET(2,NSET), the points in the set. ! ! Input, integer NBIN(2), the number of cells in the horizontal and ! vertical directions. ! ! Input, real ( kind = 8 ) BIN_MIN(2), BIN_MAX(2), the minimum and ! maximum bin values. ! ! Input, integer BIN_START(NBIN(1),NBIN(2)), BIN_LAST(NBIN(1),NBIN(2)), ! indicates the index of the first and last element in the bin, ! or -1 if none. ! ! Input, integer BIN_NEXT(NSET), the index of the next element of the bin ! containing this element. ! ! Input, real ( kind = 8 ) PTEST(2), the coordinates of the test point. ! ! Output, integer I_MIN, the index of the nearest point in PSET to PTEST. ! ! Output, real ( kind = 8 ) DIST_MIN, the distance between PTEST and ! PSET(*,I_MIN). ! ! Output, integer COMPARES, the number of point-to-point comparisons. ! implicit none integer, parameter :: ndim = 2 integer nbin(ndim) integer nset integer bin(ndim) integer bin_last(nbin(1),nbin(2)) real ( kind = 8 ) bin_max(ndim) real ( kind = 8 ) bin_min(ndim) integer bin_start(nbin(1),nbin(2)) integer bin_next(nset) real ( kind = 8 ) c_max(ndim) real ( kind = 8 ) c_min(ndim) integer compares real ( kind = 8 ) dist_min real ( kind = 8 ) dist_min_sq real ( kind = 8 ) dist_sq integer i integer i_min integer ic integer j integer jc integer layer real ( kind = 8 ) layer_width logical more_bins integer node real ( kind = 8 ) pset(ndim,nset) real ( kind = 8 ) ptest(ndim) real ( kind = 8 ) search_radius compares = 0 ! ! Special cases. ! if ( nset <= 0 ) then dist_min = huge ( dist_min ) i_min = 0 return end if if ( nset == 1 ) then dist_min = sqrt ( sum ( ( ptest(1:ndim) - pset(1:ndim,1) )**2 ) ) i_min = 1 compares = 1 return end if ! ! The efficiency of the code will suffer if the data in the vector ! ! bin_max(1:ndim) - bin_min(1:ndim) / real ( nbin(1:ndim), kind = 8 ) ! ! varies significantly. ! layer_width = minval ( & bin_max(1:ndim) - bin_min(1:ndim) / real ( nbin(1:ndim), kind = 8 ) ) dist_min_sq = huge ( dist_min_sq ) i_min = 0 search_radius = 0.0D+00 ! ! Determine the bin coordinates of the point P. ! call r82_to_bin_even3 ( nbin, bin_min, bin_max, ptest, bin ) ! ! Determine the radius of the ball of space that will be completely ! searched after this first bin, "layer 0", is done. ! call bin_to_r82_even3 ( nbin, bin, bin_min, bin_max, c_min, c_max ) ! ! Set the central bin of the layers. ! ic = bin(1) jc = bin(2) layer = 0 ! ! Search all legal bins in layer LAYER. ! do more_bins = .false. call index_box2_next_2d ( layer, layer, ic, jc, i, j, more_bins ) ! ! In layer LAYER, search each BIN I, J. ! do if ( 1 <= i .and. i <= nbin(1) .and. 1 <= j .and. j <= nbin(2) ) then node = bin_start(i,j) do while ( 0 < node ) dist_sq = sum ( ( ptest(1:ndim) - pset(1:ndim,node) )**2 ) compares = compares + 1 if ( dist_sq < dist_min_sq ) then dist_min_sq = dist_sq i_min = node end if node = bin_next(node) end do end if ! ! We have now searched all points in bin I, J. ! ! Determine the next bin in this layer. ! ! But if it lies outside the region, discard it, and go on to the next one. ! Once you get past the last bin, exit because you're done the layer. ! do call index_box2_next_2d ( layer, layer, ic, jc, i, j, more_bins ) if ( .not. more_bins ) then exit end if if ( 1 <= i .and. i <= nbin(1) .and. & 1 <= j .and. j <= nbin(2) ) then exit end if end do if ( .not. more_bins ) then exit end if end do ! ! We've completed layer LAYER. ! Update the radius of the searched area. ! if ( layer == 0 ) then search_radius = min ( & minval ( abs ( ptest(1:ndim) - c_min(1:ndim) ) ), & minval ( abs ( ptest(1:ndim) - c_max(1:ndim) ) ) ) else search_radius = search_radius + layer_width end if ! ! We are done if: ! ! * We've found at least one neighbor; ! AND ! * We've searched all legal boxes that contain the circle ! with PTEST at the center and the nearest neighbor on the circumference. ! if ( i_min /= 0 ) then dist_min = sqrt ( dist_min_sq ) if ( dist_min <= search_radius ) then exit end if end if layer = layer + 1 end do return end subroutine points_nearest_point_del_2d ( point_num, xc, xd, nabes_first, & nabes_num, nabes_dim, nabes, nnear, dnear ) !*****************************************************************************80 ! !! POINTS_NEAREST_POINT_DEL_2D: nearest neighbor in a Delaunay triangulation. ! ! Discussion: ! ! A set of points XC is given, along with its Delaunay triangulation. ! Now a new point XD is given, and we need to know the closest point in XC. ! ! This algorithm starts at a random point in XC, and then repeatedly moves ! to a neighboring point that is closer to XD. This is guaranteed to be ! possible because the triangulation is Delaunay. Otherwise, it ! would be possible to reach a vertex which was not the closest, ! but for which all neighbors were further away. ! ! This algorithm is able to handle the case where the point XD lies ! outside the convex hull. ! ! The algorithm is very simple to code. In the most likely ! case, it should have an expected time complexity of O(sqrt(N)). ! ! Overhead occurs in the development of the vertex adjacency data ! structure. The size of this array should be roughly 6*N on average. ! Given the list of nodes that make up triangles, the vertex adjacency ! data can be constructed by storing every pair of nodes (I,J) and (J,I), ! and sorting the data into dictionary order. ! ! Modified: ! ! 20 July 2001 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer POINT_NUM, the number of points. ! ! Input, real ( kind = 8 ) XC(2,POINT_NUM), the set of points. ! ! Input, real ( kind = 8 ) XD(2), a point whose nearest neighbor is ! to be found. ! ! Input, integer NABES_FIRST(POINT_NUM), the index in NABES of the first ! neighbor in the list for each node. ! ! Input, integer NABES_NUM(POINT_NUM), the number of neighbors of each node. ! ! Input, integer NABES_DIM, the dimension of NABES. ! ! Input, integer NABES(NABES_DIM), a list of the neighbors of all the nodes. ! Neighbors of node 1 are listed first, and so on. ! ! Output, integer NNEAR, the nearest node to XD. ! ! Output, real ( kind = 8 ) DNEAR, the distance of the nearest node to XD. ! implicit none integer nabes_dim integer point_num real ( kind = 8 ) dist real ( kind = 8 ) dnear integer i integer i1 integer j integer nabes(nabes_dim) integer nabes_first(point_num) integer nabes_num(point_num) integer nnear integer nnear_old real ( kind = 8 ) x real ( kind = 8 ) x1 real ( kind = 8 ) xc(2,point_num) real ( kind = 8 ) xd(2) real ( kind = 8 ) y real ( kind = 8 ) y1 x = xd(1) y = xd(2) ! ! Select a random vertex. ! nnear = 1 x1 = xc(1,nnear) y1 = xc(2,nnear) dnear = sqrt ( ( x1 - x )**2 + ( y1 - y )**2 ) ! ! From the current vertex, consider all neighboring vertices. ! For each neighbor, compute the distance to the point. ! If no neighbor is closer, then the current vertex is the closest. ! Otherwise, set the current vertex to the neighbor that was closest, ! and repeat. ! do nnear_old = nnear j = nabes_first(nnear_old) do i = 1, nabes_num(nnear_old) i1 = nabes(j) x1 = xc(1,i1) y1 = xc(2,i1) dist = sqrt ( ( x1 - x )**2 + ( y1 - y )**2 ) if ( dist < dnear ) then dnear = dist nnear = i1 end if j = j + 1 end do ! ! If no neighbor was closer, we're done. ! if ( nnear == nnear_old ) then exit end if end do return end subroutine points_nearest_point_naive_2d ( nset, pset, ptest, i_min, dist_min ) !*****************************************************************************80 ! !! POINTS_NEAREST_POINT_NAIVE_2D finds the nearest point to a given point in 2D. ! ! Discussion: ! ! A naive algorithm is used. The distance to every point is calculated, ! in order to determine the smallest. ! ! Modified: ! ! 06 February 2001 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer NSET, the number of points in the set. ! ! Input, real ( kind = 8 ) PSET(2,NSET), the points in the set. ! ! Input, real ( kind = 8 ) PTEST(2), the point whose nearest neighbor ! is sought. ! ! Output, integer I_MIN, the index of the nearest point in PSET to P. ! ! Output, real ( kind = 8 ) DIST_MIN, the distance between P and PSET(*,I_MIN). ! implicit none integer nset integer, parameter :: ndim = 2 real ( kind = 8 ) d real ( kind = 8 ) dist_min integer i integer i_min real ( kind = 8 ) pset(ndim,nset) real ( kind = 8 ) ptest(ndim) dist_min = huge ( dist_min ) i_min = -1 do i = 1, nset d = sum ( ( ptest(1:ndim) - pset(1:ndim,i) )**2 ) if ( d < dist_min ) then dist_min = d i_min = i end if end do dist_min = sqrt ( dist_min ) return end subroutine points_nearest_point_naive_3d ( nset, pset, ptest, i_min, dist_min ) !*****************************************************************************80 ! !! POINTS_NEAREST_POINT_NAIVE_3D finds the nearest point to a given point in 3D. ! ! Discussion: ! ! A naive algorithm is used. The distance to every point is calculated, ! in order to determine the smallest. ! ! Modified: ! ! 09 February 2001 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer NSET, the number of points in the set. ! ! Input, real ( kind = 8 ) PSET(3,NSET), the points in the set. ! ! Input, real ( kind = 8 ) PTEST(3), the point whose nearest neighbor ! is sought. ! ! Output, integer I_MIN, the index of the nearest point in PSET to P. ! ! Output, real ( kind = 8 ) DIST_MIN, the distance between P and PSET(*,I_MIN). ! implicit none integer nset integer, parameter :: ndim = 3 real ( kind = 8 ) d real ( kind = 8 ) dist_min integer i integer i_min real ( kind = 8 ) pset(ndim,nset) real ( kind = 8 ) ptest(ndim) dist_min = huge ( dist_min ) i_min = -1 do i = 1, nset d = sum ( ( ptest(1:ndim) - pset(1:ndim,i) )**2 ) if ( d < dist_min ) then dist_min = d i_min = i end if end do dist_min = sqrt ( dist_min ) return end subroutine points_nearest_points_bins_2d ( nset, pset, nbin, bin_min, bin_max, & bin_start, bin_last, bin_next, ntest, ptest, i_min, dist_min, compares ) !*****************************************************************************80 ! !! POINTS_NEAREST_POINTS_BINS_2D finds the nearest point to given points in 2D. ! ! Discussion: ! ! A set of NSET points with coordinates PSET is assumed to lie within a ! rectangle. The limits of this rectangle are given in BIN_MIN and BIN_MAX. ! The rectangle is divided up into an NBIN by NBIN regular grid of cells. ! ! The cells are ordered lexicographically, as suggested by the following ! diagram for NBIN = 5: ! ! *----*----*----*----*----* ! | | | | | | ! | 21 | 22 | 23 | 24 | 25 | ! | | | | | | ! *----*----*----*----*----* ! | | | | | | ! | 16 | 17 | 18 | 19 | 20 | ! | | | | | | ! *----*----*----*----*----* ! | | | | | | ! | 11 | 12 | 13 | 14 | 15 | ! | | | | | | ! *----*----*----*----*----* ! | | | | | | ! | 6 | 7 | 8 | 9 | 10 | ! | | | | | | ! *----*----*----*----*----* ! | | | | | | ! | 1 | 2 | 3 | 4 | 5 | ! | | | | | | ! *----*----*----*----*----* ! ! The points in PSET are ordered by cell, and within a cell they ! are ordered lexicographically. Thus, for points P1 and P2, ! ! P1 < P2 implies that ! * P1 is in a lower ordered cell than P2, or ! * P1 is in the same cell as P2, but P1.X < P2.X, or ! * P1 is in the same cell as P2, P1.X = P2.X, but P1.Y < P2.Y. ! ! The arrays BIN_START(I,J) and BIN_LAST(I,J) are given the coordinates ! I, J of a cell, and return the lowest and highest index in PSET of a ! point in the I, J cell. All indices in between also belong to ! this cell. If the cell has no points, then both arrays are -1. ! ! ! After all this preprocessing, the algorithm for finding the nearest ! point goes as follows: ! ! 1) for a point PTEST, determine the cell it belongs to. ! Note that this algorithm will NOT be valid if PTEST lies outside ! the bin limits. ! ! 2) Search for a cell that has at least one member of PSET in it. ! We start at the cell containing PTEST, but if there are no members ! there, we spiral out from the cell, one layer at a time. ! ! 3) Within this cell, find the point nearest to PTEST. We now know that ! we don't need to search any cell whose points will all be further ! from PTEST than this distance. ! ! 4) Now, search in all other cells that could have a point closer to ! PTEST than what we have found so far. ! ! Modified: ! ! 08 February 2001 ! ! Author: ! ! John Burkardt ! ! Reference: ! ! Jon Bentley, Bruce Weide, Andrew Yao, ! Optimal Expected Time Algorithms for Closest Point Problems, ! ACM Transactions on Mathematical Software, ! Volume 6, Number 4, December 1980, pages 563-580. ! ! Parameters: ! ! Input, integer NSET, the number of points in the set. ! ! Input, real ( kind = 8 ) PSET(2,NSET), the points in the set. ! ! Input, integer NBIN, the number of cells. ! ! Input, real ( kind = 8 ) BIN_MIN(2), BIN_MAX(2), the minimum and ! maximum bin values. ! ! Input, integer BIN_START(NBIN,NBIN), BIN_LAST(NBIN,NBIN), indicates ! the index of the first and last element in the bin, or -1 if none. ! ! Input, integer BIN_NEXT(NSET), the index of the next element of the bin ! containing this element. ! ! Input, integer NTEST, the number of test points. ! ! Input, real ( kind = 8 ) PTEST(2,NTEST), the test points. ! ! Output, integer I_MIN(NTEST), the index of the nearest point in PSET ! to PTEST(ITEST). ! ! Output, real ( kind = 8 ) DIST_MIN(NTEST), the distance between ! PTEST(*,ITEST) and PSET(*,I_MIN). ! ! Output, integer COMPARES(NTEST), the number of point-to-point comparisons. ! implicit none integer nbin integer, parameter :: ndim = 2 integer nset integer ntest integer bin(ndim) integer bin_last(nbin,nbin) real ( kind = 8 ) bin_max(ndim) real ( kind = 8 ) bin_min(ndim) integer bin_start(nbin,nbin) integer bin_next(nset) real ( kind = 8 ) c_max(ndim) real ( kind = 8 ) c_min(ndim) integer compares(ntest) real ( kind = 8 ) dist_min(ntest) real ( kind = 8 ) dist_min_sq real ( kind = 8 ) dist_sq integer i integer i_min(ntest) integer ic integer il integer itest integer j integer jc integer jl integer layer real ( kind = 8 ) layer_width logical more_bins integer node real ( kind = 8 ) pset(ndim,nset) real ( kind = 8 ) ptest(ndim,ntest) real ( kind = 8 ) search_radius compares(1:ntest) = 0 ! ! Special cases. ! if ( nset <= 0 ) then dist_min(1:ntest) = huge ( dist_min ) i_min(1:ntest) = 0 return end if if ( nset == 1 ) then do itest = 1, ntest dist_min(itest) = sqrt ( & sum ( ( ptest(1:ndim,itest) - pset(1:ndim,1) )**2 ) ) end do compares(1:ntest) = 1 i_min(1:ntest) = 1 return end if layer_width = minval ( bin_max(1:ndim) - bin_min(1:ndim) ) & / real ( nbin, kind = 8 ) ! ! Search for each of the NTEST points. ! do itest = 1, ntest dist_min_sq = huge ( dist_min_sq ) i_min(itest) = 0 search_radius = 0.0D+00 ! ! Determine the bin coordinates of the point P. ! call r82_to_bin_even2 ( nbin, bin_min, bin_max, ptest(1,itest), bin ) ! ! Determine the radius of the ball of space that will be completely ! searched after this first bin, "layer 0", is done. ! call bin_to_r82_even2 ( nbin, bin, bin_min, bin_max, c_min, c_max ) ! ! Set the central bin of the layers. ! ic = bin(1) jc = bin(2) ! ! Set ! * the current layer, ! * the starting bin of the current layer, ! * the current bin ! layer = 0 il = ic jl = jc i = il j = jl do ! ! Search all legal bins in layer LAYER. ! do ! ! Search BIN I, J. ! if ( 1 <= i .and. i <= nbin .and. 1 <= j .and. j <= nbin ) then node = bin_start(i,j) do while ( 0 < node ) dist_sq = sum ( ( ptest(1:ndim,itest) - pset(1:ndim,node) )**2 ) compares(itest) = compares(itest) + 1 if ( dist_sq < dist_min_sq ) then dist_min_sq = dist_sq i_min(itest) = node end if node = bin_next(node) end do end if ! ! We have now searched all points in bin I, J. ! ! Determine the next bin in this layer. ! ! But if it lies outside the region, discard it, and go on to the next one. ! Once you get past the last bin, exit because you're done the layer. ! more_bins = .true. do if ( i < ic + layer .and. j == jc - layer ) then i = i + 1 else if ( i == ic + layer .and. j < jc + layer ) then j = j + 1 else if ( ic - layer < i .and. j == jc + layer ) then i = i - 1 else if ( i == ic - layer .and. jc - layer + 1 < j ) then j = j - 1 else more_bins = .false. exit end if if ( 1 <= i .and. i <= nbin .and. & 1 <= j .and. j <= nbin ) then exit end if end do if ( .not. more_bins ) then exit end if end do ! ! We've completed this layer. ! Update the radius of the searched area. ! if ( layer == 0 ) then search_radius = min ( & minval ( abs ( ptest(1:ndim,itest) - c_min(1:ndim) ) ), & minval ( abs ( ptest(1:ndim,itest) - c_max(1:ndim) ) ) ) else search_radius = search_radius + layer_width end if ! ! We are done with PTEST(*,ITEST) if: ! ! * We've found at least one neighbor; ! AND ! * We've searched all legal boxes that contain the circle ! with PTEST at the center and the nearest neighbor on the circumference. ! if ( i_min(itest) /= 0 ) then dist_min(itest) = sqrt ( dist_min_sq ) if ( dist_min(itest) <= search_radius ) then exit end if end if ! ! Prepare to search the next layer. ! layer = layer + 1 il = ic - layer jl = jc - layer i = il j = jl end do end do return end subroutine points_nearest_points_bins_2d_2 ( nset, pset, nbin, bin_min, & bin_max, bin_start, bin_last, bin_next, ntest, ptest, i_min, dist_min, compares ) !*****************************************************************************80 ! !! POINTS_NEAREST_POINTS_BINS_2D_2 finds the nearest point to given points in 2D. ! ! Discussion: ! ! This code differs from POINTS_NEAREST_POINTS_BINS_2D by calling ! a subroutine to compute the next bin index. ! ! A set of NSET points with coordinates PSET is assumed to lie within a ! rectangle. The limits of this rectangle are given in BIN_MIN and BIN_MAX. ! The rectangle is divided up into an NBIN by NBIN regular grid of cells. ! ! The cells are ordered lexicographically, as suggested by the following ! diagram for NBIN = 5: ! ! *----*----*----*----*----* ! | | | | | | ! | 21 | 22 | 23 | 24 | 25 | ! | | | | | | ! *----*----*----*----*----* ! | | | | | | ! | 16 | 17 | 18 | 19 | 20 | ! | | | | | | ! *----*----*----*----*----* ! | | | | | | ! | 11 | 12 | 13 | 14 | 15 | ! | | | | | | ! *----*----*----*----*----* ! | | | | | | ! | 6 | 7 | 8 | 9 | 10 | ! | | | | | | ! *----*----*----*----*----* ! | | | | | | ! | 1 | 2 | 3 | 4 | 5 | ! | | | | | | ! *----*----*----*----*----* ! ! The points in PSET are ordered by cell, and within a cell they ! are ordered lexicographically. Thus, for points P1 and P2, ! ! P1 < P2 implies that ! * P1 is in a lower ordered cell than P2, or ! * P1 is in the same cell as P2, but P1.X < P2.X, or ! * P1 is in the same cell as P2, P1.X = P2.X, but P1.Y < P2.Y. ! ! The arrays BIN_START(I,J) and BIN_LAST(I,J) are given the coordinates ! I, J of a cell, and return the lowest and highest index in PSET of a ! point in the I, J cell. All indices in between also belong to ! this cell. If the cell has no points, then both arrays are -1. ! ! ! After all this preprocessing, the algorithm for finding the nearest ! point goes as follows: ! ! 1) for a point PTEST, determine the cell it belongs to. ! Note that this algorithm will NOT be valid if PTEST lies outside ! the bin limits. ! ! 2) Search for a cell that has at least one member of PSET in it. ! We start at the cell containing PTEST, but if there are no members ! there, we spiral out from the cell, one layer at a time. ! ! 3) Within this cell, find the point nearest to PTEST. We now know that ! we don't need to search any cell whose points will all be further ! from PTEST than this distance. ! ! 4) Now, search in all other cells that could have a point closer to ! PTEST than what we have found so far. ! ! Modified: ! ! 08 February 2001 ! ! Author: ! ! John Burkardt ! ! Reference: ! ! Jon Bentley, Bruce Weide, Andrew Yao, ! Optimal Expected Time Algorithms for Closest Point Problems, ! ACM Transactions on Mathematical Software, ! Volume 6, Number 4, December 1980, pages 563-580. ! ! Parameters: ! ! Input, integer NSET, the number of points in the set. ! ! Input, real ( kind = 8 ) PSET(2,NSET), the points in the set. ! ! Input, integer NBIN, the number of cells. ! ! Input, real ( kind = 8 ) BIN_MIN(2), BIN_MAX(2), the minimum and ! maximum bin values. ! ! Input, integer BIN_START(NBIN,NBIN), BIN_LAST(NBIN,NBIN), indicates ! the index of the first and last element in the bin, or -1 if none. ! ! Input, integer BIN_NEXT(NSET), the index of the next element of the bin ! containing this element. ! ! Input, integer NTEST, the number of test points. ! ! Input, real ( kind = 8 ) PTEST(2,NTEST), the test points. ! ! Output, integer I_MIN(NTEST), the index of the nearest point in PSET ! to PTEST(ITEST). ! ! Output, real ( kind = 8 ) DIST_MIN(NTEST), the distance between ! PTEST(*,ITEST) and PSET(*,I_MIN). ! ! Output, integer COMPARES(NTEST), the number of point-to-point comparisons. ! implicit none integer nbin integer, parameter :: ndim = 2 integer nset integer ntest integer bin(ndim) integer bin_last(nbin,nbin) real ( kind = 8 ) bin_max(ndim) real ( kind = 8 ) bin_min(ndim) integer bin_start(nbin,nbin) integer bin_next(nset) real ( kind = 8 ) c_max(ndim) real ( kind = 8 ) c_min(ndim) integer compares(ntest) real ( kind = 8 ) dist_min(ntest) real ( kind = 8 ) dist_min_sq real ( kind = 8 ) dist_sq integer i integer i_min(ntest) integer ic integer itest integer j integer jc integer layer real ( kind = 8 ) layer_width logical more_bins integer node real ( kind = 8 ) pset(ndim,nset) real ( kind = 8 ) ptest(ndim,ntest) real ( kind = 8 ) search_radius compares(1:ntest) = 0 ! ! Special cases. ! if ( nset <= 0 ) then dist_min(1:ntest) = huge ( dist_min ) i_min(1:ntest) = 0 return end if if ( nset == 1 ) then do itest = 1, ntest dist_min(itest) = sqrt ( & sum ( ( ptest(1:ndim,itest) - pset(1:ndim,1) )**2 ) ) end do compares(1:ntest) = 1 i_min(1:ntest) = 1 return end if layer_width = minval ( bin_max(1:ndim) - bin_min(1:ndim) ) & / real ( nbin, kind = 8 ) ! ! Search for each of the NTEST points. ! do itest = 1, ntest dist_min_sq = huge ( dist_min_sq ) i_min(itest) = 0 search_radius = 0.0D+00 ! ! Determine the bin coordinates of the point P. ! call r82_to_bin_even2 ( nbin, bin_min, bin_max, ptest(1,itest), bin ) ! ! Determine the radius of the ball of space that will be completely ! searched after this first bin, "layer 0", is done. ! call bin_to_r82_even2 ( nbin, bin, bin_min, bin_max, c_min, c_max ) ! ! Set the central bin of the layers. ! ic = bin(1) jc = bin(2) layer = 0 ! ! Search all legal bins in layer LAYER. ! do more_bins = .false. call index_box2_next_2d ( layer, layer, ic, jc, i, j, more_bins ) ! ! In layer LAYER, search each BIN I, J. ! do if ( 1 <= i .and. i <= nbin .and. 1 <= j .and. j <= nbin ) then node = bin_start(i,j) do while ( 0 < node ) dist_sq = sum ( ( ptest(1:ndim,itest) - pset(1:ndim,node) )**2 ) compares(itest) = compares(itest) + 1 if ( dist_sq < dist_min_sq ) then dist_min_sq = dist_sq i_min(itest) = node end if node = bin_next(node) end do end if ! ! We have now searched all points in bin I, J. ! ! Determine the next bin in this layer. ! ! But if it lies outside the region, discard it, and go on to the next one. ! Once you get past the last bin, exit because you're done the layer. ! do call index_box2_next_2d ( layer, layer, ic, jc, i, j, more_bins ) if ( .not. more_bins ) then exit end if if ( 1 <= i .and. i <= nbin .and. & 1 <= j .and. j <= nbin ) then exit end if end do if ( .not. more_bins ) then exit end if end do ! ! We've completed layer LAYER. ! Update the radius of the searched area. ! if ( layer == 0 ) then search_radius = min ( & minval ( abs ( ptest(1:ndim,itest) - c_min(1:ndim) ) ), & minval ( abs ( ptest(1:ndim,itest) - c_max(1:ndim) ) ) ) else search_radius = search_radius + layer_width end if ! ! We are done with PTEST(*,ITEST) if: ! ! * We've found at least one neighbor; ! AND ! * We've searched all legal boxes that contain the circle ! with PTEST at the center and the nearest neighbor on the circumference. ! if ( i_min(itest) /= 0 ) then dist_min(itest) = sqrt ( dist_min_sq ) if ( dist_min(itest) <= search_radius ) then exit end if end if layer = layer + 1 end do ! ! We are now done with all the layers. ! end do ! ! We are now done with all the test points. ! return end subroutine points_nearest_points_bins_2d_3 ( nset, pset, nbin, bin_min, & bin_max, bin_start, bin_last, bin_next, ntest, ptest, i_min, dist_min, compares ) !*****************************************************************************80 ! !! POINTS_NEAREST_POINTS_BINS_2D_3 finds the nearest point to given points in 2D. ! ! Discussion: ! ! This code differs from POINTS_NEAREST_POINTS_BINS_2D by allowing the ! user to specify the number of bins in each dimension. ! ! A set of NSET points with coordinates PSET is assumed to lie within a ! rectangle. The limits of this rectangle are given in BIN_MIN and BIN_MAX. ! The rectangle is divided up into an NBIN(1) by NBIN(2) regular grid of ! cells. ! ! The cells are ordered lexicographically, as suggested by the following ! diagram for NBIN = (/ 5, 4 /) ! ! *----*----*----*----*----* ! | | | | | | ! | 16 | 17 | 18 | 19 | 20 | ! | | | | | | ! *----*----*----*----*----* ! | | | | | | ! | 11 | 12 | 13 | 14 | 15 | ! | | | | | | ! *----*----*----*----*----* ! | | | | | | ! | 6 | 7 | 8 | 9 | 10 | ! | | | | | | ! *----*----*----*----*----* ! | | | | | | ! | 1 | 2 | 3 | 4 | 5 | ! | | | | | | ! *----*----*----*----*----* ! ! The points in PSET are ordered by cell, and within a cell they ! are ordered lexicographically. Thus, for points P1 and P2, ! ! P1 < P2 implies that ! * P1 is in a lower ordered cell than P2, or ! * P1 is in the same cell as P2, but P1.X < P2.X, or ! * P1 is in the same cell as P2, P1.X = P2.X, but P1.Y < P2.Y. ! ! The arrays BIN_START(I,J) and BIN_LAST(I,J) are given the coordinates ! I, J of a cell, and return the lowest and highest index in PSET of a ! point in the I, J cell. All indices in between also belong to ! this cell. If the cell has no points, then both arrays are -1. ! ! ! After all this preprocessing, the algorithm for finding the nearest ! point goes as follows: ! ! 1) for a point PTEST, determine the cell it belongs to. ! Note that this algorithm will NOT be valid if PTEST lies outside ! the bin limits. ! ! 2) Search for a cell that has at least one member of PSET in it. ! We start at the cell containing PTEST, but if there are no members ! there, we spiral out from the cell, one layer at a time. ! ! 3) Within this cell, find the point nearest to PTEST. We now know that ! we don't need to search any cell whose points will all be further ! from PTEST than this distance. ! ! 4) Now, search in all other cells that could have a point closer to ! PTEST than what we have found so far. ! ! Modified: ! ! 18 March 2001 ! ! Author: ! ! John Burkardt ! ! Reference: ! ! Jon Bentley, Bruce Weide, Andrew Yao, ! Optimal Expected Time Algorithms for Closest Point Problems, ! ACM Transactions on Mathematical Software, ! Volume 6, Number 4, December 1980, pages 563-580. ! ! Parameters: ! ! Input, integer NSET, the number of points in the set. ! ! Input, real ( kind = 8 ) PSET(2,NSET), the points in the set. ! ! Input, integer NBIN(2), the number of cells in the horizontal and ! vertical directions. ! ! Input, real ( kind = 8 ) BIN_MIN(2), BIN_MAX(2), the minimum and ! maximum bin values. ! ! Input, integer BIN_START(NBIN(1),NBIN(2)), BIN_LAST(NBIN(1),NBIN(2)), ! indicates the index of the first and last element in the bin, ! or -1 if none. ! ! Input, integer BIN_NEXT(NSET), the index of the next element of the bin ! containing this element. ! ! Input, integer NTEST, the number of test points. ! ! Input, real ( kind = 8 ) PTEST(2,NTEST), the coordinates of the ! test points. ! ! Output, integer I_MIN(NTEST), the index of the nearest point in PSET ! to PTEST(ITEST). ! ! Output, real ( kind = 8 ) DIST_MIN(NTEST), the distance between ! PTEST(*,ITEST) and PSET(*,I_MIN). ! ! Output, integer COMPARES(NTEST), the number of point-to-point comparisons. ! implicit none integer, parameter :: ndim = 2 integer nbin(ndim) integer nset integer ntest integer bin(ndim) integer bin_last(nbin(1),nbin(2)) real ( kind = 8 ) bin_max(ndim) real ( kind = 8 ) bin_min(ndim) integer bin_start(nbin(1),nbin(2)) integer bin_next(nset) real ( kind = 8 ) c_max(ndim) real ( kind = 8 ) c_min(ndim) integer compares(ntest) real ( kind = 8 ) dist_min(ntest) real ( kind = 8 ) dist_min_sq real ( kind = 8 ) dist_sq integer i integer i_min(ntest) integer ic integer itest integer j integer jc integer layer real ( kind = 8 ) layer_width logical more_bins integer node real ( kind = 8 ) pset(ndim,nset) real ( kind = 8 ) ptest(ndim,ntest) real ( kind = 8 ) search_radius compares(1:ntest) = 0 ! ! Special cases. ! if ( nset <= 0 ) then dist_min(1:ntest) = huge ( dist_min(1) ) i_min(1:ntest) = 0 return end if if ( nset == 1 ) then do itest = 1, ntest dist_min(itest) = sqrt ( sum ( ( ptest(1:ndim,itest) - pset(1:ndim,1) )**2 ) ) end do compares(1:ntest) = 1 i_min(1:ntest) = 1 return end if ! ! The efficiency of the code will suffer if the data in the vector ! ! bin_max(1:ndim) - bin_min(1:ndim) / real ( nbin(1:ndim), kind = 8 ) ! ! varies significantly. ! layer_width = minval ( & bin_max(1:ndim) - bin_min(1:ndim) / real ( nbin(1:ndim), kind = 8 ) ) ! ! Search for each of the NTEST points. ! do itest = 1, ntest dist_min_sq = huge ( dist_min_sq ) i_min(itest) = 0 search_radius = 0.0D+00 ! ! Determine the bin coordinates of the point P. ! call r82_to_bin_even3 ( nbin, bin_min, bin_max, ptest(1,itest), bin ) ! ! Determine the radius of the ball of space that will be completely ! searched after this first bin, "layer 0", is done. ! call bin_to_r82_even3 ( nbin, bin, bin_min, bin_max, c_min, c_max ) ! ! Set the central bin of the layers. ! ic = bin(1) jc = bin(2) layer = 0 ! ! Search all legal bins in layer LAYER. ! do more_bins = .false. call index_box2_next_2d ( layer, layer, ic, jc, i, j, more_bins ) ! ! In layer LAYER, search each BIN I, J. ! do if ( 1 <= i .and. i <= nbin(1) .and. 1 <= j .and. j <= nbin(2) ) then node = bin_start(i,j) do while ( 0 < node ) dist_sq = sum ( ( ptest(1:ndim,itest) - pset(1:ndim,node) )**2 ) compares(itest) = compares(itest) + 1 if ( dist_sq < dist_min_sq ) then dist_min_sq = dist_sq i_min(itest) = node end if node = bin_next(node) end do end if ! ! We have now searched all points in bin I, J. ! ! Determine the next bin in this layer. ! ! But if it lies outside the region, discard it, and go on to the next one. ! Once you get past the last bin, exit because you're done the layer. ! do call index_box2_next_2d ( layer, layer, ic, jc, i, j, more_bins ) if ( .not. more_bins ) then exit end if if ( 1 <= i .and. i <= nbin(1) .and. & 1 <= j .and. j <= nbin(2) ) then exit end if end do if ( .not. more_bins ) then exit end if end do ! ! We've completed layer LAYER. ! Update the radius of the searched area. ! if ( layer == 0 ) then search_radius = min ( & minval ( abs ( ptest(1:ndim,itest) - c_min(1:ndim) ) ), & minval ( abs ( ptest(1:ndim,itest) - c_max(1:ndim) ) ) ) else search_radius = search_radius + layer_width end if ! ! We are done with PTEST(*,ITEST) if: ! ! * We've found at least one neighbor; ! AND ! * We've searched all legal boxes that contain the circle ! with PTEST at the center and the nearest neighbor on the circumference. ! if ( i_min(itest) /= 0 ) then dist_min(itest) = sqrt ( dist_min_sq ) if ( dist_min(itest) <= search_radius ) then exit end if end if layer = layer + 1 end do ! ! We are now done with all the layers. ! end do ! ! We are now done with all the test points. ! return end subroutine points_nearest_points_bins_2d_4 ( nset, pset, nbin, bin_min, & bin_max, bin_start, bin_last, bin_next, ntest, ptest, i_min, dist_min, compares ) !*****************************************************************************80 ! !! POINTS_NEAREST_POINTS_BINS_2D_4 finds the nearest point to given points in 2D. ! ! Discussion: ! ! This code differs from POINTS_NEAREST_POINTS_BINS_2D by allowing the ! user to specify the number of bins in each dimension. The main reason ! for doing this is to efficiently handle problems where the extent ! of the region varies widely from one dimension to another. ! ! A set of NSET points with coordinates PSET is assumed to lie within a ! rectangle. The limits of this rectangle are given in BIN_MIN and BIN_MAX. ! The rectangle is divided up into an NBIN(1) by NBIN(2) regular grid of ! cells. ! ! The cells are ordered lexicographically, as suggested by the following ! diagram for NBIN = (/ 5, 4 /) ! ! *----*----*----*----*----* ! | | | | | | ! | 16 | 17 | 18 | 19 | 20 | ! | | | | | | ! *----*----*----*----*----* ! | | | | | | ! | 11 | 12 | 13 | 14 | 15 | ! | | | | | | ! *----*----*----*----*----* ! | | | | | | ! | 6 | 7 | 8 | 9 | 10 | ! | | | | | | ! *----*----*----*----*----* ! | | | | | | ! | 1 | 2 | 3 | 4 | 5 | ! | | | | | | ! *----*----*----*----*----* ! ! The points in PSET are ordered by cell, and within a cell they ! are ordered lexicographically. Thus, for points P1 and P2, ! ! P1 < P2 implies that ! * P1 is in a lower ordered cell than P2, or ! * P1 is in the same cell as P2, but P1.X < P2.X, or ! * P1 is in the same cell as P2, P1.X = P2.X, but P1.Y < P2.Y. ! ! The arrays BIN_START(I,J) and BIN_LAST(I,J) are given the coordinates ! I, J of a cell, and return the lowest and highest index in PSET of a ! point in the I, J cell. All indices in between also belong to ! this cell. If the cell has no points, then both arrays are -1. ! ! ! After all this preprocessing, the algorithm for finding the nearest ! point goes as follows: ! ! 1) for a point PTEST, determine the cell it belongs to. ! Note that this algorithm will NOT be valid if PTEST lies outside ! the bin limits. ! ! 2) Search for a cell that has at least one member of PSET in it. ! We start at the cell containing PTEST, but if there are no members ! there, we spiral out from the cell, one layer at a time. ! ! 3) Within this cell, find the point nearest to PTEST. We now know that ! we don't need to search any cell whose points will all be further ! from PTEST than this distance. ! ! 4) Now, search in all other cells that could have a point closer to ! PTEST than what we have found so far. ! ! Modified: ! ! 18 March 2001 ! ! Author: ! ! John Burkardt ! ! Reference: ! ! Jon Bentley, Bruce Weide, Andrew Yao, ! Optimal Expected Time Algorithms for Closest Point Problems, ! ACM Transactions on Mathematical Software, ! Volume 6, Number 4, December 1980, pages 563-580. ! ! Parameters: ! ! Input, integer NSET, the number of points in the set. ! ! Input, real ( kind = 8 ) PSET(2,NSET), the points in the set. ! ! Input, integer NBIN(2), the number of cells in the horizontal and ! vertical directions. ! ! Input, real ( kind = 8 ) BIN_MIN(2), BIN_MAX(2), the minimum and ! maximum bin values. ! ! Input, integer BIN_START(NBIN(1),NBIN(2)), BIN_LAST(NBIN(1),NBIN(2)), ! indicates the index of the first and last element in the bin, ! or -1 if none. ! ! Input, integer BIN_NEXT(NSET), the index of the next element of the bin ! containing this element. ! ! Input, integer NTEST, the number of test points. ! ! Input, real ( kind = 8 ) PTEST(2,NTEST), the test points. ! ! Output, integer I_MIN(NTEST), the index of the nearest point in PSET ! to PTEST(ITEST). ! ! Output, real ( kind = 8 ) DIST_MIN(NTEST), the distance between ! PTEST(*,ITEST) and PSET(*,I_MIN). ! ! Output, integer COMPARES(NTEST), the number of point-to-point comparisons. ! implicit none integer, parameter :: ndim = 2 integer nbin(ndim) integer nset integer ntest integer bin(ndim) integer bin_last(nbin(1),nbin(2)) real ( kind = 8 ) bin_max(ndim) real ( kind = 8 ) bin_min(ndim) integer bin_start(nbin(1),nbin(2)) integer bin_next(nset) real ( kind = 8 ) c_max(ndim) real ( kind = 8 ) c_min(ndim) real ( kind = 8 ) cell_width_i integer compares(ntest) real ( kind = 8 ) dist_min(ntest) real ( kind = 8 ) dist_min_sq real ( kind = 8 ) disc real ( kind = 8 ) first_dist logical found_a_neighbor integer i integer i_min(ntest) integer itest integer k real ( kind = 8 ) pset(ndim,nset) real ( kind = 8 ) ptest(ndim,ntest) integer rank integer search1(ndim) integer search2(ndim) logical searched_everywhere integer searched_hi(ndim) integer searched_hi_new(ndim) integer searched_lo(ndim) integer searched_lo_new(ndim) real ( kind = 8 ) w real ( kind = 8 ) wall real ( kind = 8 ) wall_dist real ( kind = 8 ) width real ( kind = 8 ) width_i compares(1:ntest) = 0 ! ! Special cases. ! if ( nset <= 0 ) then dist_min(1:ntest) = huge ( dist_min(1) ) i_min(1:ntest) = 0 return end if if ( nset == 1 ) then do itest = 1, ntest dist_min(itest) = sqrt ( & sum ( ( ptest(1:ndim,itest) - pset(1:ndim,1) )**2 ) ) end do compares(1:ntest) = 1 i_min(1:ntest) = 1 return end if ! ! Search for each of the NTEST points. ! do itest = 1, ntest ! ! PART ONE: Initialize data. ! ! Determine the bin coordinates of the test point. ! Determine the limits of the bin containing P. ! Search the bin. ! Set the indices of the searched region to the index of this bin. ! dist_min_sq = huge ( dist_min_sq ) i_min(itest) = 0 found_a_neighbor = .false. searched_everywhere = .false. call r82_to_bin_even3 ( nbin, bin_min, bin_max, ptest(1,itest), bin ) call bin_to_r82_even3 ( nbin, bin, bin_min, bin_max, c_min, c_max ) call bin_search_one_2d ( bin, nset, pset, nbin, bin_start, bin_next, & ptest(1,itest), found_a_neighbor, i_min(itest), dist_min_sq, & compares(itest) ) if ( found_a_neighbor) then dist_min(itest) = sqrt ( dist_min_sq ) end if searched_lo(1:ndim) = bin(1:ndim) searched_hi(1:ndim) = bin(1:ndim) if ( all ( nbin(1:ndim) == 1 ) ) then if ( found_a_neighbor ) then dist_min(itest) = sqrt ( dist_min_sq ) else dist_min(itest) = huge ( dist_min(1) ) end if searched_everywhere = .true. cycle end if ! ! PART TWO. Look for a neighbor. ! ! Expand the search area in each dimension. ! Consider push the search limits out in every direction by 1 index. ! Determine the maximum width of the search region achieved in this way. ! Where possible, move all smaller indices out to the maximum width. ! ! Organize the search of the annexed region by dimension. ! Decrease and increase the first dimension to its new limits. ! Repeat for each dimension. ! ! Jump to next section ALMOST as soon as you have found a neighbor; ! just finish up the search in that coordinate direction, so that the ! search can be picked up cleanly in the final section. ! do while ( .not. found_a_neighbor ) do i = 1, ndim searched_hi_new(i) = min ( searched_hi(i) + 1, nbin(i) ) searched_lo_new(i) = max ( searched_lo(i) - 1, 1 ) end do width = 0.0D+00 do i = 1, ndim width_i = ( searched_hi_new(i) + 1 - searched_lo_new(i) ) * & ( bin_max(i) - bin_min(i) ) / nbin(i) width = max ( width, width_i ) end do do i = 1, ndim cell_width_i = ( bin_max(i) - bin_min(i) ) / nbin(i) width_i = ( searched_hi_new(i) + 1 - searched_lo_new(i) ) * cell_width_i disc = width - width_i if ( 2.0D+00 * cell_width_i <= disc ) then k = int ( disc / ( 2.0D+00 * cell_width_i ) ) searched_hi_new(i) = searched_hi_new(i) + k searched_lo_new(i) = searched_lo_new(i) - k searched_hi_new(i) = min ( searched_hi_new(i), nbin(i) ) searched_lo_new(i) = max ( searched_lo_new(i), 1 ) end if end do searched_everywhere = .true. do i = 1, ndim if ( searched_hi(i) < searched_hi_new(i) ) then searched_everywhere = .false. exit end if if ( searched_lo_new(i) < searched_lo(i) ) then searched_everywhere = .false. exit end if end do if ( searched_everywhere ) then exit end if do i = 1, ndim if ( searched_lo_new(i) < searched_lo(i) ) then search1(1:ndim) = searched_lo(1:ndim) search1(i) = searched_lo_new(i) search2(1:i-1) = searched_hi(1:i-1) search2(i) = searched_lo(i) - 1 search2(i+1:ndim) = searched_lo(i+1:ndim) rank = 0 do call tuple_next2 ( ndim, search1, search2, bin, rank ) if ( rank == 0 ) then exit end if call bin_search_one_2d ( bin, nset, pset, nbin, bin_start, & bin_next, ptest(1,itest), found_a_neighbor, i_min(itest), & dist_min_sq, compares(itest) ) end do searched_lo(i) = searched_lo_new(i) if ( found_a_neighbor ) then exit end if end if if ( searched_hi(i) < searched_hi_new(i) ) then search1(1:i-1) = searched_lo(1:ndim) search1(i) = searched_hi(i) + 1 search1(i+1:ndim) = searched_hi(1:ndim) search2(1:ndim) = searched_hi(1:ndim) search2(i) = searched_hi_new(i) rank = 0 do call tuple_next2 ( ndim, search1, search2, bin, rank ) if ( rank == 0 ) then exit end if call bin_search_one_2d ( bin, nset, pset, nbin, bin_start, & bin_next, ptest(1,itest), found_a_neighbor, i_min(itest), & dist_min_sq, compares(itest) ) end do searched_hi(i) = searched_hi_new(i) if ( found_a_neighbor ) then exit end if end if end do dist_min(itest) = sqrt ( dist_min_sq ) end do ! ! PART THREE: Final search ! ! You have found a neighbor in PSET to PTEST. ! If the neighbor is closer than the nearest wall which might have ! something on the other side, you're done. ! Otherwise, expand the region enough in each direction so that, once ! it is searched, we are sure to be done. ! wall_dist = huge ( wall_dist ) do i = 1, ndim if ( 1 < searched_lo(i) ) then wall = & ( real ( nbin(i) - searched_lo(i) + 1, kind = 8 ) * bin_min(i) & + real ( searched_lo(i) - 1, kind = 8 ) * bin_max(i) ) & / real ( nbin(i), kind = 8 ) wall_dist = min ( wall_dist, ptest(i,itest) - wall ) end if if ( searched_hi(i) < nbin(i) ) then wall = & ( real ( nbin(i) - searched_hi(i), kind = 8 ) * bin_min(i) & + real ( searched_hi(i), kind = 8 ) * bin_max(i) ) & / real ( nbin(i), kind = 8 ) wall_dist = min ( wall_dist, wall - ptest(i,itest) ) end if end do first_dist = dist_min(itest) if ( first_dist < wall_dist ) then cycle end if do i = 1, ndim if ( 1 < searched_lo(i) ) then ! ! Solve for SEARCH_NEW(I) so that FIRST_DIST < PTEST(I,ITEST) - WALL. ! w = ( ptest(i,itest) - first_dist - bin_min(i) ) & / real ( bin_max(i) - bin_min(i), kind = 8 ) k = int ( real ( nbin(i), kind = 8 ) * w ) k = max ( k, 1 ) search1(1:ndim) = searched_lo(1:ndim) search1(i) = k search2(1:ndim) = searched_hi(1:ndim) search2(i) = searched_lo(i) - 1 if ( search1(i) <= search2(i) ) then rank = 0 do call tuple_next2 ( ndim, search1, search2, bin, rank ) if ( rank == 0 ) then exit end if call bin_search_one_2d ( bin, nset, pset, nbin, bin_start, & bin_next, ptest(1,itest), found_a_neighbor, i_min(itest), & dist_min_sq, compares(itest) ) end do searched_lo(i) = k end if end if if ( searched_hi(i) < nbin(i) ) then ! ! Solve for SEARCH_NEW(I) so that FIRST_DIST < WALL - PTEST(I,ITEST). ! w = ( ptest(i,itest) + first_dist - bin_min(i) ) & / real ( bin_max(i) - bin_min(i), kind = 8 ) k = 1 + int ( real ( nbin(i), kind = 8 ) * w ) k = min ( k, nbin(i) ) search1(1:ndim) = searched_lo(1:ndim) search1(i) = searched_hi(i)+1 search2(1:ndim) = searched_hi(1:ndim) search2(i) = k if ( search1(i) <= search2(i) ) then rank = 0 do call tuple_next2 ( ndim, search1, search2, bin, rank ) if ( rank == 0 ) then exit end if call bin_search_one_2d ( bin, nset, pset, nbin, bin_start, & bin_next, ptest(1,itest), found_a_neighbor, i_min(itest), & dist_min_sq, compares(itest) ) end do searched_hi(i) = k end if end if end do dist_min(itest) = sqrt ( dist_min_sq ) end do return end subroutine points_nearest_points_bins_3d_2 ( nset, pset, nbin, bin_min, & bin_max, bin_start, bin_last, bin_next, ntest, ptest, i_min, dist_min, compares ) !*****************************************************************************80 ! !! POINTS_NEAREST_POINTS_BINS_3D_2 finds the nearest point to given points in 3D. ! ! Discussion: ! ! A set of NSET points with coordinates PSET is assumed to lie within a ! rectangle. The limits of this rectangle are given in BIN_MIN and BIN_MAX. ! The rectangle is divided up into an NBIN by NBIN regular grid of cells. ! ! The cells are ordered lexicographically, as suggested by the following ! diagram for NBIN = 5: ! ! *----*----*----*----*----* ! | | | | | | ! | 21 | 22 | 23 | 24 | 25 | ! | | | | | | ! *----*----*----*----*----* ! | | | | | | ! | 16 | 17 | 18 | 19 | 20 | ! | | | | | | ! *----*----*----*----*----* ! | | | | | | ! | 11 | 12 | 13 | 14 | 15 | ! | | | | | | ! *----*----*----*----*----* ! | | | | | | ! | 6 | 7 | 8 | 9 | 10 | ! | | | | | | ! *----*----*----*----*----* ! | | | | | | ! | 1 | 2 | 3 | 4 | 5 | ! | | | | | | ! *----*----*----*----*----* ! ! The points in PSET are ordered by cell, and within a cell they ! are ordered lexicographically. Thus, for points P1 and P2, ! ! P1 < P2 implies that ! * P1 is in a lower ordered cell than P2, or ! * P1 is in the same cell as P2, but P1.X < P2.X, or ! * P1 is in the same cell as P2, P1.X = P2.X, but P1.Y < P2.Y. ! ! The arrays BIN_START(I,J,K) and BIN_LAST(I,J,K) are given the coordinates ! I, J, K of a cell, and return the lowest and highest index in PSET of a ! point in the I, J, K cell. All indices in between also belong to ! this cell. If the cell has no points, then both arrays are -1. ! ! ! After all this preprocessing, the algorithm for finding the nearest ! point goes as follows: ! ! 1) for a point PTEST, determine the cell it belongs to. ! Note that this algorithm will NOT be valid if PTEST lies outside ! the bin limits. ! ! 2) Search for a cell that has at least one member of PSET in it. ! We start at the cell containing PTEST, but if there are no members ! there, we spiral out from the cell, one layer at a time. ! ! 3) Within this cell, find the point nearest to PTEST. We now know that ! we don't need to search any cell whose points will all be further ! from PTEST than this distance. ! ! 4) Now, search in all other cells that could have a point closer to ! PTEST than what we have found so far. ! ! Modified: ! ! 15 February 2001 ! ! Author: ! ! John Burkardt ! ! Reference: ! ! Jon Bentley, Bruce Weide, Andrew Yao, ! Optimal Expected Time Algorithms for Closest Point Problems, ! ACM Transactions on Mathematical Software, ! Volume 6, Number 4, December 1980, pages 563-580. ! ! Parameters: ! ! Input, integer NSET, the number of points in the set. ! ! Input, real ( kind = 8 ) PSET(3,NSET), the points in the set. ! ! Input, integer NBIN, the number of cells. NBIN must be at least 3. ! ! Input, real ( kind = 8 ) BIN_MIN(3), BIN_MAX(3), the minimum and ! maximum bin values. ! ! Input, integer BIN_START(NBIN,NBIN,NBIN), BIN_LAST(NBIN,NBIN,NBIN), ! the index of the first and last element in the bin, or -1 if none. ! ! Input, integer BIN_NEXT(NSET), the index of the next element of the bin ! containing this element. ! ! Input, integer NTEST, the number of test points. ! ! Input, real ( kind = 8 ) PTEST(3,NTEST), the test points. ! ! Output, integer I_MIN(NTEST), the index of the nearest point in PSET ! to PTEST(ITEST). ! ! Output, real ( kind = 8 ) DIST_MIN(NTEST), the distance between ! PTEST(*,ITEST) and PSET(*,I_MIN). ! ! Output, integer COMPARES(NTEST), the number of point-to-point comparisons. ! implicit none integer nbin integer, parameter :: ndim = 3 integer nset integer ntest integer bin(ndim) integer bin_last(nbin,nbin,nbin) real ( kind = 8 ) bin_max(ndim) real ( kind = 8 ) bin_min(ndim) integer bin_start(nbin,nbin,nbin) integer bin_next(nset) real ( kind = 8 ) c_max(ndim) real ( kind = 8 ) c_min(ndim) integer compares(ntest) real ( kind = 8 ) dist_min(ntest) real ( kind = 8 ) dist_min_sq real ( kind = 8 ) dist_sq integer i integer i_min(ntest) integer ic integer itest integer j integer jc integer k integer kc integer layer real ( kind = 8 ) layer_width logical more_bins integer node real ( kind = 8 ) pset(ndim,nset) real ( kind = 8 ) ptest(ndim,ntest) real ( kind = 8 ) search_radius compares(1:ntest) = 0 ! ! Special cases. ! if ( nset <= 0 ) then dist_min(1:ntest) = huge ( dist_min ) i_min(1:ntest) = 0 return end if if ( nset == 1 ) then do itest = 1, ntest dist_min(itest) = sqrt ( & sum ( ( ptest(1:ndim,itest) - pset(1:ndim,1) )**2 ) ) end do compares(1:ntest) = 1 i_min(1:ntest) = 1 return end if layer_width = minval ( bin_max(1:ndim) - bin_min(1:ndim) ) & / real ( nbin, kind = 8 ) ! ! Search for each of the NTEST points. ! do itest = 1, ntest dist_min_sq = huge ( dist_min_sq ) i_min(itest) = 0 search_radius = 0.0D+00 ! ! Determine the bin coordinates of the point P. ! call r83_to_bin_even2 ( nbin, bin_min, bin_max, ptest(1,itest), bin ) ! ! Determine the radius of the ball of space that will be completely ! searched after this first bin, "layer 0", is done. ! call bin_to_r83_even2 ( nbin, bin, bin_min, bin_max, c_min, c_max ) ! ! Set the central bin of the layers. ! ic = bin(1) jc = bin(2) kc = bin(3) layer = 0 ! ! Search all legal bins in layer LAYER. ! do more_bins = .false. call index_box2_next_3d ( layer, layer, layer, ic, jc, kc, i, j, k, & more_bins ) ! ! In layer LAYER, search each BIN I, J, K. ! do if ( 1 <= i .and. i <= nbin .and. & 1 <= j .and. j <= nbin .and. & 1 <= k .and. k <= nbin ) then node = bin_start(i,j,k) do while ( 0 < node ) dist_sq = sum ( ( ptest(1:ndim,itest) - pset(1:ndim,node) )**2 ) compares(itest) = compares(itest) + 1 if ( dist_sq < dist_min_sq ) then dist_min_sq = dist_sq i_min(itest) = node end if node = bin_next(node) end do end if ! ! We have now searched all points in bin I, J, K. ! ! Determine the next bin in this layer. ! ! But if it lies outside the region, discard it, and go on to the next one. ! Once you get past the last bin, exit because you're done the layer. ! do call index_box2_next_3d ( layer, layer, layer, ic, jc, kc, & i, j, k, more_bins ) if ( .not. more_bins ) then exit end if if ( 1 <= i .and. i <= nbin .and. & 1 <= j .and. j <= nbin .and. & 1 <= k .and. k <= nbin ) then exit end if end do if ( .not. more_bins ) then exit end if end do ! ! We've completed layer LAYER. ! Update the radius of the searched area. ! if ( layer == 0 ) then search_radius = min ( & minval ( abs ( ptest(1:ndim,itest) - c_min(1:ndim) ) ), & minval ( abs ( ptest(1:ndim,itest) - c_max(1:ndim) ) ) ) else search_radius = search_radius + layer_width end if ! ! We are done with PTEST(*,ITEST) if: ! ! * We've found at least one neighbor; ! AND ! * We've searched all legal boxes that contain the circle ! with PTEST at the center and the nearest neighbor on the circumference. ! if ( i_min(itest) /= 0 ) then dist_min(itest) = sqrt ( dist_min_sq ) if ( dist_min(itest) <= search_radius ) then exit end if end if ! ! Prepare to search the next layer. ! layer = layer + 1 end do end do return end subroutine points_nearest_points_bins_3d_3 ( nset, pset, nbin, bin_min, & bin_max, bin_start, bin_last, bin_next, ntest, ptest, i_min, dist_min, compares ) !*****************************************************************************80 ! !! POINTS_NEAREST_POINTS_BINS_3D_3 finds the nearest point to given points in 3D. ! ! Discussion: ! ! This code differs from POINTS_NEAREST_POINTS_BINS_3D by allowing the ! user to specify the number of bins in each dimension. ! ! A set of NSET points with coordinates PSET is assumed to lie within a ! box. The limits of this rectangle are given in BIN_MIN and BIN_MAX. ! The rectangle is divided up into an NBIN(1) by NBIN(2) by NBIN(3) ! regular grid of cells. ! ! The cells are ordered lexicographically, as suggested by the following ! diagram for NBIN = (/ 5, 4, 2 /) ! ! Z LAYER 1 Z LAYER 2 ! *----*----*----*----*----* *----*----*----*----*----* ! | | | | | | | | | | | | ! | 16 | 17 | 18 | 19 | 20 | | 36 | 37 | 38 | 39 | 40 | ! | | | | | | | | | | | | ! *----*----*----*----*----* *----*----*----*----*----* ! | | | | | | | | | | | | ! | 11 | 12 | 13 | 14 | 15 | | 31 | 32 | 33 | 34 | 35 | ! | | | | | | | | | | | | ! *----*----*----*----*----* *----*----*----*----*----* ! | | | | | | | | | | | | ! | 6 | 7 | 8 | 9 | 10 | | 26 | 27 | 28 | 29 | 30 | ! | | | | | | | | | | | | ! *----*----*----*----*----* *----*----*----*----*----* ! | | | | | | | | | | | | ! | 1 | 2 | 3 | 4 | 5 | | 21 | 22 | 23 | 24 | 25 | ! | | | | | | | | | | | | ! *----*----*----*----*----* *----*----*----*----*----* ! ! The points in PSET are ordered by cell, and within a cell they ! are ordered lexicographically. Thus, for points P1 and P2, ! ! P1 < P2 implies that ! * P1 is in a lower ordered cell than P2, or ! * P1 is in the same cell as P2, but P1.X < P2.X, or ! * P1 is in the same cell as P2, P1.X = P2.X, but P1.Y < P2.Y. ! * P1 is in the same cell as P2, P1.X = P2.X, P1.Y = P2.Y, but P1.Z < P2.Z ! ! The arrays BIN_START(I,J,K) and BIN_LAST(I,J,K) are given the coordinates ! I, J, K of a cell, and return the lowest and highest index in PSET of a ! point in the I, J, K cell. All indices in between also belong to ! this cell. If the cell has no points, then both arrays are -1. ! ! ! After all this preprocessing, the algorithm for finding the nearest ! point goes as follows: ! ! 1) for a point PTEST, determine the cell it belongs to. ! Note that this algorithm will NOT be valid if PTEST lies outside ! the bin limits. ! ! 2) Search for a cell that has at least one member of PSET in it. ! We start at the cell containing PTEST, but if there are no members ! there, we spiral out from the cell, one layer at a time. ! ! 3) Within this cell, find the point nearest to PTEST. We now know that ! we don't need to search any cell whose points will all be f