1 /*=========================================================================
2
3 Program: Visualization Toolkit
4 Module: vtkTriangle.cxx
5
6 Copyright (c) Ken Martin, Will Schroeder, Bill Lorensen
7 All rights reserved.
8 See Copyright.txt or http://www.kitware.com/Copyright.htm for details.
9
10 This software is distributed WITHOUT ANY WARRANTY; without even
11 the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
12 PURPOSE. See the above copyright notice for more information.
13
14 =========================================================================*/
15 #include "vtkTriangle.h"
16
17 #include "vtkCellArray.h"
18 #include "vtkCellData.h"
19 #include "vtkLine.h"
20 #include "vtkMath.h"
21 #include "vtkObjectFactory.h"
22 #include "vtkPlane.h"
23 #include "vtkPointData.h"
24 #include "vtkIncrementalPointLocator.h"
25 #include "vtkPoints.h"
26 #include "vtkPolygon.h"
27 #include "vtkQuadric.h"
28
29 vtkStandardNewMacro(vtkTriangle);
30
31 //----------------------------------------------------------------------------
32 // Construct the triangle with three points.
vtkTriangle()33 vtkTriangle::vtkTriangle()
34 {
35 this->Points->SetNumberOfPoints(3);
36 this->PointIds->SetNumberOfIds(3);
37 for (int i = 0; i < 3; i++)
38 {
39 this->Points->SetPoint(i, 0.0, 0.0, 0.0);
40 this->PointIds->SetId(i,0);
41 }
42 this->Line = vtkLine::New();
43 }
44
45 //----------------------------------------------------------------------------
~vtkTriangle()46 vtkTriangle::~vtkTriangle()
47 {
48 this->Line->Delete();
49 }
50
51 //----------------------------------------------------------------------------
52 // This function simply calls the static function:
53 // vtkTriangle::TriangleArea(double p1[3], double p2[3], double p3[3])
54 // with the appropriate parameters from the instantiated vtkTriangle.
ComputeArea()55 double vtkTriangle::ComputeArea()
56 {
57 double p0[3];
58 double p1[3];
59 double p2[3];
60 this->GetPoints()->GetPoint(0, p0);
61 this->GetPoints()->GetPoint(1, p1);
62 this->GetPoints()->GetPoint(2, p2);
63 return vtkTriangle::TriangleArea(p0, p1, p2);
64 }
65
66
67 //----------------------------------------------------------------------------
68 // Create a new cell and copy this triangle's information into the cell.
69 // Returns a poiner to the new cell created.
EvaluatePosition(double x[3],double * closestPoint,int & subId,double pcoords[3],double & dist2,double * weights)70 int vtkTriangle::EvaluatePosition(double x[3], double* closestPoint,
71 int& subId, double pcoords[3],
72 double& dist2, double *weights)
73 {
74 int i, j;
75 double pt1[3], pt2[3], pt3[3], n[3], fabsn;
76 double rhs[2], c1[2], c2[2];
77 double det;
78 double maxComponent;
79 int idx=0, indices[2];
80 double dist2Point, dist2Line1, dist2Line2;
81 double *closest, closestPoint1[3], closestPoint2[3], cp[3];
82
83 subId = 0;
84
85 // Get normal for triangle, only the normal direction is needed, i.e. the
86 // normal need not be normalized (unit length)
87 //
88 this->Points->GetPoint(1, pt1);
89 this->Points->GetPoint(2, pt2);
90 this->Points->GetPoint(0, pt3);
91
92 vtkTriangle::ComputeNormalDirection(pt1, pt2, pt3, n);
93
94 // Project point to plane
95 //
96 vtkPlane::GeneralizedProjectPoint(x,pt1,n,cp);
97
98 // Construct matrices. Since we have over determined system, need to find
99 // which 2 out of 3 equations to use to develop equations. (Any 2 should
100 // work since we've projected point to plane.)
101 //
102 for (maxComponent=0.0, i=0; i<3; i++)
103 {
104 // trying to avoid an expensive call to fabs()
105 if (n[i] < 0)
106 {
107 fabsn = -n[i];
108 }
109 else
110 {
111 fabsn = n[i];
112 }
113 if (fabsn > maxComponent)
114 {
115 maxComponent = fabsn;
116 idx = i;
117 }
118 }
119 for (j=0, i=0; i<3; i++)
120 {
121 if ( i != idx )
122 {
123 indices[j++] = i;
124 }
125 }
126
127 for (i=0; i<2; i++)
128 {
129 rhs[i] = cp[indices[i]] - pt3[indices[i]];
130 c1[i] = pt1[indices[i]] - pt3[indices[i]];
131 c2[i] = pt2[indices[i]] - pt3[indices[i]];
132 }
133
134 if ( (det = vtkMath::Determinant2x2(c1,c2)) == 0.0 )
135 {
136 pcoords[0] = pcoords[1] = pcoords[2] = 0.0;
137 return -1;
138 }
139
140 pcoords[0] = vtkMath::Determinant2x2(rhs,c2) / det;
141 pcoords[1] = vtkMath::Determinant2x2(c1,rhs) / det;
142 pcoords[2] = 1.0 - (pcoords[0] + pcoords[1]);
143
144 // Okay, now find closest point to element
145 //
146 weights[0] = pcoords[2];
147 weights[1] = pcoords[0];
148 weights[2] = pcoords[1];
149
150 if ( pcoords[0] >= 0.0 && pcoords[0] <= 1.0 &&
151 pcoords[1] >= 0.0 && pcoords[1] <= 1.0 &&
152 pcoords[2] >= 0.0 && pcoords[2] <= 1.0 )
153 {
154 //projection distance
155 if (closestPoint)
156 {
157 dist2 = vtkMath::Distance2BetweenPoints(cp,x);
158 closestPoint[0] = cp[0];
159 closestPoint[1] = cp[1];
160 closestPoint[2] = cp[2];
161 }
162 return 1;
163 }
164 else
165 {
166 double t;
167 if (closestPoint)
168 {
169 if ( pcoords[0] < 0.0 && pcoords[1] < 0.0 )
170 {
171 dist2Point = vtkMath::Distance2BetweenPoints(x,pt3);
172 dist2Line1 = vtkLine::DistanceToLine(x,pt1,pt3,t,closestPoint1);
173 dist2Line2 = vtkLine::DistanceToLine(x,pt3,pt2,t,closestPoint2);
174 if (dist2Point < dist2Line1)
175 {
176 dist2 = dist2Point;
177 closest = pt3;
178 }
179 else
180 {
181 dist2 = dist2Line1;
182 closest = closestPoint1;
183 }
184 if (dist2Line2 < dist2)
185 {
186 dist2 = dist2Line2;
187 closest = closestPoint2;
188 }
189 for (i=0; i<3; i++)
190 {
191 closestPoint[i] = closest[i];
192 }
193 }
194 else if ( pcoords[1] < 0.0 && pcoords[2] < 0.0 )
195 {
196 dist2Point = vtkMath::Distance2BetweenPoints(x,pt1);
197 dist2Line1 = vtkLine::DistanceToLine(x,pt1,pt3,t,closestPoint1);
198 dist2Line2 = vtkLine::DistanceToLine(x,pt1,pt2,t,closestPoint2);
199 if (dist2Point < dist2Line1)
200 {
201 dist2 = dist2Point;
202 closest = pt1;
203 }
204 else
205 {
206 dist2 = dist2Line1;
207 closest = closestPoint1;
208 }
209 if (dist2Line2 < dist2)
210 {
211 dist2 = dist2Line2;
212 closest = closestPoint2;
213 }
214 for (i=0; i<3; i++)
215 {
216 closestPoint[i] = closest[i];
217 }
218 }
219 else if ( pcoords[0] < 0.0 && pcoords[2] < 0.0 )
220 {
221 dist2Point = vtkMath::Distance2BetweenPoints(x,pt2);
222 dist2Line1 = vtkLine::DistanceToLine(x,pt2,pt3,t,closestPoint1);
223 dist2Line2 = vtkLine::DistanceToLine(x,pt1,pt2,t,closestPoint2);
224 if (dist2Point < dist2Line1)
225 {
226 dist2 = dist2Point;
227 closest = pt2;
228 }
229 else
230 {
231 dist2 = dist2Line1;
232 closest = closestPoint1;
233 }
234 if (dist2Line2 < dist2)
235 {
236 dist2 = dist2Line2;
237 closest = closestPoint2;
238 }
239 for (i=0; i<3; i++)
240 {
241 closestPoint[i] = closest[i];
242 }
243 }
244 else if ( pcoords[0] < 0.0 )
245 {
246 dist2 = vtkLine::DistanceToLine(x,pt2,pt3,t,closestPoint);
247 }
248 else if ( pcoords[1] < 0.0 )
249 {
250 dist2 = vtkLine::DistanceToLine(x,pt1,pt3,t,closestPoint);
251 }
252 else if ( pcoords[2] < 0.0 )
253 {
254 dist2 = vtkLine::DistanceToLine(x,pt1,pt2,t,closestPoint);
255 }
256 }
257 return 0;
258 }
259 }
260
261 //----------------------------------------------------------------------------
EvaluateLocation(int & vtkNotUsed (subId),double pcoords[3],double x[3],double * weights)262 void vtkTriangle::EvaluateLocation(int& vtkNotUsed(subId), double pcoords[3],
263 double x[3], double *weights)
264 {
265 double u3;
266 double pt0[3], pt1[3], pt2[3];
267 int i;
268
269 this->Points->GetPoint(0, pt0);
270 this->Points->GetPoint(1, pt1);
271 this->Points->GetPoint(2, pt2);
272
273 u3 = 1.0 - pcoords[0] - pcoords[1];
274
275 for (i=0; i<3; i++)
276 {
277 x[i] = pt0[i]*u3 + pt1[i]*pcoords[0] + pt2[i]*pcoords[1];
278 }
279
280 weights[0] = u3;
281 weights[1] = pcoords[0];
282 weights[2] = pcoords[1];
283 }
284
285 //----------------------------------------------------------------------------
286 // Compute iso-parametric interpolation functions
287 //
InterpolationFunctions(double pcoords[3],double sf[3])288 void vtkTriangle::InterpolationFunctions(double pcoords[3], double sf[3])
289 {
290 sf[0] = 1. - pcoords[0] - pcoords[1];
291 sf[1] = pcoords[0];
292 sf[2] = pcoords[1];
293 }
294
295 //----------------------------------------------------------------------------
InterpolationDerivs(double *,double derivs[6])296 void vtkTriangle::InterpolationDerivs(double *, double derivs[6])
297 {
298 //r-derivatives
299 derivs[0] = -1;
300 derivs[1] = 1;
301 derivs[2] = 0;
302
303 //s-derivatives
304 derivs[3] = -1;
305 derivs[4] = 0;
306 derivs[5] = 1;
307 }
308
309 //----------------------------------------------------------------------------
CellBoundary(int vtkNotUsed (subId),double pcoords[3],vtkIdList * pts)310 int vtkTriangle::CellBoundary(int vtkNotUsed(subId), double pcoords[3],
311 vtkIdList *pts)
312 {
313 double t1=pcoords[0]-pcoords[1];
314 double t2=0.5*(1.0-pcoords[0])-pcoords[1];
315 double t3=2.0*pcoords[0]+pcoords[1]-1.0;
316
317 pts->SetNumberOfIds(2);
318
319 // compare against three lines in parametric space that divide element
320 // into three pieces
321 if ( t1 >= 0.0 && t2 >= 0.0 )
322 {
323 pts->SetId(0,this->PointIds->GetId(0));
324 pts->SetId(1,this->PointIds->GetId(1));
325 }
326
327 else if ( t2 < 0.0 && t3 >= 0.0 )
328 {
329 pts->SetId(0,this->PointIds->GetId(1));
330 pts->SetId(1,this->PointIds->GetId(2));
331 }
332
333 else //( t1 < 0.0 && t3 < 0.0 )
334 {
335 pts->SetId(0,this->PointIds->GetId(2));
336 pts->SetId(1,this->PointIds->GetId(0));
337 }
338
339 if ( pcoords[0] < 0.0 || pcoords[1] < 0.0 ||
340 pcoords[0] > 1.0 || pcoords[1] > 1.0 ||
341 (1.0 - pcoords[0] - pcoords[1]) < 0.0 )
342 {
343 return 0;
344 }
345 else
346 {
347 return 1;
348 }
349
350 }
351
352 //----------------------------------------------------------------------------
353 //
354 // Marching triangles
355 //
356 typedef int EDGE_LIST;
357 typedef struct {
358 EDGE_LIST edges[3];
359 } LINE_CASES;
360
361 static LINE_CASES lineCases[] = {
362 {{-1, -1, -1}},
363 {{0, 2, -1}},
364 {{1, 0, -1}},
365 {{1, 2, -1}},
366 {{2, 1, -1}},
367 {{0, 1, -1}},
368 {{2, 0, -1}},
369 {{-1, -1, -1}}
370 };
371
372 static int edges[3][2] = { {0,1}, {1,2}, {2,0} };
373
374 //----------------------------------------------------------------------------
GetEdgeArray(int edgeId)375 int *vtkTriangle::GetEdgeArray(int edgeId)
376 {
377 return edges[edgeId];
378 }
379
380 //----------------------------------------------------------------------------
Contour(double value,vtkDataArray * cellScalars,vtkIncrementalPointLocator * locator,vtkCellArray * verts,vtkCellArray * lines,vtkCellArray * vtkNotUsed (polys),vtkPointData * inPd,vtkPointData * outPd,vtkCellData * inCd,vtkIdType cellId,vtkCellData * outCd)381 void vtkTriangle::Contour(double value, vtkDataArray *cellScalars,
382 vtkIncrementalPointLocator *locator,
383 vtkCellArray *verts,
384 vtkCellArray *lines,
385 vtkCellArray *vtkNotUsed(polys),
386 vtkPointData *inPd, vtkPointData *outPd,
387 vtkCellData *inCd, vtkIdType cellId,
388 vtkCellData *outCd)
389 {
390 static int CASE_MASK[3] = {1,2,4};
391 LINE_CASES *lineCase;
392 EDGE_LIST *edge;
393 int i, j, index, *vert;
394 vtkIdType pts[2];
395 int e1, e2, newCellId;
396 double t, x1[3], x2[3], x[3], deltaScalar;
397 vtkIdType offset = verts->GetNumberOfCells();
398
399 // Build the case table
400 for ( i=0, index = 0; i < 3; i++)
401 {
402 if (cellScalars->GetComponent(i,0) >= value)
403 {
404 index |= CASE_MASK[i];
405 }
406 }
407
408 lineCase = lineCases + index;
409 edge = lineCase->edges;
410
411 for ( ; edge[0] > -1; edge += 2 )
412 {
413 for (i=0; i<2; i++) // insert line
414 {
415 vert = edges[edge[i]];
416 // calculate a preferred interpolation direction
417 deltaScalar = (cellScalars->GetComponent(vert[1],0)
418 - cellScalars->GetComponent(vert[0],0));
419 if (deltaScalar > 0)
420 {
421 e1 = vert[0]; e2 = vert[1];
422 }
423 else
424 {
425 e1 = vert[1]; e2 = vert[0];
426 deltaScalar = -deltaScalar;
427 }
428
429 // linear interpolation
430 if (deltaScalar == 0.0)
431 {
432 t = 0.0;
433 }
434 else
435 {
436 t = (value - cellScalars->GetComponent(e1,0)) / deltaScalar;
437 }
438
439 this->Points->GetPoint(e1, x1);
440 this->Points->GetPoint(e2, x2);
441
442 for (j=0; j<3; j++)
443 {
444 x[j] = x1[j] + t * (x2[j] - x1[j]);
445 }
446 if ( locator->InsertUniquePoint(x, pts[i]) )
447 {
448 if ( outPd )
449 {
450 vtkIdType p1 = this->PointIds->GetId(e1);
451 vtkIdType p2 = this->PointIds->GetId(e2);
452 outPd->InterpolateEdge(inPd,pts[i],p1,p2,t);
453 }
454 }
455 }
456 // check for degenerate line
457 if ( pts[0] != pts[1] )
458 {
459 newCellId = offset + lines->InsertNextCell(2,pts);
460 outCd->CopyData(inCd,cellId,newCellId);
461 }
462 }
463 }
464
465 //----------------------------------------------------------------------------
466 // Get the edge specified by edgeId (range 0 to 2) and return that edge's
467 // coordinates.
GetEdge(int edgeId)468 vtkCell *vtkTriangle::GetEdge(int edgeId)
469 {
470 int edgeIdPlus1 = (edgeId > 1 ? 0 : (edgeId+1) );
471
472 // load point id's
473 this->Line->PointIds->SetId(0,this->PointIds->GetId(edgeId));
474 this->Line->PointIds->SetId(1,this->PointIds->GetId(edgeIdPlus1));
475
476 // load coordinates
477 this->Line->Points->SetPoint(0,this->Points->GetPoint(edgeId));
478 this->Line->Points->SetPoint(1,this->Points->GetPoint(edgeIdPlus1));
479
480 return this->Line;
481 }
482
483 //----------------------------------------------------------------------------
484 // Plane intersection plus in/out test on triangle. The in/out test is
485 // performed using tol as the tolerance.
IntersectWithLine(double p1[3],double p2[3],double tol,double & t,double x[3],double pcoords[3],int & subId)486 int vtkTriangle::IntersectWithLine(double p1[3], double p2[3], double tol,
487 double& t, double x[3], double pcoords[3],
488 int& subId)
489 {
490 double pt1[3], pt2[3], pt3[3], n[3];
491 double tol2 = tol*tol;
492 double closestPoint[3];
493 double dist2, weights[3];
494
495 subId = 0;
496
497 // Get normal for triangle
498 //
499 this->Points->GetPoint(1, pt1);
500 this->Points->GetPoint(2, pt2);
501 this->Points->GetPoint(0, pt3);
502
503 vtkTriangle::ComputeNormal (pt1, pt2, pt3, n);
504
505 // Intersect plane of triangle with line
506 //
507 if ( ! vtkPlane::IntersectWithLine(p1,p2,n,pt1,t,x) )
508 {
509 pcoords[0] = pcoords[1] = pcoords[2] = 0.0;
510 return 0;
511 }
512
513 // Evaluate position
514 //
515 int inside;
516 if ( (inside = this->EvaluatePosition(x, closestPoint, subId, pcoords,
517 dist2, weights)) >= 0)
518 {
519 if ( dist2 <= tol2 )
520 {
521 pcoords[2] = 0.0;
522 return 1;
523 }
524 return inside;
525 }
526
527 // so the easy test failed. The line is not intersecting the triangle.
528 // Let's now do the 3d case check to see how close the line comes.
529 // basically we just need to test against the three lines of the triangle
530 this->Line->PointIds->InsertId(0,0);
531 this->Line->PointIds->InsertId(1,1);
532
533 if (pcoords[2] < 0.0)
534 {
535 this->Line->Points->InsertPoint(0,pt1);
536 this->Line->Points->InsertPoint(1,pt2);
537 if (this->Line->IntersectWithLine(p1,p2,tol,t,x,pcoords,subId))
538 {
539 pcoords[2] = 0.0;
540 return 1;
541 }
542 }
543
544 if (pcoords[0] < 0.0)
545 {
546 this->Line->Points->InsertPoint(0,pt2);
547 this->Line->Points->InsertPoint(1,pt3);
548 if (this->Line->IntersectWithLine(p1,p2,tol,t,x,pcoords,subId))
549 {
550 pcoords[2] = 0.0;
551 return 1;
552 }
553 }
554
555 if (pcoords[1] < 0.0)
556 {
557 this->Line->Points->InsertPoint(0,pt3);
558 this->Line->Points->InsertPoint(1,pt1);
559 if (this->Line->IntersectWithLine(p1,p2,tol,t,x,pcoords,subId))
560 {
561 pcoords[2] = 0.0;
562 return 1;
563 }
564 }
565
566 pcoords[0] = pcoords[1] = pcoords[2] = 0.0;
567 return 0;
568 }
569
570 //----------------------------------------------------------------------------
Triangulate(int vtkNotUsed (index),vtkIdList * ptIds,vtkPoints * pts)571 int vtkTriangle::Triangulate(int vtkNotUsed(index), vtkIdList *ptIds,
572 vtkPoints *pts)
573 {
574 pts->Reset();
575 ptIds->Reset();
576
577 for ( int i=0; i < 3; i++ )
578 {
579 ptIds->InsertId(i,this->PointIds->GetId(i));
580 pts->InsertPoint(i,this->Points->GetPoint(i));
581 }
582
583 return 1;
584 }
585
586 //----------------------------------------------------------------------------
587 // Used a staged computation: first compute derivatives in local x'-y'
588 // coordinate system; then convert into x-y-z modelling system.
Derivatives(int vtkNotUsed (subId),double vtkNotUsed (pcoords)[3],double * values,int dim,double * derivs)589 void vtkTriangle::Derivatives(int vtkNotUsed(subId), double vtkNotUsed(pcoords)[3],
590 double *values, int dim, double *derivs)
591 {
592 double v0[2], v1[2], v2[2], v[3], v10[3], v20[3], lenX;
593 double x0[3], x1[3], x2[3], n[3];
594 double *J[2], J0[2], J1[2];
595 double *JI[2], JI0[2], JI1[2];
596 double functionDerivs[6], sum[2], dBydx, dBydy;
597 int i, j;
598
599 // Project points of triangle into 2D system
600 this->Points->GetPoint(0, x0);
601 this->Points->GetPoint(1, x1);
602 this->Points->GetPoint(2, x2);
603 vtkTriangle::ComputeNormal (x0, x1, x2, n);
604
605 for (i=0; i < 3; i++)
606 {
607 v10[i] = x1[i] - x0[i];
608 v[i] = x2[i] - x0[i];
609 }
610
611 vtkMath::Cross(n,v10,v20); //creates local y' axis
612
613 if ( (lenX=vtkMath::Normalize(v10)) <= 0.0
614 || vtkMath::Normalize(v20) <= 0.0 ) //degenerate
615 {
616 for ( j=0; j < dim; j++ )
617 {
618 for ( i=0; i < 3; i++ )
619 {
620 derivs[j*dim + i] = 0.0;
621 }
622 }
623 return;
624 }
625
626 v0[0] = v0[1] = 0.0; //convert points to 2D (i.e., local system)
627 v1[0] = lenX; v1[1] = 0.0;
628 v2[0] = vtkMath::Dot(v,v10);
629 v2[1] = vtkMath::Dot(v,v20);
630
631 // Compute interpolation function derivatives
632 vtkTriangle::InterpolationDerivs(NULL,functionDerivs);
633
634 // Compute Jacobian: Jacobian is constant for a triangle.
635 J[0] = J0; J[1] = J1;
636 JI[0] = JI0; JI[1] = JI1;
637
638 J[0][0] = v1[0] - v0[0];
639 J[1][0] = v2[0] - v0[0];
640 J[0][1] = v1[1] - v0[1];
641 J[1][1] = v2[1] - v0[1];
642
643 // Compute inverse Jacobian
644 vtkMath::InvertMatrix(J,JI,2);
645
646 // Loop over "dim" derivative values. For each set of values, compute
647 // derivatives in local system and then transform into modelling system.
648 // First compute derivatives in local x'-y' coordinate system
649 for ( j=0; j < dim; j++ )
650 {
651 sum[0] = sum[1] = 0.0;
652 for ( i=0; i < 3; i++) //loop over interp. function derivatives
653 {
654 sum[0] += functionDerivs[i] * values[dim*i + j];
655 sum[1] += functionDerivs[3 + i] * values[dim*i + j];
656 }
657 dBydx = sum[0]*JI[0][0] + sum[1]*JI[0][1];
658 dBydy = sum[0]*JI[1][0] + sum[1]*JI[1][1];
659
660 // Transform into global system (dot product with global axes)
661 derivs[3*j] = dBydx * v10[0] + dBydy * v20[0];
662 derivs[3*j + 1] = dBydx * v10[1] + dBydy * v20[1];
663 derivs[3*j + 2] = dBydx * v10[2] + dBydy * v20[2];
664 }
665 }
666
667 //----------------------------------------------------------------------------
668 // Compute the triangle normal from a points list, and a list of point ids
669 // that index into the points list.
ComputeNormal(vtkPoints * p,int vtkNotUsed (numPts),vtkIdType * pts,double n[3])670 void vtkTriangle::ComputeNormal(vtkPoints *p, int vtkNotUsed(numPts),
671 vtkIdType *pts, double n[3])
672 {
673 double v1[3], v2[3], v3[3];
674
675 p->GetPoint(pts[0],v1);
676 p->GetPoint(pts[1],v2);
677 p->GetPoint(pts[2],v3);
678
679 vtkTriangle::ComputeNormal(v1,v2,v3,n);
680 }
681
682 //----------------------------------------------------------------------------
683 // Compute the circumcenter (center[3]) and radius squared (method
684 // return value) of a triangle defined by the three points x1, x2, and
685 // x3. (Note that the coordinates are 2D. 3D points can be used but
686 // the z-component will be ignored.)
Circumcircle(double x1[2],double x2[2],double x3[2],double center[2])687 double vtkTriangle::Circumcircle(double x1[2], double x2[2], double x3[2],
688 double center[2])
689 {
690 double n12[2], n13[2], x12[2], x13[2];
691 double *A[2], rhs[2], sum, diff;
692 int i;
693
694 // calculate normals and intersection points of bisecting planes.
695 //
696 for (i=0; i<2; i++)
697 {
698 n12[i] = x2[i] - x1[i];
699 n13[i] = x3[i] - x1[i];
700 x12[i] = (x2[i] + x1[i])/2.0;
701 x13[i] = (x3[i] + x1[i])/2.0;
702 }
703
704 // Compute solutions to the intersection of two bisecting lines
705 // (2-eqns. in 2-unknowns).
706 //
707 // form system matrices
708 //
709 A[0] = n12;
710 A[1] = n13;
711
712 rhs[0] = vtkMath::Dot2D(n12,x12);
713 rhs[1] = vtkMath::Dot2D(n13,x13);
714
715 // Solve system of equations
716 //
717 if ( vtkMath::SolveLinearSystem(A,rhs,2) == 0 )
718 {
719 center[0] = center[1] = 0.0;
720 return VTK_DOUBLE_MAX;
721 }
722 else
723 {
724 center[0] = rhs[0]; center[1] = rhs[1];
725 }
726
727 //determine average value of radius squared
728 for (sum=0, i=0; i<2; i++)
729 {
730 diff = x1[i] - center[i];
731 sum += diff*diff;
732 diff = x2[i] - center[i];
733 sum += diff*diff;
734 diff = x3[i] - center[i];
735 sum += diff*diff;
736 }
737
738 if ( (sum /= 3.0) > VTK_DOUBLE_MAX )
739 {
740 return VTK_DOUBLE_MAX;
741 }
742 else
743 {
744 return sum;
745 }
746 }
747
748 //----------------------------------------------------------------------------
749 // Given a 2D point x[2], determine the barycentric coordinates of the point.
750 // Barycentric coordinates are a natural coordinate system for simplices that
751 // express a position as a linear combination of the vertices. For a
752 // triangle, there are three barycentric coordinates (because there are
753 // fourthree vertices), and the sum of the coordinates must equal 1. If a
754 // point x is inside a simplex, then all three coordinates will be strictly
755 // positive. If two coordinates are zero (so the third =1), then the
756 // point x is on a vertex. If one coordinates are zero, the point x is on an
757 // edge. In this method, you must specify the vertex coordinates x1->x3.
758 // Returns 0 if triangle is degenerate.
BarycentricCoords(double x[2],double x1[2],double x2[2],double x3[2],double bcoords[3])759 int vtkTriangle::BarycentricCoords(double x[2], double x1[2], double x2[2],
760 double x3[2], double bcoords[3])
761 {
762 double *A[3], p[3], a1[3], a2[3], a3[3];
763 int i;
764
765 // Homogenize the variables; load into arrays.
766 //
767 a1[0] = x1[0]; a1[1] = x2[0]; a1[2] = x3[0];
768 a2[0] = x1[1]; a2[1] = x2[1]; a2[2] = x3[1];
769 a3[0] = 1.0; a3[1] = 1.0; a3[2] = 1.0;
770 p[0] = x[0]; p[1] = x[1]; p[2] = 1.0;
771
772 // Now solve system of equations for barycentric coordinates
773 //
774 A[0] = a1;
775 A[1] = a2;
776 A[2] = a3;
777
778 if ( vtkMath::SolveLinearSystem(A,p,3) )
779 {
780 for (i=0; i<3; i++)
781 {
782 bcoords[i] = p[i];
783 }
784 return 1;
785 }
786 else
787 {
788 return 0;
789 }
790 }
791
792 //----------------------------------------------------------------------------
793 // Project triangle defined in 3D to 2D coordinates. Returns 0 if degenerate
794 // triangle; non-zero value otherwise. Input points are x1->x3; output 2D
795 // points are v1->v3.
ProjectTo2D(double x1[3],double x2[3],double x3[3],double v1[2],double v2[2],double v3[2])796 int vtkTriangle::ProjectTo2D(double x1[3], double x2[3], double x3[3],
797 double v1[2], double v2[2], double v3[2])
798 {
799 double n[3], v21[3], v31[3], v[3], xLen;
800
801 // Get normal for triangle
802 vtkTriangle::ComputeNormal (x1, x2, x3, n);
803
804 for (int i=0; i < 3; i++)
805 {
806 v21[i] = x2[i] - x1[i];
807 v31[i] = x3[i] - x1[i];
808 }
809
810 if ( (xLen=vtkMath::Normalize(v21)) <= 0.0 )
811 {
812 return 0;
813 }
814
815 // The first point is at (0,0); the next at (xLen,0); compute the other
816 // point relative to the first two.
817 v1[0] = v1[1] = 0.0;
818 v2[0] = xLen; v2[1] = 0.0;
819
820 vtkMath::Cross(n,v21,v);
821
822 v3[0] = vtkMath::Dot(v31,v21);
823 v3[1] = vtkMath::Dot(v31,v);
824
825 return 1;
826 }
827
828 //----------------------------------------------------------------------------
829 // Support triangle clipping. Note that the table defines triangles (three ids
830 // at a time define a triangle, -1 ends the list). Numbers in the list >= 100
831 // correspond to already existing vertices; otherwise the numbers refer to edge
832 // ids.
833 typedef int TRIANGLE_EDGE_LIST;
834 typedef struct {
835 TRIANGLE_EDGE_LIST edges[7];
836 } TRIANGLE_CASES;
837
838 static TRIANGLE_CASES triangleCases[] = {
839 {{-1, -1, -1, -1, -1, -1, -1}}, // 0
840 {{0, 2, 100, -1, -1, -1, -1}}, // 1
841 {{1, 0, 101, -1, -1, -1, -1}}, // 2
842 {{1, 2, 100, 1, 100, 101, -1}}, // 3
843 {{2, 1, 102, -1, -1, -1, -1}}, // 4
844 {{0, 1, 102, 102, 100, 0, -1}}, // 5
845 {{0, 101, 2, 2, 101, 102, -1}}, // 6
846 {{100, 101, 102, -1, -1, -1, -1}} // 7
847 };
848
849 //----------------------------------------------------------------------------
850 // Clip this triangle using scalar value provided. Like contouring, except
851 // that it cuts the triangle to produce other triangles.
Clip(double value,vtkDataArray * cellScalars,vtkIncrementalPointLocator * locator,vtkCellArray * tris,vtkPointData * inPd,vtkPointData * outPd,vtkCellData * inCd,vtkIdType cellId,vtkCellData * outCd,int insideOut)852 void vtkTriangle::Clip(double value, vtkDataArray *cellScalars,
853 vtkIncrementalPointLocator *locator, vtkCellArray *tris,
854 vtkPointData *inPd, vtkPointData *outPd,
855 vtkCellData *inCd, vtkIdType cellId, vtkCellData *outCd,
856 int insideOut)
857 {
858 static int CASE_MASK[3] = {1,2,4};
859 TRIANGLE_CASES *triangleCase;
860 TRIANGLE_EDGE_LIST *edge;
861 int i, j, index, *vert;
862 int e1, e2, newCellId;
863 vtkIdType pts[3];
864 int vertexId;
865 double t, x1[3], x2[3], x[3], deltaScalar;
866
867 // Build the case table
868 if ( insideOut )
869 {
870 for ( i=0, index = 0; i < 3; i++)
871 {
872 if (cellScalars->GetComponent(i,0) <= value)
873 {
874 index |= CASE_MASK[i];
875 }
876 }
877 }
878 else
879 {
880 for ( i=0, index = 0; i < 3; i++)
881 {
882 if (cellScalars->GetComponent(i,0) > value)
883 {
884 index |= CASE_MASK[i];
885 }
886 }
887 }
888
889 // Select the case based on the index and get the list of edges for this case
890 triangleCase = triangleCases + index;
891 edge = triangleCase->edges;
892
893 // generate each triangle
894 for ( ; edge[0] > -1; edge += 3 )
895 {
896 for (i=0; i<3; i++) // insert triangle
897 {
898 // vertex exists, and need not be interpolated
899 if (edge[i] >= 100)
900 {
901 vertexId = edge[i] - 100;
902 this->Points->GetPoint(vertexId, x);
903 if ( locator->InsertUniquePoint(x, pts[i]) )
904 {
905 outPd->CopyData(inPd,this->PointIds->GetId(vertexId),pts[i]);
906 }
907 }
908
909 else //new vertex, interpolate
910 {
911 vert = edges[edge[i]];
912
913 // calculate a preferred interpolation direction
914 deltaScalar = (cellScalars->GetComponent(vert[1],0) -
915 cellScalars->GetComponent(vert[0],0));
916 if (deltaScalar > 0)
917 {
918 e1 = vert[0]; e2 = vert[1];
919 }
920 else
921 {
922 e1 = vert[1]; e2 = vert[0];
923 deltaScalar = -deltaScalar;
924 }
925
926 // linear interpolation
927 if (deltaScalar == 0.0)
928 {
929 t = 0.0;
930 }
931 else
932 {
933 t = (value - cellScalars->GetComponent(e1,0)) / deltaScalar;
934 }
935
936 this->Points->GetPoint(e1, x1);
937 this->Points->GetPoint(e2, x2);
938
939 for (j=0; j<3; j++)
940 {
941 x[j] = x1[j] + t * (x2[j] - x1[j]);
942 }
943 if ( locator->InsertUniquePoint(x, pts[i]) )
944 {
945 vtkIdType p1 = this->PointIds->GetId(e1);
946 vtkIdType p2 = this->PointIds->GetId(e2);
947 outPd->InterpolateEdge(inPd,pts[i],p1,p2,t);
948 }
949 }
950 }
951 // check for degenerate tri's
952 if (pts[0] == pts[1] || pts[0] == pts[2] || pts[1] == pts[2])
953 {
954 continue;
955 }
956
957 newCellId = tris->InsertNextCell(3,pts);
958 outCd->CopyData(inCd,cellId,newCellId);
959 }
960 }
961
962 //----------------------------------------------------------------------------
963 // Given a point x, determine whether it is inside (within the
964 // tolerance squared, tol2) the triangle defined by the three
965 // coordinate values p1, p2, p3. Method is via comparing dot products.
966 // (Note: in current implementation the tolerance only works in the
967 // neighborhood of the three vertices of the triangle.
PointInTriangle(double x[3],double p1[3],double p2[3],double p3[3],double tol2)968 int vtkTriangle::PointInTriangle(double x[3], double p1[3], double p2[3],
969 double p3[3], double tol2)
970 {
971 double x1[3], x2[3], x3[3], v13[3], v21[3], v32[3];
972 double n1[3], n2[3], n3[3];
973 int i;
974
975 // Compute appropriate vectors
976 //
977 for (i=0; i<3; i++)
978 {
979 x1[i] = x[i] - p1[i];
980 x2[i] = x[i] - p2[i];
981 x3[i] = x[i] - p3[i];
982 v13[i] = p1[i] - p3[i];
983 v21[i] = p2[i] - p1[i];
984 v32[i] = p3[i] - p2[i];
985 }
986
987 // See whether intersection point is within tolerance of a vertex.
988 //
989 if ( (x1[0]*x1[0] + x1[1]*x1[1] + x1[2]*x1[2]) <= tol2 ||
990 (x2[0]*x2[0] + x2[1]*x2[1] + x2[2]*x2[2]) <= tol2 ||
991 (x3[0]*x3[0] + x3[1]*x3[1] + x3[2]*x3[2]) <= tol2 )
992 {
993 return 1;
994 }
995
996 // If not near a vertex, check whether point is inside of triangular face.
997 //
998 // Obtain normal off of triangular face
999 //
1000 vtkMath::Cross (x1, v13, n1);
1001 vtkMath::Cross (x2, v21, n2);
1002 vtkMath::Cross (x3, v32, n3);
1003
1004 // Check whether ALL the three normals go in same direction
1005 //
1006 if ( (vtkMath::Dot(n1,n2) >= 0.0) &&
1007 (vtkMath::Dot(n2,n3) >= 0.0) &&
1008 (vtkMath::Dot(n1,n3) >= 0.0) )
1009 {
1010 return 1;
1011 }
1012 else
1013 {
1014 return 0;
1015 }
1016 }
1017
1018 //----------------------------------------------------------------------------
GetParametricDistance(double pcoords[3])1019 double vtkTriangle::GetParametricDistance(double pcoords[3])
1020 {
1021 int i;
1022 double pDist, pDistMax=0.0;
1023 double pc[3];
1024
1025 pc[0] = pcoords[0];
1026 pc[1] = pcoords[1];
1027 pc[2] = 1.0 - pcoords[0] - pcoords[1];
1028
1029 for (i=0; i<3; i++)
1030 {
1031 if ( pc[i] < 0.0 )
1032 {
1033 pDist = -pc[i];
1034 }
1035 else if ( pc[i] > 1.0 )
1036 {
1037 pDist = pc[i] - 1.0;
1038 }
1039 else //inside the cell in the parametric direction
1040 {
1041 pDist = 0.0;
1042 }
1043 if ( pDist > pDistMax )
1044 {
1045 pDistMax = pDist;
1046 }
1047 }
1048
1049 return pDistMax;
1050 }
1051
1052
1053 //----------------------------------------------------------------------------
ComputeQuadric(double x1[3],double x2[3],double x3[3],double quadric[4][4])1054 void vtkTriangle::ComputeQuadric(double x1[3], double x2[3], double x3[3],
1055 double quadric[4][4])
1056 {
1057 double crossX1X2[3], crossX2X3[3], crossX3X1[3];
1058 double determinantABC;
1059 double ABCx[3][3];
1060 double n[4];
1061 int i, j;
1062
1063 for (i = 0; i < 3; i++)
1064 {
1065 ABCx[0][i] = x1[i];
1066 ABCx[1][i] = x2[i];
1067 ABCx[2][i] = x3[i];
1068 }
1069
1070 vtkMath::Cross(x1, x2, crossX1X2);
1071 vtkMath::Cross(x2, x3, crossX2X3);
1072 vtkMath::Cross(x3, x1, crossX3X1);
1073 determinantABC = vtkMath::Determinant3x3(ABCx);
1074
1075 n[0] = crossX1X2[0] + crossX2X3[0] + crossX3X1[0];
1076 n[1] = crossX1X2[1] + crossX2X3[1] + crossX3X1[1];
1077 n[2] = crossX1X2[2] + crossX2X3[2] + crossX3X1[2];
1078 n[3] = -determinantABC;
1079
1080 for (i = 0; i < 4; i++)
1081 {
1082 for (j = 0; j < 4; j++)
1083 {
1084 quadric[i][j] = n[i] * n[j];
1085 }
1086 }
1087 }
1088
1089 //----------------------------------------------------------------------------
ComputeQuadric(double x1[3],double x2[3],double x3[3],vtkQuadric * quadric)1090 void vtkTriangle::ComputeQuadric(double x1[3], double x2[3], double x3[3],
1091 vtkQuadric *quadric)
1092 {
1093 double quadricMatrix[4][4];
1094
1095 ComputeQuadric(x1, x2, x3, quadricMatrix);
1096 quadric->SetCoefficients(quadricMatrix[0][0], quadricMatrix[1][1],
1097 quadricMatrix[2][2], 2*quadricMatrix[0][1],
1098 2*quadricMatrix[1][2], 2*quadricMatrix[0][2],
1099 2*quadricMatrix[0][3], 2*quadricMatrix[1][3],
1100 2*quadricMatrix[2][3], quadricMatrix[3][3]);
1101 }
1102
1103 //----------------------------------------------------------------------------
1104 static double vtkTriangleCellPCoords[9] =
1105 {0.0,0.0,0.0, 1.0,0.0,0.0, 0.0,1.0,0.0};
GetParametricCoords()1106 double *vtkTriangle::GetParametricCoords()
1107 {
1108 return vtkTriangleCellPCoords;
1109 }
1110
1111 //----------------------------------------------------------------------------
PrintSelf(ostream & os,vtkIndent indent)1112 void vtkTriangle::PrintSelf(ostream& os, vtkIndent indent)
1113 {
1114 this->Superclass::PrintSelf(os,indent);
1115
1116 os << indent << "Line:\n";
1117 this->Line->PrintSelf(os,indent.GetNextIndent());
1118 }
1119