11 # include "geompack.H"
13 //******************************************************************************
15 double d_epsilon ( void )
17 //******************************************************************************
21 // D_EPSILON returns the round off unit for double precision arithmetic.
25 // D_EPSILON is a number R which is a power of 2 with the property that,
26 // to the precision of the computer's arithmetic,
41 // Output, double D_EPSILON, the floating point round-off unit.
55 //*********************************************************************
57 double d_max ( double x, double y )
59 //*********************************************************************
63 // D_MAX returns the maximum of two real values.
75 // Input, double X, Y, the quantities to compare.
77 // Output, double D_MAX, the maximum of X and Y.
89 //*********************************************************************
91 double d_min ( double x, double y )
93 //*********************************************************************
97 // D_MIN returns the minimum of two real values.
109 // Input, double X, Y, the quantities to compare.
111 // Output, double D_MIN, the minimum of X and Y.
123 //******************************************************************************
125 void d2vec_part_quick_a ( int n, double a[], int *l, int *r )
127 //******************************************************************************
131 // D2VEC_PART_QUICK_A reorders an R2 vector as part of a quick sort.
135 // The routine reorders the entries of A. Using A(1:2,1) as a
136 // key, all entries of A that are less than or equal to the key will
137 // precede the key, which precedes all entries that are greater than the key.
145 // A = ( (2,4), (8,8), (6,2), (0,2), (10,6), (10,0), (0,6), (4,8) )
151 // A = ( (0,2), (0,6), (2,4), (8,8), (6,2), (10,6), (10,0), (4,8) )
152 // ----------- ----------------------------------
165 // Input, int N, the number of entries of A.
167 // Input/output, double A[N*2]. On input, the array to be checked.
168 // On output, A has been reordered as described above.
170 // Output, int *L, *R, the indices of A that define the three segments.
171 // Let KEY = the input value of A(1:2,1). Then
172 // I <= L A(1:2,I) < KEY;
173 // L < I < R A(1:2,I) = KEY;
174 // R <= I A(1:2,I) > KEY.
187 cout << "D2VEC_PART_QUICK_A - Fatal error!\n";
203 // The elements of unknown size have indices between L+1 and R-1.
208 for ( i = 2; i <= n; i++ )
210 if ( dvec_gt ( 2, a+2*ll, key ) )
213 dvec_swap ( 2, a+2*(rr-1), a+2*ll );
215 else if ( dvec_eq ( 2, a+2*ll, key ) )
218 dvec_swap ( 2, a+2*(m-1), a+2*ll );
221 else if ( dvec_lt ( 2, a+2*ll, key ) )
228 // Now shift small elements to the left, and KEY elements to center.
230 for ( i = 0; i < ll - m; i++ )
232 for ( j = 0; j < 2; j++ )
234 a[2*i+j] = a[2*(i+m)+j];
240 for ( i = ll; i < ll+m; i++ )
242 for ( j = 0; j < 2; j++ )
253 //*******************************************************************************
255 void d2vec_permute ( int n, double a[], int p[] )
257 //*******************************************************************************
261 // D2VEC_PERMUTE permutes an R2 vector in place.
265 // This routine permutes an array of real "objects", but the same
266 // logic can be used to permute an array of objects of any arithmetic
267 // type, or an array of objects of any complexity. The only temporary
268 // storage required is enough to store a single object. The number
269 // of data movements made is N + the number of cycles of order 2 or more,
270 // which is never more than N + N/2.
277 // P = ( 2, 4, 5, 1, 3 )
278 // A = ( 1.0, 2.0, 3.0, 4.0, 5.0 )
279 // (11.0, 22.0, 33.0, 44.0, 55.0 )
283 // A = ( 2.0, 4.0, 5.0, 1.0, 3.0 )
284 // ( 22.0, 44.0, 55.0, 11.0, 33.0 ).
296 // Input, int N, the number of objects.
298 // Input/output, double A[2*N], the array to be permuted.
300 // Input, int P[N], the permutation. P(I) = J means
301 // that the I-th element of the output array should be the J-th
302 // element of the input array. P must be a legal permutation
303 // of the integers from 1 to N, otherwise the algorithm will
304 // fail catastrophically.
313 if ( !perm_check ( n, p ) )
316 cout << "D2VEC_PERMUTE - Fatal error!\n";
317 cout << " The input array does not represent\n";
318 cout << " a proper permutation.\n";
322 // Search for the next element of the permutation that has not been used.
324 for ( istart = 1; istart <= n; istart++ )
326 if ( p[istart-1] < 0 )
330 else if ( p[istart-1] == istart )
332 p[istart-1] = -p[istart-1];
337 a_temp[0] = a[0+(istart-1)*2];
338 a_temp[1] = a[1+(istart-1)*2];
341 // Copy the new value into the vacated entry.
348 p[iput-1] = -p[iput-1];
350 if ( iget < 1 || n < iget )
353 cout << "D2VEC_PERMUTE - Fatal error!\n";
357 if ( iget == istart )
359 a[0+(iput-1)*2] = a_temp[0];
360 a[1+(iput-1)*2] = a_temp[1];
363 a[0+(iput-1)*2] = a[0+(iget-1)*2];
364 a[1+(iput-1)*2] = a[1+(iget-1)*2];
369 // Restore the signs of the entries.
371 for ( i = 0; i < n; i++ )
378 //******************************************************************************
380 int *d2vec_sort_heap_index_a ( int n, double a[] )
382 //******************************************************************************
386 // D2VEC_SORT_HEAP_INDEX_A does an indexed heap ascending sort of an R2 vector.
390 // The sorting is not actually carried out. Rather an index array is
391 // created which defines the sorting. This array may be used to sort
392 // or index the array, or to sort or index related arrays keyed on the
395 // Once the index array is computed, the sorting can be carried out
398 // A(1:2,INDX(I)), I = 1 to N is sorted,
400 // or explicitly, by the call
402 // call D2VEC_PERMUTE ( N, A, INDX )
404 // after which A(1:2,I), I = 1 to N is sorted.
416 // Input, int N, the number of entries in the array.
418 // Input, double A[2*N], an array to be index-sorted.
420 // Output, int D2VEC_SORT_HEAP_INDEX_A[N], the sort index. The
421 // I-th element of the sorted array is A(0:1,D2VEC_SORT_HEAP_INDEX_A(I-1)).
444 indx = ivec_indicator ( n );
455 aval[0] = a[0+(indxt-1)*2];
456 aval[1] = a[1+(indxt-1)*2];
461 aval[0] = a[0+(indxt-1)*2];
462 aval[1] = a[1+(indxt-1)*2];
463 indx[ir-1] = indx[0];
481 if ( a[0+(indx[j-1]-1)*2] < a[0+(indx[j]-1)*2] ||
482 ( a[0+(indx[j-1]-1)*2] == a[0+(indx[j]-1)*2] &&
483 a[1+(indx[j-1]-1)*2] < a[1+(indx[j]-1)*2] ) )
489 if ( aval[0] < a[0+(indx[j-1]-1)*2] ||
490 ( aval[0] == a[0+(indx[j-1]-1)*2] &&
491 aval[1] < a[1+(indx[j-1]-1)*2] ) )
493 indx[i-1] = indx[j-1];
507 //*****************************************************************************
509 void d2vec_sort_quick_a ( int n, double a[] )
511 //*****************************************************************************
515 // D2VEC_SORT_QUICK_A ascending sorts an R2 vector using quick sort.
519 // The data structure is a set of N pairs of real numbers.
520 // These values are stored in a one dimensional array, by pairs.
532 // Input, int N, the number of entries in the array.
534 // Input/output, double A[N*2].
535 // On input, the array to be sorted.
536 // On output, the array has been sorted.
539 # define LEVEL_MAX 25
545 int rsave[LEVEL_MAX];
551 cout << "D2VEC_SORT_QUICK_A - Fatal error!\n";
562 rsave[level-1] = n + 1;
566 while ( 0 < n_segment )
569 // Partition the segment.
571 d2vec_part_quick_a ( n_segment, a+2*(base-1)+0, &l_segment, &r_segment );
573 // If the left segment has more than one element, we need to partition it.
577 if ( LEVEL_MAX < level )
580 cout << "D2VEC_SORT_QUICK_A - Fatal error!\n";
581 cout << " Exceeding recursion maximum of " << LEVEL_MAX << "\n";
586 n_segment = l_segment;
587 rsave[level-1] = r_segment + base - 1;
590 // The left segment and the middle segment are sorted.
591 // Must the right segment be partitioned?
593 else if ( r_segment < n_segment )
595 n_segment = n_segment + 1 - r_segment;
596 base = base + r_segment - 1;
599 // Otherwise, we back up a level if there is an earlier one.
611 base = rsave[level-1];
612 n_segment = rsave[level-2] - rsave[level-1];
625 //******************************************************************************
627 int diaedg ( double x0, double y0, double x1, double y1, double x2, double y2,
628 double x3, double y3 )
630 //******************************************************************************
634 // DIAEDG chooses a diagonal edge.
638 // The routine determines whether 0--2 or 1--3 is the diagonal edge
639 // that should be chosen, based on the circumcircle criterion, where
640 // (X0,Y0), (X1,Y1), (X2,Y2), (X3,Y3) are the vertices of a simple
641 // quadrilateral in counterclockwise order.
650 // Department of Computing Science,
651 // University of Alberta,
652 // Edmonton, Alberta, Canada T6G 2H1
657 // GEOMPACK - a software package for the generation of meshes
658 // using geometric algorithms,
659 // Advances in Engineering Software,
660 // Volume 13, pages 325-331, 1991.
664 // Input, double X0, Y0, X1, Y1, X2, Y2, X3, Y3, the coordinates of the
665 // vertices of a quadrilateral, given in counter clockwise order.
667 // Output, int DIAEDG, chooses a diagonal:
668 // +1, if diagonal edge 02 is chosen;
669 // -1, if diagonal edge 13 is chosen;
670 // 0, if the four vertices are cocircular.
689 tol = 100.0 * d_epsilon ( );
700 tola = tol * d_max ( fabs ( dx10 ),
701 d_max ( fabs ( dy10 ),
702 d_max ( fabs ( dx30 ), fabs ( dy30 ) ) ) );
704 tolb = tol * d_max ( fabs ( dx12 ),
705 d_max ( fabs ( dy12 ),
706 d_max ( fabs ( dx32 ), fabs ( dy32 ) ) ) );
708 ca = dx10 * dx30 + dy10 * dy30;
709 cb = dx12 * dx32 + dy12 * dy32;
711 if ( tola < ca && tolb < cb )
715 else if ( ca < -tola && cb < -tolb )
721 tola = d_max ( tola, tolb );
722 s = ( dx10 * dy30 - dx30 * dy10 ) * cb
723 + ( dx32 * dy12 - dx12 * dy32 ) * ca;
729 else if ( s < -tola )
742 //******************************************************************************
744 void dmat_transpose_print ( int m, int n, double a[], const char *title )
746 //******************************************************************************
750 // DMAT_TRANSPOSE_PRINT prints a real matrix, transposed.
762 // Input, int M, N, the number of rows and columns.
764 // Input, double A[M*N], an M by N matrix to be printed.
766 // Input, const char *TITLE, an optional title.
769 dmat_transpose_print_some ( m, n, a, 1, 1, m, n, title );
773 //******************************************************************************
775 void dmat_transpose_print_some ( int m, int n, double a[], int ilo, int jlo,
776 int ihi, int jhi, const char *title )
778 //******************************************************************************
782 // DMAT_TRANSPOSE_PRINT_SOME prints some of a real matrix, transposed.
794 // Input, int M, N, the number of rows and columns.
796 // Input, double A[M*N], an M by N matrix to be printed.
798 // Input, int ILO, JLO, the first row and column to print.
800 // Input, int IHI, JHI, the last row and column to print.
802 // Input, const char *TITLE, an optional title.
816 if ( 0 < s_len_trim ( title ) )
819 cout << title << "\n";
822 for ( i2lo = i_max ( ilo, 1 ); i2lo <= i_min ( ihi, m ); i2lo = i2lo + INCX )
824 i2hi = i2lo + INCX - 1;
825 i2hi = i_min ( i2hi, m );
826 i2hi = i_min ( i2hi, ihi );
828 inc = i2hi + 1 - i2lo;
832 for ( i = i2lo; i <= i2hi; i++ )
834 cout << setw(7) << i << " ";
840 j2lo = i_max ( jlo, 1 );
841 j2hi = i_min ( jhi, n );
843 for ( j = j2lo; j <= j2hi; j++ )
845 cout << setw(5) << j << " ";
846 for ( i2 = 1; i2 <= inc; i2++ )
849 cout << setw(14) << a[(i-1)+(j-1)*m];
859 //******************************************************************************
861 void dmat_uniform ( int m, int n, double b, double c, int *seed, double r[] )
863 //******************************************************************************
867 // DMAT_UNIFORM fills a double precision array with scaled pseudorandom values.
871 // This routine implements the recursion
873 // seed = 16807 * seed mod ( 2**31 - 1 )
874 // unif = seed / ( 2**31 - 1 )
876 // The integer arithmetic never requires more than 32 bits,
877 // including a sign bit.
889 // Paul Bratley, Bennett Fox, Linus Schrage,
890 // A Guide to Simulation,
891 // Springer Verlag, pages 201-202, 1983.
895 // Implementation and Relative Efficiency of Quasirandom
896 // Sequence Generators,
897 // ACM Transactions on Mathematical Software,
898 // Volume 12, Number 4, pages 362-376, 1986.
900 // Peter Lewis, Allen Goodman, James Miller,
901 // A Pseudo-Random Number Generator for the System/360,
902 // IBM Systems Journal,
903 // Volume 8, pages 136-143, 1969.
907 // Input, int M, N, the number of rows and columns.
909 // Input, double B, C, the limits of the pseudorandom values.
911 // Input/output, int *SEED, the "seed" value. Normally, this
912 // value should not be 0, otherwise the output value of SEED
913 // will still be 0, and D_UNIFORM will be 0. On output, SEED has
916 // Output, double DMAT_UNIFORM[M*N], a matrix of pseudorandom values.
923 for ( j = 0; j < n; j++ )
925 for ( i = 0; i < m; i++ )
929 *seed = 16807 * ( *seed - k * 127773 ) - k * 2836;
933 *seed = *seed + 2147483647;
936 // Although SEED can be represented exactly as a 32 bit integer,
937 // it generally cannot be represented exactly as a 32 bit real number!
939 r[i+j*m] = b + ( c - b ) * double( *seed ) * 4.656612875E-10;
945 //******************************************************************************
947 int dtris2 ( int point_num, double point_xy[], int *tri_num,
948 int tri_vert[], int tri_nabe[] )
950 //******************************************************************************
954 // DTRIS2 constructs a Delaunay triangulation of 2D vertices.
958 // The routine constructs the Delaunay triangulation of a set of 2D vertices
959 // using an incremental approach and diagonal edge swaps. Vertices are
960 // first sorted in lexicographically increasing (X,Y) order, and
961 // then are inserted one at a time from outside the convex hull.
970 // Department of Computing Science,
971 // University of Alberta,
972 // Edmonton, Alberta, Canada T6G 2H1
977 // GEOMPACK - a software package for the generation of meshes
978 // using geometric algorithms,
979 // Advances in Engineering Software,
980 // Volume 13, pages 325-331, 1991.
984 // Input, int POINT_NUM, the number of vertices.
986 // Input/output, double POINT_XY[POINT_NUM*2], the coordinates of the vertices.
987 // On output, the vertices have been sorted into dictionary order.
989 // Output, int *TRI_NUM, the number of triangles in the triangulation;
990 // TRI_NUM is equal to 2*POINT_NUM - NB - 2, where NB is the number
991 // of boundary vertices.
993 // Output, int TRI_VERT[TRI_NUM*3], the nodes that make up each triangle.
994 // The elements are indices of POINT_XY. The vertices of the triangles are
995 // in counter clockwise order.
997 // Output, int TRI_NABE[TRI_NUM*3], the triangle neighbor list.
998 // Positive elements are indices of TIL; negative elements are used for links
999 // of a counter clockwise linked list of boundary edges; LINK = -(3*I + J-1)
1000 // where I, J = triangle, edge index; TRI_NABE[I,J] refers to
1001 // the neighbor along edge from vertex J to J+1 (mod 3).
1003 // Output, int DTRIS2, is 0 for no error.
1027 stack = new int[point_num];
1029 tol = 100.0 * d_epsilon ( );
1031 // Sort the vertices by increasing (x,y).
1033 indx = d2vec_sort_heap_index_a ( point_num, point_xy );
1035 d2vec_permute ( point_num, point_xy, indx );
1037 // Make sure that the data points are "reasonably" distinct.
1041 for ( i = 2; i <= point_num; i++ )
1048 for ( j = 0; j <= 1; j++ )
1050 cmax = d_max ( fabs ( point_xy[2*(m-1)+j] ),
1051 fabs ( point_xy[2*(m1-1)+j] ) );
1053 if ( tol * ( cmax + 1.0 )
1054 < fabs ( point_xy[2*(m-1)+j] - point_xy[2*(m1-1)+j] ) )
1065 cout << "DTRIS2 - Fatal error!\n";
1066 cout << " Fails for point number I = " << i << "\n";
1067 cout << " M = " << m << "\n";
1068 cout << " M1 = " << m1 << "\n";
1069 cout << " X,Y(M) = " << point_xy[2*(m-1)+0] << " "
1070 << point_xy[2*(m-1)+1] << "\n";
1071 cout << " X,Y(M1) = " << point_xy[2*(m1-1)+0] << " "
1072 << point_xy[2*(m1-1)+1] << "\n";
1079 // Starting from points M1 and M2, search for a third point M that
1080 // makes a "healthy" triangle (M1,M2,M)
1088 if ( point_num < j )
1091 cout << "DTRIS2 - Fatal error!\n";
1098 lr = lrline ( point_xy[2*(m-1)+0], point_xy[2*(m-1)+1],
1099 point_xy[2*(m1-1)+0], point_xy[2*(m1-1)+1],
1100 point_xy[2*(m2-1)+0], point_xy[2*(m2-1)+1], 0.0 );
1111 // Set up the triangle information for (M1,M2,M), and for any other
1112 // triangles you created because points were collinear with M1, M2.
1118 tri_vert[3*0+0] = m1;
1119 tri_vert[3*0+1] = m2;
1120 tri_vert[3*0+2] = m;
1121 tri_nabe[3*0+2] = -3;
1123 for ( i = 2; i <= *tri_num; i++ )
1127 tri_vert[3*(i-1)+0] = m1;
1128 tri_vert[3*(i-1)+1] = m2;
1129 tri_vert[3*(i-1)+2] = m;
1130 tri_nabe[3*(i-1)+0] = -3 * i;
1131 tri_nabe[3*(i-1)+1] = i;
1132 tri_nabe[3*(i-1)+2] = i - 1;
1136 tri_nabe[3*(*tri_num-1)+0] = -3 * (*tri_num) - 1;
1137 tri_nabe[3*(*tri_num-1)+1] = -5;
1143 tri_vert[3*0+0] = m2;
1144 tri_vert[3*0+1] = m1;
1145 tri_vert[3*0+2] = m;
1146 tri_nabe[3*0+0] = -4;
1148 for ( i = 2; i <= *tri_num; i++ )
1152 tri_vert[3*(i-1)+0] = m2;
1153 tri_vert[3*(i-1)+1] = m1;
1154 tri_vert[3*(i-1)+2] = m;
1155 tri_nabe[3*(i-2)+2] = i;
1156 tri_nabe[3*(i-1)+0] = -3 * i - 3;
1157 tri_nabe[3*(i-1)+1] = i - 1;
1160 tri_nabe[3*(*tri_num-1)+2] = -3 * (*tri_num);
1161 tri_nabe[3*0+1] = -3 * (*tri_num) - 2;
1166 // Insert the vertices one at a time from outside the convex hull,
1167 // determine visible boundary edges, and apply diagonal edge swaps until
1168 // Delaunay triangulation of vertices (so far) is obtained.
1172 for ( i = j+1; i <= point_num; i++ )
1175 m1 = tri_vert[3*(ltri-1)+ledg-1];
1179 m2 = tri_vert[3*(ltri-1)+ledg];
1183 m2 = tri_vert[3*(ltri-1)+0];
1186 lr = lrline ( point_xy[2*(m-1)+0], point_xy[2*(m-1)+1],
1187 point_xy[2*(m1-1)+0], point_xy[2*(m1-1)+1],
1188 point_xy[2*(m2-1)+0], point_xy[2*(m2-1)+1], 0.0 );
1198 l = -tri_nabe[3*(ltri-1)+ledg-1];
1203 vbedg ( point_xy[2*(m-1)+0], point_xy[2*(m-1)+1], point_num,
1204 point_xy, *tri_num, tri_vert, tri_nabe, <ri, &ledg, &rtri, &redg );
1207 l = -tri_nabe[3*(ltri-1)+ledg-1];
1213 l = -tri_nabe[3*(t-1)+e-1];
1214 m2 = tri_vert[3*(t-1)+e-1];
1218 m1 = tri_vert[3*(t-1)+e];
1222 m1 = tri_vert[3*(t-1)+0];
1225 *tri_num = *tri_num + 1;
1226 tri_nabe[3*(t-1)+e-1] = *tri_num;
1227 tri_vert[3*(*tri_num-1)+0] = m1;
1228 tri_vert[3*(*tri_num-1)+1] = m2;
1229 tri_vert[3*(*tri_num-1)+2] = m;
1230 tri_nabe[3*(*tri_num-1)+0] = t;
1231 tri_nabe[3*(*tri_num-1)+1] = *tri_num - 1;
1232 tri_nabe[3*(*tri_num-1)+2] = *tri_num + 1;
1235 if ( point_num < top )
1238 cout << "DTRIS2 - Fatal error!\n";
1239 cout << " Stack overflow.\n";
1244 stack[top-1] = *tri_num;
1246 if ( t == rtri && e == redg )
1253 tri_nabe[3*(ltri-1)+ledg-1] = -3 * n - 1;
1254 tri_nabe[3*(n-1)+1] = -3 * (*tri_num) - 2;
1255 tri_nabe[3*(*tri_num-1)+2] = -l;
1259 error = swapec ( m, &top, <ri, &ledg, point_num, point_xy, *tri_num,
1260 tri_vert, tri_nabe, stack );
1265 cout << "DTRIS2 - Fatal error!\n";
1266 cout << " Error return from SWAPEC.\n";
1273 // Now account for the sorting that we did.
1275 for ( i = 0; i < 3; i++ )
1277 for ( j = 0; j < *tri_num; j++ )
1279 tri_vert[i+j*3] = indx [ tri_vert[i+j*3] - 1 ];
1283 perm_inv ( point_num, indx );
1285 d2vec_permute ( point_num, point_xy, indx );
1292 //******************************************************************************
1294 bool dvec_eq ( int n, double a1[], double a2[] )
1296 //******************************************************************************
1300 // DVEC_EQ is true if every pair of entries in two vectors is equal.
1312 // Input, int N, the number of entries in the vectors.
1314 // Input, double A1[N], A2[N], two vectors to compare.
1316 // Output, bool DVEC_EQ.
1317 // DVEC_EQ is TRUE if every pair of elements A1(I) and A2(I) are equal,
1318 // and FALSE otherwise.
1323 for ( i = 0; i < n; i++ )
1325 if ( a1[i] != a2[i] )
1333 //******************************************************************************
1335 bool dvec_gt ( int n, double a1[], double a2[] )
1337 //******************************************************************************
1341 // DVEC_GT == ( A1 > A2 ) for real vectors.
1345 // The comparison is lexicographic.
1347 // A1 > A2 <=> A1(1) > A2(1) or
1348 // ( A1(1) == A2(1) and A1(2) > A2(2) ) or
1350 // ( A1(1:N-1) == A2(1:N-1) and A1(N) > A2(N)
1362 // Input, int N, the dimension of the vectors.
1364 // Input, double A1[N], A2[N], the vectors to be compared.
1366 // Output, bool DVEC_GT, is TRUE if and only if A1 > A2.
1371 for ( i = 0; i < n; i++ )
1374 if ( a2[i] < a1[i] )
1378 else if ( a1[i] < a2[i] )
1387 //******************************************************************************
1389 bool dvec_lt ( int n, double a1[], double a2[] )
1391 //******************************************************************************
1395 // DVEC_LT == ( A1 < A2 ) for real vectors.
1399 // The comparison is lexicographic.
1401 // A1 < A2 <=> A1(1) < A2(1) or
1402 // ( A1(1) == A2(1) and A1(2) < A2(2) ) or
1404 // ( A1(1:N-1) == A2(1:N-1) and A1(N) < A2(N)
1416 // Input, int N, the dimension of the vectors.
1418 // Input, double A1[N], A2[N], the vectors to be compared.
1420 // Output, bool DVEC_LT, is TRUE if and only if A1 < A2.
1425 for ( i = 0; i < n; i++ )
1427 if ( a1[i] < a2[i] )
1431 else if ( a2[i] < a1[i] )
1440 //********************************************************************
1442 void dvec_print ( int n, double a[], const char *title )
1444 //********************************************************************
1448 // DVEC_PRINT prints a double precision vector.
1460 // Input, int N, the number of components of the vector.
1462 // Input, double A[N], the vector to be printed.
1464 // Input, const char *TITLE, a title to be printed first.
1465 // TITLE may be blank.
1470 if ( 0 < s_len_trim ( title ) )
1473 cout << title << "\n";
1477 for ( i = 0; i <= n-1; i++ )
1479 cout << setw(6) << i + 1 << " "
1480 << setw(14) << a[i] << "\n";
1485 //******************************************************************************
1487 void dvec_swap ( int n, double a1[], double a2[] )
1489 //******************************************************************************
1493 // DVEC_SWAP swaps the entries of two real vectors.
1505 // Input, int N, the number of entries in the arrays.
1507 // Input/output, double A1[N], A2[N], the vectors to swap.
1513 for ( i = 0; i < n; i++ )
1522 //****************************************************************************
1524 int i_max ( int i1, int i2 )
1526 //****************************************************************************
1530 // I_MAX returns the maximum of two integers.
1542 // Input, int I1, I2, are two integers to be compared.
1544 // Output, int I_MAX, the larger of I1 and I2.
1557 //****************************************************************************
1559 int i_min ( int i1, int i2 )
1561 //****************************************************************************
1565 // I_MIN returns the smaller of two integers.
1577 // Input, int I1, I2, two integers to be compared.
1579 // Output, int I_MIN, the smaller of I1 and I2.
1592 //*********************************************************************
1594 int i_modp ( int i, int j )
1596 //*********************************************************************
1600 // I_MODP returns the nonnegative remainder of integer division.
1605 // NREM = I_MODP ( I, J )
1606 // NMULT = ( I - NREM ) / J
1608 // I = J * NMULT + NREM
1609 // where NREM is always nonnegative.
1613 // The MOD function computes a result with the same sign as the
1614 // quantity being divided. Thus, suppose you had an angle A,
1615 // and you wanted to ensure that it was between 0 and 360.
1616 // Then mod(A,360) would do, if A was positive, but if A
1617 // was negative, your result would be between -360 and 0.
1619 // On the other hand, I_MODP(A,360) is between 0 and 360, always.
1623 // I J MOD I_MODP I_MODP Factorization
1625 // 107 50 7 7 107 = 2 * 50 + 7
1626 // 107 -50 7 7 107 = -2 * -50 + 7
1627 // -107 50 -7 43 -107 = -3 * 50 + 43
1628 // -107 -50 -7 43 -107 = 3 * -50 + 43
1640 // Input, int I, the number to be divided.
1642 // Input, int J, the number that divides I.
1644 // Output, int I_MODP, the nonnegative remainder when I is
1653 cout << "I_MODP - Fatal error!\n";
1654 cout << " I_MODP ( I, J ) called with J = " << j << "\n";
1662 value = value + abs ( j );
1667 //********************************************************************
1669 int i_sign ( int i )
1671 //********************************************************************
1675 // I_SIGN returns the sign of an integer.
1679 // The sign of 0 and all positive integers is taken to be +1.
1680 // The sign of all negative integers is -1.
1692 // Input, int I, the integer whose sign is desired.
1694 // Output, int I_SIGN, the sign of I.
1706 //*******************************************************************************
1708 int i_wrap ( int ival, int ilo, int ihi )
1710 //*******************************************************************************
1714 // I_WRAP forces an integer to lie between given limits by wrapping.
1750 // Input, int IVAL, an integer value.
1752 // Input, int ILO, IHI, the desired bounds for the integer value.
1754 // Output, int I_WRAP, a "wrapped" version of IVAL.
1762 jlo = i_min ( ilo, ihi );
1763 jhi = i_max ( ilo, ihi );
1765 wide = jhi + 1 - jlo;
1773 value = jlo + i_modp ( ival - jlo, wide );
1778 //******************************************************************************
1780 void imat_transpose_print ( int m, int n, int a[], const char *title )
1782 //******************************************************************************
1786 // IMAT_TRANSPOSE_PRINT prints an integer matrix, transposed.
1798 // Input, int M, the number of rows in A.
1800 // Input, int N, the number of columns in A.
1802 // Input, int A[M*N], the M by N matrix.
1804 // Input, const char *TITLE, a title to be printed.
1807 imat_transpose_print_some ( m, n, a, 1, 1, m, n, title );
1811 //******************************************************************************
1813 void imat_transpose_print_some ( int m, int n, int a[], int ilo, int jlo,
1814 int ihi, int jhi, const char *title )
1816 //******************************************************************************
1820 // IMAT_TRANSPOSE_PRINT_SOME prints some of an integer matrix, transposed.
1832 // Input, int M, the number of rows of the matrix.
1833 // M must be positive.
1835 // Input, int N, the number of columns of the matrix.
1836 // N must be positive.
1838 // Input, int A[M*N], the matrix.
1840 // Input, int ILO, JLO, IHI, JHI, designate the first row and
1841 // column, and the last row and column to be printed.
1843 // Input, const char *TITLE, a title for the matrix.
1854 if ( 0 < s_len_trim ( title ) )
1857 cout << title << "\n";
1860 // Print the columns of the matrix, in strips of INCX.
1862 for ( i2lo = ilo; i2lo <= ihi; i2lo = i2lo + INCX )
1864 i2hi = i2lo + INCX - 1;
1865 i2hi = i_min ( i2hi, m );
1866 i2hi = i_min ( i2hi, ihi );
1870 // For each row I in the current range...
1872 // Write the header.
1875 for ( i = i2lo; i <= i2hi; i++ )
1877 cout << setw(7) << i << " ";
1883 // Determine the range of the rows in this strip.
1885 j2lo = i_max ( jlo, 1 );
1886 j2hi = i_min ( jhi, n );
1888 for ( j = j2lo; j <= j2hi; j++ )
1891 // Print out (up to INCX) entries in column J, that lie in the current strip.
1893 cout << setw(5) << j << " ";
1894 for ( i = i2lo; i <= i2hi; i++ )
1896 cout << setw(6) << a[i-1+(j-1)*m] << " ";
1908 //********************************************************************
1910 void ivec_heap_d ( int n, int a[] )
1912 /*********************************************************************
1916 // IVEC_HEAP_D reorders an array of integers into a descending heap.
1920 // A heap is an array A with the property that, for every index J,
1921 // A[J] >= A[2*J+1] and A[J] >= A[2*J+2], (as long as the indices
1922 // 2*J+1 and 2*J+2 are legal).
1930 // A(3) A(4) A(5) A(6)
1932 // A(7) A(8) A(9) A(10)
1936 // Albert Nijenhuis, Herbert Wilf,
1937 // Combinatorial Algorithms,
1938 // Academic Press, 1978, second edition,
1939 // ISBN 0-12-519260-6.
1951 // Input, int N, the size of the input array.
1953 // Input/output, int A[N].
1954 // On input, an unsorted array.
1955 // On output, the array has been reordered into a heap.
1963 // Only nodes (N/2)-1 down to 0 can be "parent" nodes.
1965 for ( i = (n/2)-1; 0 <= i; i-- )
1968 // Copy the value out of the parent node.
1969 // Position IFREE is now "open".
1977 // Positions 2*IFREE + 1 and 2*IFREE + 2 are the descendants of position
1978 // IFREE. (One or both may not exist because they equal or exceed N.)
1982 // Does the first position exist?
1991 // Does the second position exist?
1996 // If both positions exist, take the larger of the two values,
1997 // and update M if necessary.
1999 if ( a[m] < a[m+1] )
2005 // If the large descendant is larger than KEY, move it up,
2006 // and update IFREE, the location of the free position, and
2007 // consider the descendants of THIS position.
2023 // When you have stopped shifting items up, return the item you
2024 // pulled out back to the heap.
2032 //******************************************************************************
2034 int *ivec_indicator ( int n )
2036 //******************************************************************************
2040 // IVEC_INDICATOR sets an integer vector to the indicator vector.
2052 // Input, int N, the number of elements of A.
2054 // Output, int IVEC_INDICATOR(N), the initialized array.
2062 for ( i = 0; i < n; i++ )
2069 //********************************************************************
2071 void ivec_sort_heap_a ( int n, int a[] )
2073 //********************************************************************
2077 // IVEC_SORT_HEAP_A ascending sorts an array of integers using heap sort.
2081 // Albert Nijenhuis, Herbert Wilf,
2082 // Combinatorial Algorithms,
2083 // Academic Press, 1978, second edition,
2084 // ISBN 0-12-519260-6.
2096 // Input, int N, the number of entries in the array.
2098 // Input/output, int A[N].
2099 // On input, the array to be sorted;
2100 // On output, the array has been sorted.
2111 // 1: Put A into descending heap form.
2113 ivec_heap_d ( n, a );
2117 // The largest object in the heap is in A[0].
2118 // Move it to position A[N-1].
2124 // Consider the diminished heap of size N1.
2126 for ( n1 = n-1; 2 <= n1; n1-- )
2129 // Restore the heap structure of the initial N1 entries of A.
2131 ivec_heap_d ( n1, a );
2133 // Take the largest object from A[0] and move it to A[N1-1].
2143 //******************************************************************************
2145 void ivec_sorted_unique ( int n, int a[], int *nuniq )
2147 //******************************************************************************
2151 // IVEC_SORTED_UNIQUE finds unique elements in a sorted integer array.
2155 // 02 September 2003
2163 // Input, int N, the number of elements in A.
2165 // Input/output, int A[N]. On input, the sorted
2166 // integer array. On output, the unique elements in A.
2168 // Output, int *NUNIQ, the number of unique elements in A.
2182 for ( i = 1; i < n; i++ )
2184 if ( a[i] != a[*nuniq] )
2186 *nuniq = *nuniq + 1;
2194 //******************************************************************************
2196 int lrline ( double xu, double yu, double xv1, double yv1, double xv2,
2197 double yv2, double dv )
2199 //******************************************************************************
2203 // LRLINE determines where a point lies in relation to a directed line.
2207 // LRLINE determines whether a point is to the left of, right of,
2208 // or on a directed line parallel to a line through given points.
2217 // Department of Computing Science,
2218 // University of Alberta,
2219 // Edmonton, Alberta, Canada T6G 2H1
2224 // GEOMPACK - a software package for the generation of meshes
2225 // using geometric algorithms,
2226 // Advances in Engineering Software,
2227 // Volume 13, pages 325-331, 1991.
2231 // Input, double XU, YU, XV1, YV1, XV2, YV2, are vertex coordinates; the
2232 // directed line is parallel to and at signed distance DV to the left of
2233 // the directed line from (XV1,YV1) to (XV2,YV2); (XU,YU) is the vertex for
2234 // which the position relative to the directed line is to be determined.
2236 // Input, double DV, the signed distance, positive for left.
2238 // Output, int LRLINE, is +1, 0, or -1 depending on whether (XU,YU) is
2239 // to the right of, on, or left of the directed line. LRLINE is 0 if
2240 // the line degenerates to a point.
2248 double tol = 0.0000001;
2257 tolabs = tol * d_max ( fabs ( dx ),
2258 d_max ( fabs ( dy ),
2259 d_max ( fabs ( dxu ),
2260 d_max ( fabs ( dyu ), fabs ( dv ) ) ) ) );
2262 t = dy * dxu - dx * dyu + dv * sqrt ( dx * dx + dy * dy );
2268 else if ( -tolabs <= t )
2272 else if ( t < -tolabs )
2279 //******************************************************************************
2281 bool perm_check ( int n, int p[] )
2283 //******************************************************************************
2287 // PERM_CHECK checks that a vector represents a permutation.
2291 // The routine verifies that each of the integers from 1
2292 // to N occurs among the N entries of the permutation.
2304 // Input, int N, the number of entries.
2306 // Input, int P[N], the array to check.
2308 // Output, bool PERM_CHECK, is TRUE if the permutation is OK.
2315 for ( seek = 1; seek <= n; seek++ )
2319 for ( i = 0; i < n; i++ )
2337 //******************************************************************************
2339 void perm_inv ( int n, int p[] )
2341 //******************************************************************************
2345 // PERM_INV inverts a permutation "in place".
2353 // Input, int N, the number of objects being permuted.
2355 // Input/output, int P[N], the permutation, in standard index form.
2356 // On output, P describes the inverse permutation
2368 cout << "PERM_INV - Fatal error!\n";
2369 cout << " Input value of N = " << n << "\n";
2373 if ( !perm_check ( n, p ) )
2376 cout << "PERM_INV - Fatal error!\n";
2377 cout << " The input array does not represent\n";
2378 cout << " a proper permutation.\n";
2384 for ( i = 1; i <= n; i++ )
2395 is = - i_sign ( p[i-1] );
2396 p[i-1] = i_sign ( is ) * abs ( p[i-1] );
2399 for ( i = 1; i <= n; i++ )
2424 //********************************************************************
2426 int *points_delaunay_naive_2d ( int n, double p[], int *ntri )
2428 //********************************************************************
2432 // POINTS_DELAUNAY_NAIVE_2D computes the Delaunay triangulation in 2D.
2436 // A naive and inefficient (but extremely simple) method is used.
2438 // This routine is only suitable as a demonstration code for small
2439 // problems. Its running time is of order N^4. Much faster algorithms
2442 // Given a set of nodes in the plane, a triangulation is a set of
2443 // triples of distinct nodes, forming triangles, so that every
2444 // point with the convex hull of the set of nodes is either one
2445 // of the nodes, or lies on an edge of one or more triangles,
2446 // or lies within exactly one triangle.
2459 // Computational Geometry,
2460 // Cambridge University Press,
2461 // Second Edition, 1998, page 187.
2465 // Input, int N, the number of nodes. N must be at least 3.
2467 // Input, double P[2*N], the coordinates of the nodes.
2469 // Output, int *NTRI, the number of triangles.
2471 // Output, int POINTS_DELAUNAY_NAIVE_2D[3*NTRI], the indices of the
2472 // nodes making each triangle.
2490 z = new double [ n ];
2492 for ( i = 0; i < n; i++ )
2494 z[i] = p[0+i*2] * p[0+i*2] + p[1+i*2] * p[1+i*2];
2497 // First pass counts triangles,
2498 // Second pass allocates triangles and sets them.
2500 for ( pass = 1; pass <= 2; pass++ )
2504 tri = new int[3*count];
2508 // For each triple (I,J,K):
2510 for ( i = 0; i < n - 2; i++ )
2512 for ( j = i+1; j < n; j++ )
2514 for ( k = i+1; k < n; k++ )
2518 xn = ( p[1+j*2] - p[1+i*2] ) * ( z[k] - z[i] )
2519 - ( p[1+k*2] - p[1+i*2] ) * ( z[j] - z[i] );
2520 yn = ( p[0+k*2] - p[0+i*2] ) * ( z[j] - z[i] )
2521 - ( p[0+j*2] - p[0+i*2] ) * ( z[k] - z[i] );
2522 zn = ( p[0+j*2] - p[0+i*2] ) * ( p[1+k*2] - p[1+i*2] )
2523 - ( p[0+k*2] - p[0+i*2] ) * ( p[1+j*2] - p[1+i*2] );
2529 for ( m = 0; m < n; m++ )
2531 flag = flag && ( ( p[0+m*2] - p[0+i*2] ) * xn
2532 + ( p[1+m*2] - p[1+i*2] ) * yn
2533 + ( z[m] - z[i] ) * zn <= 0 );
2559 //******************************************************************************
2561 int s_len_trim ( const char *s )
2563 //******************************************************************************
2567 // S_LEN_TRIM returns the length of a string to the last nonblank.
2579 // Input, const char *S, a pointer to a string.
2581 // Output, int S_LEN_TRIM, the length of the string to the last nonblank.
2582 // If S_LEN_TRIM is 0, then the string is entirely blank.
2589 t = const_cast<char*>(s) + n - 1;
2603 //******************************************************************************
2605 int swapec ( int i, int *top, int *btri, int *bedg, int point_num,
2606 double point_xy[], int tri_num, int tri_vert[], int tri_nabe[],
2609 //******************************************************************************
2613 // SWAPEC swaps diagonal edges until all triangles are Delaunay.
2617 // The routine swaps diagonal edges in a 2D triangulation, based on
2618 // the empty circumcircle criterion, until all triangles are Delaunay,
2619 // given that I is the index of the new vertex added to the triangulation.
2623 // 03 September 2003
2628 // Department of Computing Science,
2629 // University of Alberta,
2630 // Edmonton, Alberta, Canada T6G 2H1
2635 // GEOMPACK - a software package for the generation of meshes
2636 // using geometric algorithms,
2637 // Advances in Engineering Software,
2638 // Volume 13, pages 325-331, 1991.
2642 // Input, int I, the index of the new vertex.
2644 // Input/output, int *TOP, the index of the top of the stack.
2645 // On output, TOP is zero.
2647 // Input/output, int *BTRI, *BEDG; on input, if positive, are the
2648 // triangle and edge indices of a boundary edge whose updated indices
2649 // must be recorded. On output, these may be updated because of swaps.
2651 // Input, int POINT_NUM, the number of points.
2653 // Input, double POINT_XY[POINT_NUM*2], the coordinates of the points.
2655 // Input, int TRI_NUM, the number of triangles.
2657 // Input/output, int TRI_VERT[TRI_NUM*3], the triangle incidence list.
2658 // May be updated on output because of swaps.
2660 // Input/output, int TRI_NABE[TRI_NUM*3], the triangle neighbor list;
2661 // negative values are used for links of the counter-clockwise linked
2662 // list of boundary edges; May be updated on output because of swaps.
2664 // LINK = -(3*I + J-1) where I, J = triangle, edge index.
2666 // Workspace, int STACK[MAXST]; on input, entries 1 through TOP
2667 // contain the indices of initial triangles (involving vertex I)
2668 // put in stack; the edges opposite I should be in interior; entries
2669 // TOP+1 through MAXST are used as a stack.
2671 // Output, int SWAPEC, is set to 8 for abnormal return.
2694 // Determine whether triangles in stack are Delaunay, and swap
2695 // diagonal edge of convex quadrilateral if not.
2697 x = point_xy[2*(i-1)+0];
2698 y = point_xy[2*(i-1)+1];
2707 t = stack[(*top)-1];
2710 if ( tri_vert[3*(t-1)+0] == i )
2713 b = tri_vert[3*(t-1)+2];
2715 else if ( tri_vert[3*(t-1)+1] == i )
2718 b = tri_vert[3*(t-1)+0];
2723 b = tri_vert[3*(t-1)+1];
2726 a = tri_vert[3*(t-1)+e-1];
2727 u = tri_nabe[3*(t-1)+e-1];
2729 if ( tri_nabe[3*(u-1)+0] == t )
2732 c = tri_vert[3*(u-1)+2];
2734 else if ( tri_nabe[3*(u-1)+1] == t )
2737 c = tri_vert[3*(u-1)+0];
2742 c = tri_vert[3*(u-1)+1];
2745 swap = diaedg ( x, y,
2746 point_xy[2*(a-1)+0], point_xy[2*(a-1)+1],
2747 point_xy[2*(c-1)+0], point_xy[2*(c-1)+1],
2748 point_xy[2*(b-1)+0], point_xy[2*(b-1)+1] );
2752 em1 = i_wrap ( e - 1, 1, 3 );
2753 ep1 = i_wrap ( e + 1, 1, 3 );
2754 fm1 = i_wrap ( f - 1, 1, 3 );
2755 fp1 = i_wrap ( f + 1, 1, 3 );
2757 tri_vert[3*(t-1)+ep1-1] = c;
2758 tri_vert[3*(u-1)+fp1-1] = i;
2759 r = tri_nabe[3*(t-1)+ep1-1];
2760 s = tri_nabe[3*(u-1)+fp1-1];
2761 tri_nabe[3*(t-1)+ep1-1] = u;
2762 tri_nabe[3*(u-1)+fp1-1] = t;
2763 tri_nabe[3*(t-1)+e-1] = s;
2764 tri_nabe[3*(u-1)+f-1] = r;
2766 if ( 0 < tri_nabe[3*(u-1)+fm1-1] )
2769 stack[(*top)-1] = u;
2774 if ( tri_nabe[3*(s-1)+0] == u )
2776 tri_nabe[3*(s-1)+0] = t;
2778 else if ( tri_nabe[3*(s-1)+1] == u )
2780 tri_nabe[3*(s-1)+1] = t;
2784 tri_nabe[3*(s-1)+2] = t;
2789 if ( point_num < *top )
2794 stack[(*top)-1] = t;
2798 if ( u == *btri && fp1 == *bedg )
2804 l = - ( 3 * t + e - 1 );
2808 while ( 0 < tri_nabe[3*(tt-1)+ee-1] )
2810 tt = tri_nabe[3*(tt-1)+ee-1];
2812 if ( tri_vert[3*(tt-1)+0] == a )
2816 else if ( tri_vert[3*(tt-1)+1] == a )
2827 tri_nabe[3*(tt-1)+ee-1] = l;
2833 if ( tri_nabe[3*(r-1)+0] == t )
2835 tri_nabe[3*(r-1)+0] = u;
2837 else if ( tri_nabe[3*(r-1)+1] == t )
2839 tri_nabe[3*(r-1)+1] = u;
2843 tri_nabe[3*(r-1)+2] = u;
2848 if ( t == *btri && ep1 == *bedg )
2854 l = - ( 3 * u + f - 1 );
2858 while ( 0 < tri_nabe[3*(tt-1)+ee-1] )
2860 tt = tri_nabe[3*(tt-1)+ee-1];
2862 if ( tri_vert[3*(tt-1)+0] == b )
2866 else if ( tri_vert[3*(tt-1)+1] == b )
2875 tri_nabe[3*(tt-1)+ee-1] = l;
2881 //**********************************************************************
2883 void timestamp ( void )
2885 //**********************************************************************
2889 // TIMESTAMP prints the current YMDHMS date as a time stamp.
2893 // May 31 2001 09:45:54 AM
2908 # define TIME_SIZE 29
2910 static char time_buffer[TIME_SIZE];
2911 const struct tm *tm;
2915 now = time ( NULL );
2916 tm = localtime ( &now );
2918 len = strftime ( time_buffer, TIME_SIZE, "%d %B %Y %I:%M:%S %p", tm );
2922 cout << time_buffer << "\n";
2928 //**********************************************************************
2930 char *timestring ( void )
2932 //**********************************************************************
2936 // TIMESTRING returns the current YMDHMS date as a string.
2940 // May 31 2001 09:45:54 AM
2952 // Output, char *TIMESTRING, a string containing the current YMDHMS date.
2955 # define TIME_SIZE 29
2957 const struct tm *tm;
2962 now = time ( NULL );
2963 tm = localtime ( &now );
2965 s = new char[TIME_SIZE];
2967 len = strftime ( s, TIME_SIZE, "%d %B %Y %I:%M:%S %p", tm );
2972 //******************************************************************************
2974 double *triangle_circumcenter_2d ( double t[] )
2976 //******************************************************************************
2980 // TRIANGLE_CIRCUMCENTER_2D computes the circumcenter of a triangle in 2D.
2984 // The circumcenter of a triangle is the center of the circumcircle, the
2985 // circle that passes through the three vertices of the triangle.
2987 // The circumcircle contains the triangle, but it is not necessarily the
2988 // smallest triangle to do so.
2990 // If all angles of the triangle are no greater than 90 degrees, then
2991 // the center of the circumscribed circle will lie inside the triangle.
2992 // Otherwise, the center will lie outside the circle.
2994 // The circumcenter is the intersection of the perpendicular bisectors
2995 // of the sides of the triangle.
2997 // In geometry, the circumcenter of a triangle is often symbolized by "O".
3009 // Input, double T[2*3], the triangle vertices.
3011 // Output, double *X, *Y, the coordinates of the circumcenter of the triangle.
3023 center = new double[DIM_NUM];
3025 asq = ( t[0+1*2] - t[0+0*2] ) * ( t[0+1*2] - t[0+0*2] )
3026 + ( t[1+1*2] - t[1+0*2] ) * ( t[1+1*2] - t[1+0*2] );
3028 csq = ( t[0+2*2] - t[0+0*2] ) * ( t[0+2*2] - t[0+0*2] )
3029 + ( t[1+2*2] - t[1+0*2] ) * ( t[1+2*2] - t[1+0*2] );
3031 top1 = ( t[1+1*2] - t[1+0*2] ) * csq - ( t[1+2*2] - t[1+0*2] ) * asq;
3032 top2 = ( t[0+1*2] - t[0+0*2] ) * csq - ( t[0+2*2] - t[0+0*2] ) * asq;
3034 bot = ( t[1+1*2] - t[1+0*2] ) * ( t[0+2*2] - t[0+0*2] )
3035 - ( t[1+2*2] - t[1+0*2] ) * ( t[0+1*2] - t[0+0*2] );
3037 center[0] = t[0+0*2] + 0.5 * top1 / bot;
3038 center[1] = t[1+0*2] + 0.5 * top2 / bot;
3044 //******************************************************************************
3046 bool triangulation_plot_eps ( const char *file_out_name, int g_num, double g_xy[],
3047 int tri_num, int nod_tri[] )
3049 //******************************************************************************
3053 // TRIANGULATION_PLOT_EPS plots a triangulation of a pointset.
3057 // The triangulation is most usually a Delaunay triangulation,
3058 // but this is not necessary.
3060 // The data can be generated by calling DTRIS2, but this is not
3065 // 08 September 2003
3073 // Input, const char *FILE_OUT_NAME, the name of the output file.
3075 // Input, int G_NUM, the number of points.
3077 // Input, double G_XY[G_NUM,2], the coordinates of the points.
3079 // Input, int TRI_NUM, the number of triangles.
3081 // Input, int NOD_TRI[3,TRI_NUM], lists, for each triangle,
3082 // the indices of the points that form the vertices of the triangle.
3084 // Output, bool TRIANGULATION_PLOT_EPS, is TRUE for success.
3098 int x_ps_max_clip = 594;
3100 int x_ps_min_clip = 18;
3105 int y_ps_max_clip = 684;
3107 int y_ps_min_clip = 108;
3109 date_time = timestring ( );
3111 x_max = g_xy[0+0*2];
3112 x_min = g_xy[0+0*2];
3113 y_max = g_xy[1+0*2];
3114 y_min = g_xy[1+0*2];
3116 for ( g = 0; g < g_num; g++ )
3118 x_max = d_max ( x_max, g_xy[0+g*2] );
3119 x_min = d_min ( x_min, g_xy[0+g*2] );
3120 y_max = d_max ( y_max, g_xy[1+g*2] );
3121 y_min = d_min ( y_min, g_xy[1+g*2] );
3124 // Plot the Delaunay triangulation.
3127 // Open the output file.
3129 file_out.open ( file_out_name );
3134 cout << "TRIANGULATION_PLOT_EPS - Fatal error!\n";
3135 cout << " Cannot open the output file \"" << file_out_name << "\".\n";
3139 file_out << "%!PS-Adobe-3.0 EPSF-3.0\n";
3140 file_out << "%%Creator: triangulation_plot_eps.cc\n";
3141 file_out << "%%Title: " << file_out_name << "\n";
3142 file_out << "%%CreationDate: " << date_time << "\n";
3143 file_out << "%%Pages: 1\n";
3144 file_out << "%%Bounding Box: " << x_ps_min << " " << y_ps_min << " "
3145 << x_ps_max << " " << y_ps_max << "\n";
3146 file_out << "%%Document-Fonts: Times-Roman\n";
3147 file_out << "%%LanguageLevel: 1\n";
3148 file_out << "%%EndComments\n";
3149 file_out << "%%BeginProlog\n";
3150 file_out << "/inch {72 mul} def\n";
3151 file_out << "%%EndProlog\n";
3152 file_out << "%%Page: 1 1\n";
3153 file_out << "save\n";
3155 file_out << "% Set the RGB line color to very light gray.\n";
3157 file_out << "0.900 0.900 0.900 setrgbcolor\n";
3159 file_out << "% Draw a gray border around the page.\n";
3161 file_out << "newpath\n";
3162 file_out << " " << x_ps_min << " " << y_ps_min << " moveto\n";
3163 file_out << " " << x_ps_max << " " << y_ps_min << " lineto\n";
3164 file_out << " " << x_ps_max << " " << y_ps_max << " lineto\n";
3165 file_out << " " << x_ps_min << " " << y_ps_max << " lineto\n";
3166 file_out << " " << x_ps_min << " " << y_ps_min << " lineto\n";
3167 file_out << "stroke\n";
3169 file_out << "% Set the RGB line color to black.\n";
3171 file_out << "0.000 0.000 0.000 setrgbcolor\n";
3173 file_out << "% Set the font and its size.\n";
3175 file_out << "/Times-Roman findfont\n";
3176 file_out << "0.50 inch scalefont\n";
3177 file_out << "setfont\n";
3179 file_out << "% Print a title.\n";
3181 file_out << "210 702 moveto\n";
3182 file_out << "(Triangulation) show\n";
3184 file_out << "% Define a clipping polygon.\n";
3186 file_out << "newpath\n";
3187 file_out << " " << x_ps_min_clip << " " << y_ps_min_clip << " moveto\n";
3188 file_out << " " << x_ps_max_clip << " " << y_ps_min_clip << " lineto\n";
3189 file_out << " " << x_ps_max_clip << " " << y_ps_max_clip << " lineto\n";
3190 file_out << " " << x_ps_min_clip << " " << y_ps_max_clip << " lineto\n";
3191 file_out << " " << x_ps_min_clip << " " << y_ps_min_clip << " lineto\n";
3192 file_out << "clip newpath\n";
3194 file_out << "% Set the RGB line color to green.\n";
3196 file_out << "0.000 0.750 0.150 setrgbcolor\n";
3198 file_out << "% Draw the nodes.\n";
3201 for ( g = 0; g < g_num; g++ )
3204 ( ( x_max - g_xy[0+g*2] ) * double( x_ps_min )
3205 + ( g_xy[0+g*2] - x_min ) * double( x_ps_max ) )
3206 / ( x_max - x_min ) );
3209 ( ( y_max - g_xy[1+g*2] ) * double( y_ps_min )
3210 + ( g_xy[1+g*2] - y_min ) * double( y_ps_max ) )
3211 / ( y_max - y_min ) );
3213 file_out << "newpath " << x_ps << " "
3214 << y_ps << " 5 0 360 arc closepath fill\n";
3218 file_out << "% Set the RGB line color to red.\n";
3220 file_out << "0.900 0.200 0.100 setrgbcolor\n";
3222 file_out << "% Draw the triangles.\n";
3225 for ( t = 1; t <= tri_num; t++ )
3227 file_out << "newpath\n";
3229 for ( j = 1; j <= 4; j++ )
3231 e = i_wrap ( j, 1, 3 );
3233 k = nod_tri[3*(t-1)+e-1];
3236 ( ( x_max - g_xy[0+(k-1)*2] ) * double( x_ps_min )
3237 + ( g_xy[0+(k-1)*2] - x_min ) * double( x_ps_max ) )
3238 / ( x_max - x_min ) );
3241 ( ( y_max - g_xy[1+(k-1)*2] ) * double( y_ps_min )
3242 + ( g_xy[1+(k-1)*2] - y_min ) * double( y_ps_max ) )
3243 / ( y_max - y_min ) );
3247 file_out << x_ps << " " << y_ps << " moveto\n";
3251 file_out << x_ps << " " << y_ps << " lineto\n";
3256 file_out << "stroke\n";
3260 file_out << "restore showpage\n";
3262 file_out << "% End of page.\n";
3264 file_out << "%%Trailer\n";
3265 file_out << "%%EOF\n";
3271 //******************************************************************************
3273 void triangulation_print ( int point_num, double xc[], int tri_num,
3274 int tri_vert[], int tri_nabe[] )
3276 //******************************************************************************
3280 // TRIANGULATION_PRINT prints information defining a Delaunay triangulation.
3284 // Triangulations created by RTRIS include extra information encoded
3285 // in the negative values of TRI_NABE.
3287 // Because some of the nodes counted in POINT_NUM may not actually be
3288 // used in the triangulation, I needed to compute the true number
3289 // of vertices. I added this calculation on 13 October 2001.
3291 // Ernest Fasse pointed out an error in the indexing of VERTEX_LIST,
3292 // which was corrected on 19 February 2004.
3304 // Input, int POINT_NUM, the number of points.
3306 // Input, double XC[2*POINT_NUM], the point coordinates.
3308 // Input, int TRI_NUM, the number of triangles.
3310 // Input, int TRI_VERT[3*TRI_NUM], the nodes that make up the triangles.
3312 // Input, int TRI_NABE[3*TRI_NUM], the triangle neighbors on each side.
3313 // If there is no triangle neighbor on a particular side, the value of
3314 // TRI_NABE should be negative. If the triangulation data was created by
3315 // DTRIS2, then there is more information encoded in the negative values.
3335 cout << "TRIANGULATION_PRINT\n";
3336 cout << " Information defining a triangulation.\n";
3338 cout << " The number of points is " << point_num << "\n";
3340 dmat_transpose_print ( DIM_NUM, point_num, xc, " Point coordinates" );
3343 cout << " The number of triangles is " << tri_num << "\n";
3345 cout << " Sets of three points are used as vertices of\n";
3346 cout << " the triangles. For each triangle, the points\n";
3347 cout << " are listed in counterclockwise order.\n";
3349 imat_transpose_print ( 3, tri_num, tri_vert, " Triangle nodes" );
3352 cout << " On each side of a given triangle, there is either\n";
3353 cout << " another triangle, or a piece of the convex hull.\n";
3354 cout << " For each triangle, we list the indices of the three\n";
3355 cout << " neighbors, or (if negative) the codes of the\n";
3356 cout << " segments of the convex hull.\n";
3358 imat_transpose_print ( 3, tri_num, tri_nabe, " Triangle neighbors" );
3360 // Determine VERTEX_NUM, the number of vertices. This is not
3361 // the same as the number of points!
3363 vertex_list = new int[3*tri_num];
3366 for ( t = 0; t < tri_num; t++ )
3368 for ( s = 0; s < 3; s++ )
3370 vertex_list[k] = tri_vert[s+t*3];
3375 ivec_sort_heap_a ( 3*tri_num, vertex_list );
3377 ivec_sorted_unique ( 3*tri_num, vertex_list, &vertex_num );
3379 delete [] vertex_list;
3381 // Determine the number of boundary points.
3383 boundary_num = 2 * vertex_num - tri_num - 2;
3386 cout << " The number of boundary points is " << boundary_num << "\n";
3388 cout << " The segments that make up the convex hull can be\n";
3389 cout << " determined from the negative entries of the triangle\n";
3390 cout << " neighbor list.\n";
3392 cout << " # Tri Side N1 N2\n";
3399 for ( i = 0; i < tri_num; i++ )
3401 for ( j = 0; j < 3; j++ )
3403 if ( tri_nabe[j+i*3] < 0 )
3405 s = -tri_nabe[j+i*3];
3408 if ( t < 1 || tri_num < t )
3411 cout << " Sorry, this data does not use the DTRIS2\n";
3412 cout << " convention for convex hull segments.\n";
3418 s2 = i_wrap ( s1+1, 1, 3 );
3420 n1 = tri_vert[s1-1+(t-1)*3];
3421 n2 = tri_vert[s2-1+(t-1)*3];
3422 cout << setw(4) << k << " "
3423 << setw(4) << t << " "
3424 << setw(4) << s1 << " "
3425 << setw(4) << n1 << " "
3426 << setw(4) << n2 << "\n";
3441 //******************************************************************************
3443 void vbedg ( double x, double y, int point_num, double point_xy[], int tri_num,
3444 int tri_vert[], int tri_nabe[], int *ltri, int *ledg, int *rtri, int *redg )
3446 //******************************************************************************
3450 // VBEDG determines which boundary edges are visible to a point.
3454 // The point (X,Y) is assumed to be outside the convex hull of the
3455 // region covered by the 2D triangulation.
3460 // Department of Computing Science,
3461 // University of Alberta,
3462 // Edmonton, Alberta, Canada T6G 2H1
3467 // GEOMPACK - a software package for the generation of meshes
3468 // using geometric algorithms,
3469 // Advances in Engineering Software,
3470 // Volume 13, pages 325-331, 1991.
3474 // 02 September 2003
3478 // Input, double X, Y, the coordinates of a point outside the convex hull
3479 // of the current triangulation.
3481 // Input, int POINT_NUM, the number of points.
3483 // Input, double POINT_XY[POINT_NUM*2], the coordinates of the vertices.
3485 // Input, int TRI_NUM, the number of triangles.
3487 // Input, int TRI_VERT[TRI_NUM*3], the triangle incidence list.
3489 // Input, int TRI_NABE[TRI_NUM*3], the triangle neighbor list; negative
3490 // values are used for links of a counter clockwise linked list of boundary
3492 // LINK = -(3*I + J-1) where I, J = triangle, edge index.
3494 // Input/output, int *LTRI, *LEDG. If LTRI != 0 then these values are
3495 // assumed to be already computed and are not changed, else they are updated.
3496 // On output, LTRI is the index of boundary triangle to the left of the
3497 // leftmost boundary triangle visible from (X,Y), and LEDG is the boundary
3498 // edge of triangle LTRI to the left of the leftmost boundary edge visible
3499 // from (X,Y). 1 <= LEDG <= 3.
3501 // Input/output, int *RTRI. On input, the index of the boundary triangle
3502 // to begin the search at. On output, the index of the rightmost boundary
3503 // triangle visible from (X,Y).
3505 // Input/output, int *REDG, the edge of triangle RTRI that is visible
3506 // from (X,Y). 1 <= REDG <= 3.
3521 // Find the rightmost visible boundary edge using links, then possibly
3522 // leftmost visible boundary edge using triangle neighbor information.
3537 l = -tri_nabe[3*((*rtri)-1)+(*redg)-1];
3540 a = tri_vert[3*(t-1)+e-1];
3544 b = tri_vert[3*(t-1)+e];
3548 b = tri_vert[3*(t-1)+0];
3551 ax = point_xy[2*(a-1)+0];
3552 ay = point_xy[2*(a-1)+1];
3554 bx = point_xy[2*(b-1)+0];
3555 by = point_xy[2*(b-1)+1];
3557 lr = lrline ( x, y, ax, ay, bx, by, 0.0 );
3579 b = tri_vert[3*(t-1)+e-1];
3580 e = i_wrap ( e-1, 1, 3 );
3582 while ( 0 < tri_nabe[3*(t-1)+e-1] )
3584 t = tri_nabe[3*(t-1)+e-1];
3586 if ( tri_vert[3*(t-1)+0] == b )
3590 else if ( tri_vert[3*(t-1)+1] == b )
3601 a = tri_vert[3*(t-1)+e-1];
3602 ax = point_xy[2*(a-1)+0];
3603 ay = point_xy[2*(a-1)+1];
3605 bx = point_xy[2*(b-1)+0];
3606 by = point_xy[2*(b-1)+1];
3608 lr = lrline ( x, y, ax, ay, bx, by, 0.0 );