· delaunay2D.php ·
1<?php
/*
'=================================================================
' Descrizione.....: Routines per la triangolazione Delaunay in 2D.
5' Nome del File...: delaunay2D.php
' Data............: 25/11/2008
' Versione........: 1.0
' Linguaggio......: PHP 4 / 5
' Tradotto da.....: Mario Spada
10' Scritto da......: F. Languasco ® in linguaggio Visual Basic 6
' E-Mail..........: ZP7061 at zpyvax.vg
' DownLoads VB6...: http://www.flanguasco.org
'=================================================================
'
15' Routines di ingresso:
' DTRIS2 NPT, VCL(), IND(), NTRI, TIL(), TNBR(), IERR
' DTRIW2 NPT, VCL(), IND(), NTRI, TIL(), TNBR(), IERR
'
' Parametri: vedi descrizione all' interno delle routines.
20'
' Nota: Alla chiamata e' necessario che vengano date le dimensioni
' TIL(1 To 3, 1 To NPT * 2) e TNBR(1 To 3, 1 To NPT * 2);
' All' uscita sono validi solo i valori TIL(1 To 3, 1 To NTRI)
' e TNBR(1 To 3, 1 To NTRI).
25'
' Nota: Tutti i vettori e le matrici di queste routines
' iniziano dall' indice 1.
'
' Nota: Tutte le routines di questo modulo sono state tradotte
30' dal FORTRAN della libreria Geompack3.
'
*/
define ("Epsilon", pow(2,-52)); // Epsilon del FORTRAN per numeri in doppia precisione.
35
function zCodiciErrore($errNo)
{
//
//
40// Codici di errore originali della libreria GEOMPACK3:
//
// Error codes: Error codes IERR from 1 to 99 indicate not enough space in some array.
// Error codes IERR >= 100 mean that the input specification is incorrect,
// there is a program bug, or the tolerance TOL is too small or large (try
45// calling INITCB with a different TOLIN value).
//
//
$errorCodes = array (
1=>"not enough space in EDGE array for routine EDGHT",
50 2=>"not enough space in HOLV array for routine DSMCPR or DSPGDC",
3=>"not enough space in 2-D VCL array",
4=>"not enough space in HVL array",
5=>"not enough space in PVL, IANG arrays",
6=>"not enough space in IWK array",
55 7=>"not enough space in WK array",
8=>"not enough space in STACK array for routine SWAPEC, DTRIS2, or DTRIW2",
9=>"not enough space in TIL array",
10=>"not enough space in CWALK array for routine INTTRI",
11=>"not enough space in FC array for 3-D",
60 12=>"not enough space in BF array for 3-D",
13=>"not enough space in SVCL, SFVL arrays for routine SHRNK3",
14=>"not enough space in 3-D VCL array",
15=>"not enough space in FVL, EANG arrays for polyhedral decomposition",
16=>"not enough space in FACEP, NRML arrays",
65 17=>"not enough space in PFL array",
18=>"not enough space in HFL array",
19=>"not enough space in BTL array",
20=>"not enough space in VM array",
21=>"not enough space in HT array",
70 22=>"not enough space in FC array for K-D",
23=>"not enough space in BF array for K-D",
100=>"higher primes need to be added to routine PRIME",
200=>"abnormal return in routine DIAM2",
201=>"abnormal return in routine WIDTH2",
75 202=>"parallel lines in routine SHRNK2",
205=>"horizontal ray to right does not intersect polygon in routine ROTIPG",
206=>"out-of-range index when popping points from stack in routine VPRGHT",
207=>"out-of-range index when scanning vertices in routine VPSCNA",
208=>"out-of-range index when scanning vertices in routine VPSCNB",
80 209=>"out-of-range index when scanning vertices in routine VPSCNC",
210=>"out-of-range index when scanning vertices in routine VPSCND",
212=>"singular matrix in routine VORNBR",
215=>"unmatched separator interface edge determined by routine DSPGDC",
216=>"all edges of an outer boundary curve are specified as hole edges in input to routine DSPGDC",
85 218=>"cannot find subregion above top vertex of hole in routine JNHOLE",
219=>"angle at hole vertex in modified region is too far from PI due to use of relative tolerance in routine JNHOLE",
222=>"separator not found in routine MFDEC2 (not likely to occur if ANGSPC <= 30 degrees and ANGTOL <= 20 degrees)",
224=>"2 vertices with identical coordinates (in floating point arithmetic) detected in routine DTRIS2 or DTRIW2",
225=>"all vertices input to DTRIS2 or DTRIW2 are collinear (in fl pt arith)",
90 226=>"cycle detected in walk by routine WALKT2",
230=>"invalid (CW-oriented) triangle created in routine TMERGE",
231=>"unsuccessful search in routine FNDTRI",
300=>"unsuccessful search by routine HTSRC",
301=>"degenerate tetrahedron detected by routine BARYTH, CCSPH, or OPSIDE",
95 302=>"2 vertices with identical coordinates (in floating point arithmetic) detected in routine FRSTET, DTRIS3, DTRIW3, or ITRIS3",
303=>"all vertices input to DTRIS3 or DTRIW3 are collinear (in fl pt arith)",
304=>"all vertices input to DTRIS3 or DTRIW3 are coplanar (in fl pt arith)",
305=>"invalid configuration (due to tolerance) detected by routine SWAPES",
306=>"no visible boundary face found in routine DTRIS3 or ITRIS3",
100 307=>"cycle detected in walk by routine WALKT3",
308=>"nontransformable or nonexistent face at top of stack in routine SWAPTF",
309=>"ALPHA(4) = 0 occurs after calling BARYTH in SWAPES, SWAPMU, etc. TOL must be decreased",
310=>"unmatched edge determined by routine DSCPH",
313=>"incorrect value for variable T in routine XPGHPL",
105 314=>"unmatched edge in shrunken polyhedron in routine SHRNK3",
315=>"2 'identical' vertices with unequal coordinates in routine SHRNK3",
318=>"vertex ITOP not found in FVL array in routine INTMVG",
321=>"face oriented same way twice in routine DSPHDC or DSPHFH",
322=>"unmatched edge determined by routine DSPHDC, DSPHFH, or DSPHIH",
110 325=>"unsuccessful search for next edge of cut polygon in routine CUTFAC",
326=>"unsuccessful search for intersecting edge of face FL in CUTFAC",
327=>"at least one reflex edge has not been resolved by routine RESEDG",
328=>"value of constant MAXEV must be increased in routine CUTFAC",
330=>"value of MAXSV must be increased before calling SHRNK3 in TRIPR3",
115 331=>"invalid vertex detected from walk or linear search in routine BCDTRI",
334=>"unsuccessful search for next edge of cut polygon in routine SEPFAC",
335=>"unsuccessful search for intersecting edge of face FL in SEPFAC",
336=>"separator face not found in routine MFDEC3",
340=>"holes are not in nondecreasing face index order in routine DSPHFH",
120 341=>"hole on non-boundary face of polyhedral region detected by DSPHFH",
342=>"face containing non-coplanar vertices detected by DSPHFH",
344=>"cannot find subregion above top (or below bottom) vertex of hole in routine SPDECF",
346=>"hole polyhedron not connected to outer boundary by routine DSPHIH",
347=>"starting edge not found by routine RESHOL (because PT not in polyh)",
125 348=>"starting edge not found by routine RESHOL (because ray meets vertex)",
349=>"unsuccessful search for intersecting edge of face FMIN in RESHOL",
400=>"unsuccessful search by routine HTSRCK or HTSDLK",
401=>"degenerate simplex detected by routine BARYCK",
402=>"2 vertices with identical coordinates (in floating point arithmetic) detected in routine FRSMPX, DTRISK, DTRIWK, or DTRIMK",
130 403=>"all vertices input to DTRISK, DTRIWK, or DTRIMK are in same hyperplane (in floating point arithmetic)",
405=>"invalid configuration (due to tolerance) detected by routine SWAPHS",
406=>"no visible boundary face found in routine DTRISK",
407=>"cycle detected in walk by routine WALKTK"
);
135 return $errorCodes[$errNo];
}
//
//
function WALKT2($X, $Y, $NTRI, $Vcl, $TIL, $TNBR, &$ITRI, &$IEDG, &$IERR)
140{
//
//******************************************************************************
//
// WALKT2 walks through a 2D triangulation searching for a point.
145//
//
// Purpose:
//
// Walk through neighboring triangles of 2D (Delaunay)
150// triangulation until a triangle is found containing point (X,Y)
// or (X,Y) is found to be outside the convex hull. Search is
// guaranteed to terminate for a Delaunay triangulation, else a
// cycle may occur.
//
155// Author:
//
// Barry Joe,
// Department of Computing Science,
// University of Alberta,
160// Edmonton, Alberta, Canada T6G 2H1
// Email: oneel at pf.hnyoregn.pn
//
// Parameters:
//
165// Input, X, Y - 2D point.
//
// Input, NTRI - number of triangles in triangulation; used to detect cycle.
//
// Input, VCL(1:2,1:*) - coordinates of 2D vertices.
170//
// Input, TIL(1:3,1:*) - triangle incidence list.
//
// Input, TNBR(1:3,1:*) - triangle neighbor list.
//
175// Input/output, ITRI. On input, the index of triangle to begin search at.
// On output, the index of triangle that search ends at.
//
// Output, IEDG - 0 if (X,Y) is in the interior of triangle ITRI; I = 1,
// 2, or 3 if (X,Y) is on interior of edge I of ITRI;
180// I = 4, 5, or 6 if (X,Y) is (nearly) vertex I-3 of ITRI;
// I = -1, -2, or -3 if (X,Y) is outside convex hull due
// to walking past edge -I of triangle ITRI.
//
// Output, integer IERR, error flag, which is zero unless an error occurred.
185//
// Dim I&, a&, b&, c&, cnt&
// Dim Alpha#, Beta#, Gamma#, det#, Dx#, dXa#, dXb#, Dy#, dYa#, dYb#, tol#
//
// Use barycentric coordinates to determine where (X,Y) is located.
190//
$IERR = 0;
$tol = 100 * Epsilon; //(Tol)
$cnt = 0;
$IEDG = 0;
195//
//10 CONTINUE
//
while($cnt <= $NTRI){
$cnt++;
200//
$a = $TIL[1][$ITRI];
$b = $TIL[2][$ITRI];
$c = $TIL[3][$ITRI];
$dXa = $Vcl[1][$a] - $Vcl[1][$c];
205 $dYa = $Vcl[2][$a] - $Vcl[2][$c];
$dXb = $Vcl[1][$b] - $Vcl[1][$c];
$dYb = $Vcl[2][$b] - $Vcl[2][$c];
$Dx = $X - $Vcl[1][$c];
$Dy = $Y - $Vcl[2][$c];
210 $det = $dXa * $dYb - $dYa * $dXb;
$Alpha = ($Dx * $dYb - $Dy * $dXb) / $det;
$Beta = ($dXa * $Dy - $dYa * $Dx) / $det;
$Gamma = 1 - $Alpha - $Beta;
//
215 if ($Alpha > $tol && $Beta > $tol && $Gamma > $tol) {
return true;
}
//
elseif ($Alpha < ($tol * -1)) {
220 $I = $TNBR[2][$ITRI];
if ($I <= 0) {
$IEDG = -2;
return true;
}
225 }
//
elseif ($Beta < ($tol * -1)) {
$I = $TNBR[3][$ITRI];
if ($I <= 0) {
230 $IEDG = -3;
return true;
}
}
//
235 elseif ($Gamma < ($tol * -1)) {
$I = $TNBR[1][$ITRI];
if ($I <= 0) {
$IEDG = -1;
return true;
240 }
}
//
elseif ($Alpha <= $tol) {
if ($Beta <= $tol){
245 $IEDG = 6;
}
elseif ($Gamma <= $tol){
$IEDG = 5;
}
250 else{
$IEDG = 2;
}
return true;
}
255//
elseif ($Beta <= $tol) {
if ($Gamma <= $tol) {
$IEDG = 4;
}
260 else{
$IEDG = 3;
}
return true;
}
265//
else{
$IEDG = 1;
return true;
}
270
$ITRI = $I;
} //GoTo 10
//
$IERR = 226;
275 return true;
//
}
280
function DTRIW2 ($NPT, $Vcl, $Ind, &$NTRI, &$TIL, &$TNBR, &$IERR)
{
//
//******************************************************************************
285//
// DTRIW2 constructs a Delaunay triangulation of vertices in 2D.
//
//
// Purpose:
290//
// Construct Delaunay triangulation of 2D vertices using
// incremental approach and diagonal edge swaps. Vertices are
// inserted one at a time in order given by IND array. The initial
// triangles created due to a new vertex are obtained by a walk
295// through the triangulation until location of vertex is known.
//
// Author:
//
// Barry Joe,
300// Department of Computing Science,
// University of Alberta,
// Edmonton, Alberta, Canada T6G 2H1
// Email: oneel at pf.hnyoregn.pn
//
305// Parameters:
//
// Input, NPT - number of 2D points (vertices).
//
// Input, MAXST - maximum size available for STACK array; should be about
310// NPT to be safe, but MAX(10,2*LOG2(NPT)) usually enough.
//
// Input, VCL(1:2,1:*) - coordinates of 2D vertices.
//
// Input, IND(1:NPT) - indices in VCL of vertices to be triangulated;
315// vertices are inserted in order given by this array
//
// Output, NTRI - number of triangles in triangulation; equal to
// 2*NPT - NB - 2 where NB = number of boundary vertices.
//
320// Output, TIL(1:3,1:NTRI) - triangle incidence list; elements are indices
// of VCL; vertices of triangles are in counter clockwise order.
//
// Output, TNBR(1:3,1:NTRI) - triangle neighbor list; positive elements
// are indices of TIL; negative elements are used for links
325// of counter clockwise linked list of boundary edges; LINK = -(3*I + J-1)
// where I, J = triangle, edge index; TNBR(J,I) refers to
// the neighbor along edge from vertex J to J+1 (mod 3).
//
// Workspace, STACK(1:MAXST) - used for stack of triangles for which
330// circumcircle test must be made
//
// Output, integer IERR, error flag, which is zero unless an error occurred.
//
// Dim cmax#, tol#
335// Dim I&, I3&, J&, e&, EM1&, EP1&, M&, M1&, M2&, M3&, MAXST&, N&
// Dim L&, LEDG&, LR&, LTRI&, REDG&, RTRI&, T&, Top&
//
$MAXST = $NPT;
$STACK = array();
340// ReDim STACK&(1 To MAXST)
//
// Const msglvl& = 0
define ("msglvl",0);
345
//
// Determine initial triangle.
//
$IERR = 0;
350 $tol = 100 * Epsilon; //(Tol)
//
$M1 = $Ind[1];
$M2 = $Ind[2];
//
355 for ($J=1;$J<=2;$J++){
$cmax = DMAX1(abs($Vcl[$J][$M1]),abs($Vcl[$J][$M2]));
If (abs($Vcl[$J][$M1] - $Vcl[$J][$M2]) > $tol * $cmax && $cmax > $tol) {
break;
360 }
else{
$IERR = 224;
return true;
}
365 }
//
//20 //CONTINUE
//
$I3 = 3;
370 $LR = 0;
//30 //CONTINUE
//
while ($LR == 0){
375 if ($I3 > $NPT){
$IERR = 225;
return true;
}
//
380 $M = $Ind[$I3];
$LR = LRLINE($Vcl[1][$M], $Vcl[2][$M], $Vcl[1][$M1], $Vcl[2][$M1], $Vcl[1][$M2], $Vcl[2][$M2], 0);
//
if ($LR == 0) {
$I3++; // I3 = I3 + 1
385// GoTo 30
}
}
390//
if ($I3 <> 3) {
$Ind[$I3] = $Ind[3];
$Ind[3] = $M;
}
395//
$NTRI = 1;
//
if ($LR == -1) {
$TIL[1][1] = $M1;
400 $TIL[2][1] = $M2;
}
else {
$TIL[1][1] = $M2;
$TIL[2][1] = $M1;
405 }
//
$TIL[3][1] = $M;
$TNBR[1][1] = -4;
$TNBR[2][1] = -5;
410 $TNBR[3][1] = -3;
//
// If (msglvl = 4) Then
415// write ( *,600) 1,vcl(1,m1),vcl(2,m1),vcl(1,m2),vcl(2,m2)
// write ( *,600) 1,vcl(1,m2),vcl(2,m2),vcl(1,m),vcl(2,m)
// write ( *,600) 1,vcl(1,m),vcl(2,m),vcl(1,m1),vcl(2,m1)
// End If
//
420// Insert vertices one at a time from anywhere, walk through
// triangulation to determine location of new vertex, and apply
// diagonal edge swaps until Delaunay triangulation of vertices
// (so far) is obtained.
//
425 $Top = 0;
//
for ($I=4;$I<=$NPT;$I++) {
430// If (msglvl = 4) Then
// write ( *,600) i
// End If
//
$M = $Ind[$I];
435 $RTRI = $NTRI;
$res = WALKT2($Vcl[1][$M], $Vcl[2][$M], $NTRI, $Vcl, $TIL, $TNBR, $RTRI, $REDG, $IERR);
//
if ($REDG == 0) {
$M1 = $TIL[1][$RTRI];
440 $M2 = $TIL[2][$RTRI];
$M3 = $TIL[3][$RTRI];
$TIL[3][$RTRI] = $M;
//
if ($TNBR[1][$RTRI] > 0) {
445 $Top = 1;
$STACK[$Top] = $RTRI;
}
//
$NTRI++;
450 $TIL[1][$NTRI] = $M2;
$TIL[2][$NTRI] = $M3;
$TIL[3][$NTRI] = $M;
$N = $TNBR[2][$RTRI];
$TNBR[1][$NTRI] = $N;
455//
if ($N > 0) {
If ($TNBR[1][$N] == $RTRI) {
$TNBR[1][$N] = $NTRI;
}
460 elseif ($TNBR[2][$N] == $RTRI) {
$TNBR[2][$N] = $NTRI;
}
else {
$TNBR[3][$N] = $NTRI;
465 }
//
$Top++;
$STACK[$Top] = $NTRI;
}
470//
$NTRI++;
$TIL[1][$NTRI] = $M3;
$TIL[2][$NTRI] = $M1;
$TIL[3][$NTRI] = $M;
475 $N = $TNBR[3][$RTRI];
$TNBR[1][$NTRI] = $N;
//
if ($N > 0) {
if ($TNBR[1][$N] == $RTRI) {
480 $TNBR[1][$N] = $NTRI;
}
elseif ($TNBR[2][$N] == $RTRI) {
$TNBR[2][$N] = $NTRI;
}
485 else {
$TNBR[3][$N] = $NTRI;
}
//
$Top++;
490 $STACK[$Top] = $NTRI;
}
//
$TNBR[2][$RTRI] = $NTRI - 1;
$TNBR[3][$RTRI] = $NTRI;
495 $TNBR[2][$NTRI - 1] = $NTRI;
$TNBR[3][$NTRI - 1] = $RTRI;
$TNBR[2][$NTRI] = $RTRI;
$TNBR[3][$NTRI] = $NTRI - 1;
//
500 if ($TNBR[1][$NTRI - 1] <= 0) {
$T = $RTRI;
$e = 1;
//
//40 //CONTINUE
505
while($TNBR[$e][$T] > 0){
//
// if ($TNBR[$e][$T] > 0) {
$T = $TNBR[$e][$T];
510//
if ($TIL[1][$T] == $M2) {
$e = 3;
}
elseif ($TIL[2][$T] == $M2) {
515 $e = 1;
}
else {
$e = 2;
}
520//
//GoTo 40
// }
}
//
525 $TNBR[$e][$T] = -3 * $NTRI + 3;
}
//
if ($TNBR[1][$NTRI] <= 0) {
$T = $NTRI - 1;
530 $e = 1;
//
//50 //CONTINUE
//
while($TNBR[$e][$T] > 0){
535// If (TNBR(e, T) > 0) Then
$T = $TNBR[$e][$T];
//
if ($TIL[1][$T] == $M3) {
$e = 3;
540 }
elseif ($TIL[2][$T] == $M3) {
$e = 1;
}
else {
545 $e = 2;
}
//
// GoTo 50
// End If
550 }
//
$TNBR[$e][$T] = -3 * $NTRI;
//
}
555//
// if (msglvl = 4) then
// write ( *,600) 1,vcl(1,m),vcl(2,m),vcl(1,m1),vcl(2,m1)
// write ( *,600) 1,vcl(1,m),vcl(2,m),vcl(1,m2),vcl(2,m2)
// write ( *,600) 1,vcl(1,m),vcl(2,m),vcl(1,m3),vcl(2,m3)
560// End If
//
}
elseif ($REDG < 0) {
$REDG = -$REDG;
565 $LTRI = 0;
//
$res = VBEDG($Vcl[1][$M], $Vcl[2][$M], $Vcl, $TIL, $TNBR, $LTRI, $LEDG, $RTRI, $REDG);
//
$N = $NTRI + 1;
570 $L = -$TNBR[$LEDG][$LTRI];
//
//60 //CONTINUE
$goto = true;
while ($goto){
575//
$T = intval($L / 3);
$e = ($L % 3) + 1;
$L = -$TNBR[$e][$T];
$M2 = $TIL[$e][$T];
580//
if ($e <= 2) {
$M1 = $TIL[$e + 1][$T];
}
else {
585 $M1 = $TIL[1][$T];
}
//
$NTRI++;
$TNBR[$e][$T] = $NTRI;
590 $TIL[1][$NTRI] = $M1;
$TIL[2][$NTRI] = $M2;
$TIL[3][$NTRI] = $M;
$TNBR[1][$NTRI] = $T;
$TNBR[2][$NTRI] = $NTRI - 1;
595 $TNBR[3][$NTRI] = $NTRI + 1;
$Top = $Top + 1;
//
if ($Top > $MAXST) {
$IERR = 8;
600 break; //GoTo 100
}
//
$STACK[$Top] = $NTRI;
//
605// If (msglvl = 4) Then
// write ( *,600) 1,vcl(1,m),vcl(2,m),vcl(1,m2),vcl(2,m2)
// End If
//
// If (T <> RTRI Or e <> REDG) Then GoTo 60
610 if ($T <> $RTRI || $e <> $REDG){
// Then GoTo 60
$goto = true;
}
else {
615 $goto = false;
}
}
if ($IERR == 8) {
620 break;
}
//
// If (msglvl = 4) Then
625// write ( *,600) 1,vcl(1,m),vcl(2,m), vcl(1,m1),vcl(2,m1)
// End If
//
$TNBR[$LEDG][$LTRI] = -3 * $N - 1;
$TNBR[2][$N] = -3 * $NTRI - 2;
630 $TNBR[3][$NTRI] = -$L;
//
}
elseif ($REDG <= 3) {
$M1 = $TIL[$REDG][$RTRI];
635//
if ($REDG == 1) {
$e = 2;
$EP1 = 3;
}
640 elseif ($REDG == 2) {
$e = 3;
$EP1 = 1;
}
else {
645 $e = 1;
$EP1 = 2;
}
//
$M2 = $TIL[$e][$RTRI];
650 $TIL[$e][$RTRI] = $M;
$M3 = $TIL[$EP1][$RTRI];
//
if ($TNBR[$EP1][$RTRI] > 0) {
$Top = 1;
655 $STACK[$Top] = $RTRI;
}
//
$NTRI++;
$TIL[1][$NTRI] = $M;
660 $TIL[2][$NTRI] = $M2;
$TIL[3][$NTRI] = $M3;
$N = $TNBR[$e][$RTRI];
$TNBR[2][$NTRI] = $N;
$TNBR[3][$NTRI] = $RTRI;
665 $TNBR[$e][$RTRI] = $NTRI;
//
if ($N > 0) {
if ($TNBR[1][$N] == $RTRI) {
$TNBR[1][$N] = $NTRI;
670 }
elseif ($TNBR[2][$N] == $RTRI) {
$TNBR[2][$N] = $NTRI;
}
else {
675 $TNBR[3][$N] = $NTRI;
}
//
$Top++;
$STACK[$Top] = $NTRI;
680 }
//
// If (msglvl = 4) Then
// write ( *,600) 1,vcl(1,m),vcl(2,m),vcl(1,m3),vcl(2,m3)
// End If
685//
$LTRI = $TNBR[$REDG][$RTRI];
//
if ($LTRI <= 0) {
$TNBR[1][$NTRI] = $LTRI;
690 $TNBR[$REDG][$RTRI] = -3 * $NTRI;
if ($TNBR[2][$NTRI] <= 0) { $TNBR[1][$NTRI] = -3 * $NTRI - 1;}
}
//
else {
695 $TNBR[1][$NTRI] = $NTRI + 1;
$TNBR[$REDG][$RTRI] = $LTRI;
//
if ($TIL[1][$LTRI] == $M2) {
$LEDG = 1;
700 $EM1 = 2;
$e = 3;
}
elseif ($TIL[2][$LTRI] == $M2) {
$LEDG = 2;
705 $EM1 = 3;
$e = 1;
}
else {
$LEDG = 3;
710 $EM1 = 1;
$e = 2;
}
//
$TIL[$LEDG][$LTRI] = $M;
715 $M3 = $TIL[$e][$LTRI];
//
if ($TNBR[$EM1][$LTRI] > 0) {
$Top++;
$STACK[$Top] = $LTRI;
720 }
//
$NTRI++;
$TIL[1][$NTRI] = $M2;
$TIL[2][$NTRI] = $M;
725 $TIL[3][$NTRI] = $M3;
$TNBR[1][$NTRI] = $NTRI - 1;
$TNBR[2][$NTRI] = $LTRI;
$N = $TNBR[$e][$LTRI];
$TNBR[3][$NTRI] = $N;
730 $TNBR[$e][$LTRI] = $NTRI;
//
if ($N > 0) {
if ($TNBR[1][$N] == $LTRI) {
$TNBR[1][$N] = $NTRI;
735 }
elseif ($TNBR[2][$N] == $LTRI) {
$TNBR[2][$N] = $NTRI;
}
else {
740 $TNBR[3][$N] = $NTRI;
}
//
$Top++;
$STACK[$Top] = $NTRI;
745 }
//
// If (msglvl = 4) Then
// write ( *,600) 1,vcl(1,m),vcl(2,m), vcl(1,m3),vcl(2,m3)
// End If
750//
if ($TNBR[2][$NTRI - 1] <= 0) {
$T = $NTRI;
$e = 3;
//
755//70 //CONTINUE
//
while ($TNBR[$e][$T] > 0) {
// If (TNBR(e, T) > 0) Then
$T = $TNBR[$e][$T];
760//
if ($TIL[1][$T] == $M2) {
$e = 3;
}
elseif ($TIL[2][$T] == $M2) {
765 $e = 1;
}
else {
$e = 2;
}
770//
// GoTo 70
// End If
}
//
775 $TNBR[$e][$T] = -3 * $NTRI + 2;
}
//
if ($TNBR[3][$NTRI] <= 0) {
$T = $LTRI;
780//
if ($LEDG <= 2) {
$e = $LEDG + 1;
}
else {
785 $e = 1;
}
//
//80 //CONTINUE
//
790 while ($TNBR[$e][$T] > 0) {
// If (TNBR(e, T) > 0) Then
$T = $TNBR[$e][$T];
//
if ($TIL[1][$T] == $M3) {
795 $e = 3;
}
elseif ($TIL[2][$T] == $M3) {
$e = 1;
}
800 else {
$e = 2;
}
//
// GoTo 80
805//
// End If
}
//
$TNBR[$e][$T] = -3 * $NTRI - 2;
810 }
}
}
//
else {
815 $IERR = 224;
break;
// GoTo 100
}
//
820 $BTRI = 0;
$BEDG = 0;
$res = SWAPEC($M, $Top, $MAXST, $BTRI, $BEDG, $Vcl, $TIL, $TNBR, $STACK, $IERR);
//
if ($IERR <> 0) {
825 break;
// Exit For
}
}
//
830
//100 //CONTINUE
//
if ($I3 <> 3) {
$T = $Ind[$I3];
835 $Ind[$I3] = $Ind[3];
$Ind[3] = $T;
}
//
// if (msglvl = 4) then
840// write ( *,600) npt+1
// End If
//600 format (1x,i7,4f15.7)
//
845//
//
//end function
}
function DIAEDG($X0, $Y0, $X1, $Y1, $X2, $Y2, $X3, $Y3)
850{
//
//******************************************************************************
//
// DIAEDG determines which diagonal to use in a quadrilateral.
855//
//
// Purpose:
//
// Determine whether 02 or 13 is the diagonal edge chosen
860// based on the circumcircle criterion, where (X0,Y0), (X1,Y1),
// (X2,Y2), (X3,Y3) are the vertices of a simple quadrilateral
// in counterclockwise order.
//
// Author:
865//
// Barry Joe,
// Department of Computing Science,
// University of Alberta,
// Edmonton, Alberta, Canada T6G 2H1
870// Email: oneel at pf.hnyoregn.pn
//
// Parameters:
//
// Input, X0, Y0, X1, Y1, X2, Y2, X3, Y3 - vertex coordinates.
875//
// Output, DIAEDG -
// 1, if diagonal edge 02 is chosen, i.e. 02 is inside.
// quadrilateral + vertex 3 is outside circumcircle 012
// -1, if diagonal edge 13 is chosen, i.e. 13 is inside.
880// quadrilateral + vertex 0 is outside circumcircle 123
// 0, if four vertices are cocircular.
//
// Dim ca#, cb#, dx10#, dx12#, dx30#, dx32#, dy10#, dy12#, dy30#, dy32#
// Dim s#, tol#, tola#, tolb#
885//
$tol = 100 * Epsilon; //(tol)
$dx10 = $X1 - $X0;
$dy10 = $Y1 - $Y0;
$dx12 = $X1 - $X2;
890 $dy12 = $Y1 - $Y2;
$dx30 = $X3 - $X0;
$dy30 = $Y3 - $Y0;
$dx32 = $X3 - $X2;
$dy32 = $Y3 - $Y2;
895 $tola = $tol * DMAX1(abs($dx10), abs($dy10), abs($dx30), abs($dy30));
$tolb = $tol * DMAX1(abs($dx12), abs($dy12), abs($dx32), abs($dy32));
$ca = $dx10 * $dx30 + $dy10 * $dy30;
$cb = $dx12 * $dx32 + $dy12 * $dy32;
//
900 if ($ca > $tola && $cb > $tolb) {
return -1;
}
elseif ($ca < -$tola && $cb < -$tolb) {
return 1;
905 }
else {
$tola = DMAX1($tola, $tolb);
$s = ($dx10 * $dy30 - $dx30 * $dy10) * $cb + ($dx32 * $dy12 - $dx12 * $dy32) * $ca;
910 if ($s > $tola) {
return -1;
}
elseif ($s < -$tola) {
return 1;
915 }
else {
return 0;
}
}
920//
//
//
}
function SWAPEC($I, &$Top, $MAXST, &$BTRI, &$BEDG, $Vcl, &$TIL, &$TNBR, $STACK, &$IERR)
925{
//
//******************************************************************************
//
// SWAPEC swaps diagonal edges in a 2D triangulation
930//
//
// Purpose:
//
// Swap diagonal edges in 2D triangulation based on empty
935// circumcircle criterion until all triangles are Delaunay, given
// that I is index of new vertex added to triangulation.
//
// Author:
//
940// Barry Joe,
// Department of Computing Science,
// University of Alberta,
// Edmonton, Alberta, Canada T6G 2H1
// Email: oneel at pf.hnyoregn.pn
945//
// Parameters:
//
// Input, I - index in VCL of new vertex.
//
950// Input/output, TOP - index of top of stack, >= 0.
//
// Input, MAXST - maximum size available for STACK array.
//
// Input/output, BTRI, BEDG - if positive, these are triangle and edge
955// index of a boundary edge whose updated indices must be recorded.
//
// Input, VCL(1:2,1:*) - coordinates of 2D vertices.
//
// Input/output, TIL(1:3,1:*) - triangle incidence list.
960//
// Input/output, TNBR(1:3,1:*) - triangle neighbor list; negative values are
// used for links of counter clockwise linked list of boundary edges;
// LINK = -(3*I + J-1) where I, J = triangle, edge index.
//
965// Input, STACK(1:TOP) - index of initial triangles (involving vertex I)
// put in stack; the edges opposite I should be in interior.
//
// Workspace, STACK(TOP+1:MAXST) - used as stack.
//
970// Output, integer IERR, error flag, which is zero unless an error occurred.
//
// Dim L&, a&, b&, c&, e&, ee&, EM1&, EP1&
// Dim F&, fm1&, fp1&, R&, s&, swap&, T&, tt&, u&
// Dim X#, Y#
975//
// Const msglvl& = 0
define ("msglvl",0);
//
// Determine whether triangles in stack are Delaunay, and swap
980// diagonal edge of convex quadrilateral if not.
//
$IERR = 0;
$X = $Vcl[1][$I];
$Y = $Vcl[2][$I];
985//
//10 //CONTINUE
while ($Top > 0) {
//
990// If (Top <= 0) Then Exit Sub
$T = $STACK[$Top];
$Top--;
//
if ($TIL[1][$T] == $I) {
995 $e = 2;
$b = $TIL[3][$T];
}
elseif ($TIL[2][$T] == $I) {
$e = 3;
1000 $b = $TIL[1][$T];
}
else {
$e = 1;
$b = $TIL[2][$T];
1005 }
//
$a = $TIL[$e][$T];
$u = $TNBR[$e][$T];
//
1010 if ($TNBR[1][$u] == $T) {
$F = 1;
$c = $TIL[3][$u];
}
elseif ($TNBR[2][$u] == $T) {
1015 $F = 2;
$c = $TIL[1][$u];
}
else {
$F = 3;
1020 $c = $TIL[2][$u];
}
//
$swap = DIAEDG($X, $Y, $Vcl[1][$a], $Vcl[2][$a], $Vcl[1][$c], $Vcl[2][$c], $Vcl[1][$b], $Vcl[2][$b]);
//
1025 if ($swap == 1) {
$EM1 = $e - 1;
if ($EM1 == 0) { $EM1 = 3; }
$EP1 = $e + 1;
if ($EP1 == 4) { $EP1 = 1; }
1030 $fm1 = $F - 1;
if ($fm1 == 0) { $fm1 = 3; }
$fp1 = $F + 1;
if ($fp1 == 4) { $fp1 = 1; }
$TIL[$EP1][$T] = $c;
1035 $TIL[$fp1][$u] = $I;
$R = $TNBR[$EP1][$T];
$s = $TNBR[$fp1][$u];
$TNBR[$EP1][$T] = $u;
$TNBR[$fp1][$u] = $T;
1040 $TNBR[$e][$T] = $s;
$TNBR[$F][$u] = $R;
//
if ($TNBR[$fm1][$u] > 0) {
$Top++;
1045 $STACK[$Top] = $u;
}
//
if ($s > 0) {
if ($TNBR[1][$s] == $u) {
1050 $TNBR[1][$s] = $T;
}
elseif ($TNBR[2][$s] == $u) {
$TNBR[2][$s] = $T;
}
1055 else {
$TNBR[3][$s] = $T;
}
//
$Top++;
1060//
if ($Top > $MAXST) {
$IERR = 8;
return true;
}
1065//
$STACK[$Top] = $T;
//
}
else {
1070 if ($u == $BTRI && $fp1 == $BEDG) {
$BTRI = $T;
$BEDG = $e;
}
//
1075 $L = -(3 * $T + $e - 1);
$tt = $T;
$ee = $EM1;
//
//20 //CONTINUE
1080 while ($TNBR[$ee][$tt] > 0) {
//
// If (TNBR(ee, tt) > 0) Then
$tt = $TNBR[$ee][$tt];
//
1085 if ($TIL[1][$tt] == $a) {
$ee = 3;
}
elseif ($TIL[2][$tt] == $a) {
$ee = 1;
1090 }
else {
$ee = 2;
}
//
1095// GoTo 20
// End If
}
//
$TNBR[$ee][$tt] = $L;
1100 }
//
if ($R > 0) {
if ($TNBR[1][$R] == $T) {
$TNBR[1][$R] = $u;
1105 }
elseif ($TNBR[2][$R] == $T) {
$TNBR[2][$R] = $u;
}
else {
1110 $TNBR[3][$R] = $u;
}
//
}
else {
1115 if ($T == $BTRI && $EP1 == $BEDG) {
$BTRI = $u;
$BEDG = $F;
}
//
1120 $L = -(3 * $u + $F - 1);
$tt = $u;
$ee = $fm1;
//
//30 //CONTINUE
1125//
while ($TNBR[$ee][$tt] > 0) {
// If (TNBR(ee, tt) > 0) Then
$tt = $TNBR[$ee][$tt];
//
1130 if ($TIL[1][$tt] == $b) {
$ee = 3;
}
elseif ($TIL[2][$tt] == $b) {
$ee = 1;
1135 }
else {
$ee = 2;
}
//
1140// GoTo 30
// End If
}
//
$TNBR[$ee][$tt] = $L;
1145 }
//
// If (msglvl = 4) Then
// write ( *,600) 2,vcl(1,a),vcl(2,a), &
// vcl(1,b),vcl(2,b),x,y,vcl(1,c),vcl(2,c)
1150// End If
//600 format (1x,i7,4f15.7/8x,4f15.7)
//
}
//
1155// GoTo 10
}
return true;
//
}
1160function VBEDG($X, $Y, $Vcl, $TIL, $TNBR, &$LTRI, &$LEDG, &$RTRI, &$REDG)
{
//
//******************************************************************************
//
1165// VBEDG determines the boundary edges of a 2D triangulation.
//
//
// Purpose:
//
1170// Determine boundary edges of 2D triangulation which are
// visible from point (X,Y) outside convex hull.
//
// Author:
//
1175// Barry Joe,
// Department of Computing Science,
// University of Alberta,
// Edmonton, Alberta, Canada T6G 2H1
// Email: oneel at pf.hnyoregn.pn
1180//
// Parameters:
//
// Input, X, Y - 2D point outside convex hull.
//
1185// Input, VCL(1:2,1:*) - coordinates of 2D vertices.
//
// Input, TIL(1:3,1:*) - triangle incidence list.
//
// Input, TNBR(1:3,1:*) - triangle neighbor list; negative values are
1190// used for links of counter clockwise linked list of boundary edges;
// LINK = -(3*I + J-1) where I, J = triangle, edge index.
//
// Input, LTRI, LEDG - if LTRI /= 0 then they are assumed to be as defined
// below and are not changed, else they are updated. LTRI is
1195// the index of boundary triangle to left of leftmost boundary
// triangle visible from (X,Y). LEDG is the boundary edge of triangle
// LTRI to left of leftmost boundary edge visible from (X,Y)
//
// Input/output, RTRI. On input, the index of boundary triangle to begin
1200// search at. On output, the index of rightmost boundary triangle
// visible from (X,Y)
//
// Input/output, REDG. On input, the edge of triangle RTRI that is
// visible from (X,Y). On output, the edge of triangle RTRI that is
1205// visible from (X,Y)
//
// Dim a&, b&, e&, L&, LR&, T&
// Dim LDONE As Boolean
//
1210// Find rightmost visible boundary edge using links, then possibly
// leftmost visible boundary edge using triangle neighbor info.
//
if ($LTRI == 0) {
$LDONE = false;
1215 $LTRI = $RTRI;
$LEDG = $REDG;
}
else {
$LDONE = true;
1220 }
//
//10 //CONTINUE
//
$GoTo = true;
1225 while ($GoTo) {
$L = -$TNBR[$REDG][$RTRI];
$T = intval($L / 3);
$e = ($L % 3) + 1;
$a = $TIL[$e][$T];
1230//
if ($e <= 2) {
$b = $TIL[$e + 1][$T];
}
else {
1235 $b = $TIL[1][$T];
}
//
$LR = LRLINE($X, $Y, $Vcl[1][$a], $Vcl[2][$a], $Vcl[1][$b], $Vcl[2][$b], 0);
//
1240 if ($LR > 0) {
$RTRI = $T;
$REDG = $e;
$GoTo = true;
// GoTo 10
1245 }
else {
$GoTo = false;
}
}
1250//
if ($LDONE) { return true; }
//
$T = $LTRI;
$e = $LEDG;
1255//
//20 //CONTINUE
$GoTo = true;
while ($GoTo) {
//
1260 $b = $TIL[$e][$T];
//
if ($e >= 2) {
$e--;
}
1265 else {
$e = 3;
}
//
while ($TNBR[$e][$T] > 0) {
1270 $T = $TNBR[$e][$T];
//
if ($TIL[1][$T] == $b) {
$e = 3;
}
1275 elseif ($TIL[2][$T] == $b) {
$e = 1;
}
else {
$e = 2;
1280 }
}
//
$a = $TIL[$e][$T];
$LR = LRLINE($X, $Y, $Vcl[1][$a], $Vcl[2][$a], $Vcl[1][$b], $Vcl[2][$b], 0);
1285
if ($LR > 0) {
$GoTo = true;
// GoTo 20
}
1290 else {
$GoTo = false;
}
}
//
1295 $LTRI = $T;
$LEDG = $e;
//
//
//
1300}
function LRLINE($XU, $YU, $XV1, $YV1, $XV2, $YV2, $dV)
{
//
//******************************************************************************
1305//
// LRLINE determines whether a point is left, right, or on a directed line.
//
//
// Purpose:
1310//
// Determine whether a point is to the left of, right of,
// or on a directed line parallel to a line through given points.
//
// Author:
1315//
// Barry Joe,
// Department of Computing Science,
// University of Alberta,
// Edmonton, Alberta, Canada T6G 2H1
1320// Email: oneel at pf.hnyoregn.pn
//
// Parameters:
//
// Input, XU, YU, XV1, YV1, XV2, YV2 - vertex coordinates; the directed
1325// 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.
//
1330// Input, DV - signed distance (positive for left).
//
// Output, LRLINE - +1, 0, or -1 depending on whether (XU,YU) is
// to the right of, on, or left of the directed line
// (0 if line degenerates to a point).
1335//
// Dim Dx#, dxu#, Dy#, dyu#, T#, tol#, tolabs#
//
$tol = 100 * Epsilon; //(tol)
$Dx = $XV2 - $XV1;
1340 $Dy = $YV2 - $YV1;
$dxu = $XU - $XV1;
$dyu = $YU - $YV1;
$tolabs = $tol * DMAX1(abs($Dx), abs($Dy), abs($dxu), abs($dyu), abs($dV));
$T = $Dy * $dxu - $Dx * $dyu;
1345//
if ($dV <> 0) {
$T = $T + $dV * sqrt($Dx * $Dx + $Dy * $Dy);
}
//
1350// LRLINE = Int(SIGN(1#, T))
$result = intval(SIGN(1, $T));
//
if (abs($T) <= $tolabs) {
// LRLINE = 0
1355 $result = 0;
}
//
return $result;
//
1360//
}
function DMIN1()
{
$arg_list = func_get_args();
1365// Implementa la funzione DMIN1(D1, D2, ...DN) del FORTRAN:
// Dim J&, J1&, J2&, vDMin As Variant
/*
J1 = LBound(Vd)
1370 J2 = UBound(Vd)
vDMin = Vd(J1)
For J = J1 + 1 To J2
If vDMin > Vd(J) Then vDMin = Vd(J)
Next J
1375'
DMIN1 = CDbl(vDMin)
'
*/
return min($arg_list);
1380
}
function DMAX1()
{
1385//
// Implementa la funzione DMAX1(D1, D2, ...DN) del FORTRAN:
//
// Dim J&, J1&, J2&, vDMax As Variant
//
1390
$arg_list = func_get_args();
// J1 = LBound(Vd)
// J2 = UBound(Vd)
1395// vDMax = Vd(J1)
// For J = J1 + 1 To J2
// If vDMax < Vd(J) Then vDMax = Vd(J)
// Next J
//
1400// DMAX1 = CDbl(vDMax)
//
return max($arg_list);
}
function SIGN($dV, $DS)
1405{
//
// Ritorna il valore assoluto di dV con il segno di dS.
// Implementa la funzione SIGN del FORTRAN:
//
1410 if ($DS < 0) {
return -1 * abs($dV);
}
else {
return abs($dV);
1415 }
}
?>
· delaunay2D.php ends ·