1 /*=========================================================================
2
3 Program: Visualization Toolkit
4 Module: vtkQuad.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 "vtkQuad.h"
16
17 #include "vtkCellArray.h"
18 #include "vtkCellData.h"
19 #include "vtkIncrementalPointLocator.h"
20 #include "vtkLine.h"
21 #include "vtkMath.h"
22 #include "vtkObjectFactory.h"
23 #include "vtkPlane.h"
24 #include "vtkPointData.h"
25 #include "vtkPoints.h"
26 #include "vtkTriangle.h"
27
28 vtkStandardNewMacro(vtkQuad);
29
30 static const double VTK_DIVERGED = 1.e6;
31
32 //------------------------------------------------------------------------------
33 // Construct the quad with four points.
vtkQuad()34 vtkQuad::vtkQuad()
35 {
36 this->Points->SetNumberOfPoints(4);
37 this->PointIds->SetNumberOfIds(4);
38 for (int i = 0; i < 4; i++)
39 {
40 this->Points->SetPoint(i, 0.0, 0.0, 0.0);
41 this->PointIds->SetId(i, 0);
42 }
43 this->Line = vtkLine::New();
44 this->Triangle = vtkTriangle::New();
45 }
46
47 //------------------------------------------------------------------------------
~vtkQuad()48 vtkQuad::~vtkQuad()
49 {
50 this->Line->Delete();
51 this->Triangle->Delete();
52 }
53
54 //------------------------------------------------------------------------------
55 static const int VTK_QUAD_MAX_ITERATION = 20;
56 static const double VTK_QUAD_CONVERGED = 1.e-04;
57
ComputeNormal(vtkQuad * self,double pt1[3],double pt2[3],double pt3[3],double n[3])58 inline static void ComputeNormal(
59 vtkQuad* self, double pt1[3], double pt2[3], double pt3[3], double n[3])
60 {
61 vtkTriangle::ComputeNormal(pt1, pt2, pt3, n);
62
63 // If first three points are co-linear, then use fourth point
64 //
65 double pt4[3];
66 if (n[0] == 0.0 && n[1] == 0.0 && n[2] == 0.0)
67 {
68 self->Points->GetPoint(3, pt4);
69 vtkTriangle::ComputeNormal(pt2, pt3, pt4, n);
70 }
71 }
72
73 //------------------------------------------------------------------------------
EvaluatePosition(const double x[3],double closestPoint[3],int & subId,double pcoords[3],double & dist2,double weights[])74 int vtkQuad::EvaluatePosition(const double x[3], double closestPoint[3], int& subId,
75 double pcoords[3], double& dist2, double weights[])
76 {
77 int i, j;
78 double pt1[3], pt2[3], pt3[3], pt[3], n[3];
79 double det;
80 double maxComponent;
81 int idx = 0, indices[2];
82 int iteration, converged;
83 double params[2];
84 double fcol[2], rcol[2], scol[2], cp[3];
85 double derivs[8];
86
87 subId = 0;
88 pcoords[0] = pcoords[1] = params[0] = params[1] = 0.5;
89 pcoords[2] = 0.0;
90
91 // Get normal for quadrilateral
92 //
93 this->Points->GetPoint(0, pt1);
94 this->Points->GetPoint(1, pt2);
95 this->Points->GetPoint(2, pt3);
96 ComputeNormal(this, pt1, pt2, pt3, n);
97
98 // Project point to plane
99 //
100 vtkPlane::ProjectPoint(x, pt1, n, cp);
101
102 // Construct matrices. Since we have over determined system, need to find
103 // which 2 out of 3 equations to use to develop equations. (Any 2 should
104 // work since we've projected point to plane.)
105 //
106 for (maxComponent = 0.0, i = 0; i < 3; i++)
107 {
108 if (fabs(n[i]) > maxComponent)
109 {
110 maxComponent = fabs(n[i]);
111 idx = i;
112 }
113 }
114 for (j = 0, i = 0; i < 3; i++)
115 {
116 if (i != idx)
117 {
118 indices[j++] = i;
119 }
120 }
121
122 // Use Newton's method to solve for parametric coordinates
123 //
124 for (iteration = converged = 0; !converged && (iteration < VTK_QUAD_MAX_ITERATION); iteration++)
125 {
126 // calculate element interpolation functions and derivatives
127 //
128 vtkQuad::InterpolationFunctions(pcoords, weights);
129 vtkQuad::InterpolationDerivs(pcoords, derivs);
130
131 // calculate newton functions
132 //
133 for (i = 0; i < 2; i++)
134 {
135 fcol[i] = rcol[i] = scol[i] = 0.0;
136 }
137 for (i = 0; i < 4; i++)
138 {
139 this->Points->GetPoint(i, pt);
140 for (j = 0; j < 2; j++)
141 {
142 fcol[j] += pt[indices[j]] * weights[i];
143 rcol[j] += pt[indices[j]] * derivs[i];
144 scol[j] += pt[indices[j]] * derivs[i + 4];
145 }
146 }
147
148 for (j = 0; j < 2; j++)
149 {
150 fcol[j] -= cp[indices[j]];
151 }
152
153 // compute determinants and generate improvements
154 //
155 if ((det = vtkMath::Determinant2x2(rcol, scol)) == 0.0)
156 {
157 return -1;
158 }
159
160 pcoords[0] = params[0] - vtkMath::Determinant2x2(fcol, scol) / det;
161 pcoords[1] = params[1] - vtkMath::Determinant2x2(rcol, fcol) / det;
162
163 // check for convergence
164 //
165 if (((fabs(pcoords[0] - params[0])) < VTK_QUAD_CONVERGED) &&
166 ((fabs(pcoords[1] - params[1])) < VTK_QUAD_CONVERGED))
167 {
168 converged = 1;
169 }
170 // Test for bad divergence (S.Hirschberg 11.12.2001)
171 else if ((fabs(pcoords[0]) > VTK_DIVERGED) || (fabs(pcoords[1]) > VTK_DIVERGED))
172 {
173 return -1;
174 }
175
176 // if not converged, repeat
177 //
178 else
179 {
180 params[0] = pcoords[0];
181 params[1] = pcoords[1];
182 }
183 }
184
185 // if not converged, set the parametric coordinates to arbitrary values
186 // outside of element
187 //
188 if (!converged)
189 {
190 return -1;
191 }
192
193 vtkQuad::InterpolationFunctions(pcoords, weights);
194
195 if (pcoords[0] >= -0.001 && pcoords[0] <= 1.001 && pcoords[1] >= -0.001 && pcoords[1] <= 1.001)
196 {
197 if (closestPoint)
198 {
199 dist2 = vtkMath::Distance2BetweenPoints(cp, x); // projection distance
200 closestPoint[0] = cp[0];
201 closestPoint[1] = cp[1];
202 closestPoint[2] = cp[2];
203 }
204 return 1;
205 }
206 else
207 {
208 double t;
209 double pt4[3];
210
211 if (closestPoint)
212 {
213 this->Points->GetPoint(3, pt4);
214
215 if (pcoords[0] < 0.0 && pcoords[1] < 0.0)
216 {
217 dist2 = vtkMath::Distance2BetweenPoints(x, pt1);
218 for (i = 0; i < 3; i++)
219 {
220 closestPoint[i] = pt1[i];
221 }
222 }
223 else if (pcoords[0] > 1.0 && pcoords[1] < 0.0)
224 {
225 dist2 = vtkMath::Distance2BetweenPoints(x, pt2);
226 for (i = 0; i < 3; i++)
227 {
228 closestPoint[i] = pt2[i];
229 }
230 }
231 else if (pcoords[0] > 1.0 && pcoords[1] > 1.0)
232 {
233 dist2 = vtkMath::Distance2BetweenPoints(x, pt3);
234 for (i = 0; i < 3; i++)
235 {
236 closestPoint[i] = pt3[i];
237 }
238 }
239 else if (pcoords[0] < 0.0 && pcoords[1] > 1.0)
240 {
241 dist2 = vtkMath::Distance2BetweenPoints(x, pt4);
242 for (i = 0; i < 3; i++)
243 {
244 closestPoint[i] = pt4[i];
245 }
246 }
247 else if (pcoords[0] < 0.0)
248 {
249 dist2 = vtkLine::DistanceToLine(x, pt1, pt4, t, closestPoint);
250 }
251 else if (pcoords[0] > 1.0)
252 {
253 dist2 = vtkLine::DistanceToLine(x, pt2, pt3, t, closestPoint);
254 }
255 else if (pcoords[1] < 0.0)
256 {
257 dist2 = vtkLine::DistanceToLine(x, pt1, pt2, t, closestPoint);
258 }
259 else if (pcoords[1] > 1.0)
260 {
261 dist2 = vtkLine::DistanceToLine(x, pt3, pt4, t, closestPoint);
262 }
263 }
264 return 0;
265 }
266 }
267
268 //------------------------------------------------------------------------------
EvaluateLocation(int & vtkNotUsed (subId),const double pcoords[3],double x[3],double * weights)269 void vtkQuad::EvaluateLocation(
270 int& vtkNotUsed(subId), const double pcoords[3], double x[3], double* weights)
271 {
272 int i, j;
273 double pt[3];
274
275 vtkQuad::InterpolationFunctions(pcoords, weights);
276
277 x[0] = x[1] = x[2] = 0.0;
278 for (i = 0; i < 4; i++)
279 {
280 this->Points->GetPoint(i, pt);
281 for (j = 0; j < 3; j++)
282 {
283 x[j] += pt[j] * weights[i];
284 }
285 }
286 }
287
288 //------------------------------------------------------------------------------
289 // Compute iso-parametric interpolation functions
290 //
InterpolationFunctions(const double pcoords[3],double sf[4])291 void vtkQuad::InterpolationFunctions(const double pcoords[3], double sf[4])
292 {
293 double rm, sm;
294
295 rm = 1. - pcoords[0];
296 sm = 1. - pcoords[1];
297
298 sf[0] = rm * sm;
299 sf[1] = pcoords[0] * sm;
300 sf[2] = pcoords[0] * pcoords[1];
301 sf[3] = rm * pcoords[1];
302 }
303
304 //------------------------------------------------------------------------------
InterpolationDerivs(const double pcoords[3],double derivs[8])305 void vtkQuad::InterpolationDerivs(const double pcoords[3], double derivs[8])
306 {
307 double rm, sm;
308
309 rm = 1. - pcoords[0];
310 sm = 1. - pcoords[1];
311
312 derivs[0] = -sm;
313 derivs[1] = sm;
314 derivs[2] = pcoords[1];
315 derivs[3] = -pcoords[1];
316 derivs[4] = -rm;
317 derivs[5] = -pcoords[0];
318 derivs[6] = pcoords[0];
319 derivs[7] = rm;
320 }
321
322 //------------------------------------------------------------------------------
CellBoundary(int vtkNotUsed (subId),const double pcoords[3],vtkIdList * pts)323 int vtkQuad::CellBoundary(int vtkNotUsed(subId), const double pcoords[3], vtkIdList* pts)
324 {
325 double t1 = pcoords[0] - pcoords[1];
326 double t2 = 1.0 - pcoords[0] - pcoords[1];
327
328 pts->SetNumberOfIds(2);
329
330 // compare against two lines in parametric space that divide element
331 // into four pieces.
332 if (t1 >= 0.0 && t2 >= 0.0)
333 {
334 pts->SetId(0, this->PointIds->GetId(0));
335 pts->SetId(1, this->PointIds->GetId(1));
336 }
337
338 else if (t1 >= 0.0 && t2 < 0.0)
339 {
340 pts->SetId(0, this->PointIds->GetId(1));
341 pts->SetId(1, this->PointIds->GetId(2));
342 }
343
344 else if (t1 < 0.0 && t2 < 0.0)
345 {
346 pts->SetId(0, this->PointIds->GetId(2));
347 pts->SetId(1, this->PointIds->GetId(3));
348 }
349
350 else //( t1 < 0.0 && t2 >= 0.0 )
351 {
352 pts->SetId(0, this->PointIds->GetId(3));
353 pts->SetId(1, this->PointIds->GetId(0));
354 }
355
356 if (pcoords[0] < 0.0 || pcoords[0] > 1.0 || pcoords[1] < 0.0 || pcoords[1] > 1.0)
357 {
358 return 0;
359 }
360 else
361 {
362 return 1;
363 }
364 }
365
366 //------------------------------------------------------------------------------
367 // Marching (convex) quadrilaterals
368 //
369 namespace
370 { // required so we don't violate ODR
371 constexpr vtkIdType edges[4][2] = {
372 { 0, 1 },
373 { 1, 2 },
374 { 3, 2 },
375 { 0, 3 },
376 };
377
378 typedef int EDGE_LIST;
379 struct LINE_CASES_t
380 {
381 EDGE_LIST edges[5];
382 };
383 using LINE_CASES = struct LINE_CASES_t;
384
385 LINE_CASES lineCases[] = {
386 { { -1, -1, -1, -1, -1 } },
387 { { 0, 3, -1, -1, -1 } },
388 { { 1, 0, -1, -1, -1 } },
389 { { 1, 3, -1, -1, -1 } },
390 { { 2, 1, -1, -1, -1 } },
391 { { 0, 3, 2, 1, -1 } },
392 { { 2, 0, -1, -1, -1 } },
393 { { 2, 3, -1, -1, -1 } },
394 { { 3, 2, -1, -1, -1 } },
395 { { 0, 2, -1, -1, -1 } },
396 { { 1, 0, 3, 2, -1 } },
397 { { 1, 2, -1, -1, -1 } },
398 { { 3, 1, -1, -1, -1 } },
399 { { 0, 1, -1, -1, -1 } },
400 { { 3, 0, -1, -1, -1 } },
401 { { -1, -1, -1, -1, -1 } },
402 };
403 }
404
405 //------------------------------------------------------------------------------
GetEdgeArray(vtkIdType edgeId)406 const vtkIdType* vtkQuad::GetEdgeArray(vtkIdType edgeId)
407 {
408 return edges[edgeId];
409 }
410
411 //------------------------------------------------------------------------------
Contour(double value,vtkDataArray * cellScalars,vtkIncrementalPointLocator * locator,vtkCellArray * verts,vtkCellArray * lines,vtkCellArray * vtkNotUsed (polys),vtkPointData * inPd,vtkPointData * outPd,vtkCellData * inCd,vtkIdType cellId,vtkCellData * outCd)412 void vtkQuad::Contour(double value, vtkDataArray* cellScalars, vtkIncrementalPointLocator* locator,
413 vtkCellArray* verts, vtkCellArray* lines, vtkCellArray* vtkNotUsed(polys), vtkPointData* inPd,
414 vtkPointData* outPd, vtkCellData* inCd, vtkIdType cellId, vtkCellData* outCd)
415 {
416 static const int CASE_MASK[4] = { 1, 2, 4, 8 };
417 LINE_CASES* lineCase;
418 EDGE_LIST* edge;
419 int i, j, index;
420 const vtkIdType* vert;
421 int newCellId;
422 vtkIdType pts[2];
423 int e1, e2;
424 double t, x1[3], x2[3], x[3], deltaScalar;
425 vtkIdType offset = verts->GetNumberOfCells();
426
427 // Build the case table
428 for (i = 0, index = 0; i < 4; i++)
429 {
430 if (cellScalars->GetComponent(i, 0) >= value)
431 {
432 index |= CASE_MASK[i];
433 }
434 }
435
436 lineCase = lineCases + index;
437 edge = lineCase->edges;
438
439 for (; edge[0] > -1; edge += 2)
440 {
441 for (i = 0; i < 2; i++) // insert line
442 {
443 vert = edges[edge[i]];
444 // calculate a preferred interpolation direction
445 deltaScalar = (cellScalars->GetComponent(vert[1], 0) - cellScalars->GetComponent(vert[0], 0));
446 if (deltaScalar > 0)
447 {
448 e1 = vert[0];
449 e2 = vert[1];
450 }
451 else
452 {
453 e1 = vert[1];
454 e2 = vert[0];
455 deltaScalar = -deltaScalar;
456 }
457
458 // linear interpolation
459 if (deltaScalar == 0.0)
460 {
461 t = 0.0;
462 }
463 else
464 {
465 t = (value - cellScalars->GetComponent(e1, 0)) / deltaScalar;
466 }
467
468 this->Points->GetPoint(e1, x1);
469 this->Points->GetPoint(e2, x2);
470
471 for (j = 0; j < 3; j++)
472 {
473 x[j] = x1[j] + t * (x2[j] - x1[j]);
474 }
475 if (locator->InsertUniquePoint(x, pts[i]))
476 {
477 if (outPd)
478 {
479 vtkIdType p1 = this->PointIds->GetId(e1);
480 vtkIdType p2 = this->PointIds->GetId(e2);
481 outPd->InterpolateEdge(inPd, pts[i], p1, p2, t);
482 }
483 }
484 }
485 // check for degenerate line
486 if (pts[0] != pts[1])
487 {
488 newCellId = offset + lines->InsertNextCell(2, pts);
489 if (outCd)
490 {
491 outCd->CopyData(inCd, cellId, newCellId);
492 }
493 }
494 }
495 }
496
497 //------------------------------------------------------------------------------
GetEdge(int edgeId)498 vtkCell* vtkQuad::GetEdge(int edgeId)
499 {
500 int edgeIdPlus1 = edgeId + 1;
501
502 if (edgeIdPlus1 > 3)
503 {
504 edgeIdPlus1 = 0;
505 }
506
507 // load point id's
508 this->Line->PointIds->SetId(0, this->PointIds->GetId(edgeId));
509 this->Line->PointIds->SetId(1, this->PointIds->GetId(edgeIdPlus1));
510
511 // load coordinates
512 this->Line->Points->SetPoint(0, this->Points->GetPoint(edgeId));
513 this->Line->Points->SetPoint(1, this->Points->GetPoint(edgeIdPlus1));
514
515 return this->Line;
516 }
517
518 //------------------------------------------------------------------------------
519 // Intersect plane; see whether point is in quadrilateral. This code
520 // splits the quad into two triangles and intersects them (because the
521 // quad may be non-planar).
522 //
IntersectWithLine(const double p1[3],const double p2[3],double tol,double & t,double x[3],double pcoords[3],int & subId)523 int vtkQuad::IntersectWithLine(const double p1[3], const double p2[3], double tol, double& t,
524 double x[3], double pcoords[3], int& subId)
525 {
526 int diagonalCase;
527 double d1 = vtkMath::Distance2BetweenPoints(this->Points->GetPoint(0), this->Points->GetPoint(2));
528 double d2 = vtkMath::Distance2BetweenPoints(this->Points->GetPoint(1), this->Points->GetPoint(3));
529 subId = 0;
530
531 // Figure out how to uniquely tessellate the quad. Watch out for
532 // equivalent triangulations (i.e., the triangulation is equivalent
533 // no matter where the diagonal). In this case use the point ids as
534 // a tie breaker to ensure unique triangulation across the quad.
535 //
536 if (d1 == d2) // rare case; discriminate based on point id
537 {
538 int i, id, maxId = 0, maxIdx = 0;
539 for (i = 0; i < 4; i++) // find the maximum id
540 {
541 if ((id = this->PointIds->GetId(i)) > maxId)
542 {
543 maxId = id;
544 maxIdx = i;
545 }
546 }
547 if (maxIdx == 0 || maxIdx == 2)
548 diagonalCase = 0;
549 else
550 diagonalCase = 1;
551 }
552 else if (d1 < d2)
553 {
554 diagonalCase = 0;
555 }
556 else // d2 < d1
557 {
558 diagonalCase = 1;
559 }
560
561 // Note: in the following code the parametric coords must be adjusted to
562 // reflect the use of the triangle parametric coordinate system.
563 switch (diagonalCase)
564 {
565 case 0:
566 this->Triangle->Points->SetPoint(0, this->Points->GetPoint(0));
567 this->Triangle->Points->SetPoint(1, this->Points->GetPoint(1));
568 this->Triangle->Points->SetPoint(2, this->Points->GetPoint(2));
569 if (this->Triangle->IntersectWithLine(p1, p2, tol, t, x, pcoords, subId))
570 {
571 pcoords[0] = pcoords[0] + pcoords[1];
572 return 1;
573 }
574 this->Triangle->Points->SetPoint(0, this->Points->GetPoint(2));
575 this->Triangle->Points->SetPoint(1, this->Points->GetPoint(3));
576 this->Triangle->Points->SetPoint(2, this->Points->GetPoint(0));
577 if (this->Triangle->IntersectWithLine(p1, p2, tol, t, x, pcoords, subId))
578 {
579 pcoords[0] = 1.0 - (pcoords[0] + pcoords[1]);
580 pcoords[1] = 1.0 - pcoords[1];
581 return 1;
582 }
583 return 0;
584
585 case 1:
586 this->Triangle->Points->SetPoint(0, this->Points->GetPoint(0));
587 this->Triangle->Points->SetPoint(1, this->Points->GetPoint(1));
588 this->Triangle->Points->SetPoint(2, this->Points->GetPoint(3));
589 if (this->Triangle->IntersectWithLine(p1, p2, tol, t, x, pcoords, subId))
590 {
591 return 1;
592 }
593 this->Triangle->Points->SetPoint(0, this->Points->GetPoint(2));
594 this->Triangle->Points->SetPoint(1, this->Points->GetPoint(3));
595 this->Triangle->Points->SetPoint(2, this->Points->GetPoint(1));
596 if (this->Triangle->IntersectWithLine(p1, p2, tol, t, x, pcoords, subId))
597 {
598 pcoords[0] = 1.0 - pcoords[0];
599 pcoords[1] = 1.0 - pcoords[1];
600 return 1;
601 }
602
603 return 0;
604 }
605
606 return 0;
607 }
608
609 //------------------------------------------------------------------------------
Triangulate(int vtkNotUsed (index),vtkIdList * ptIds,vtkPoints * pts)610 int vtkQuad::Triangulate(int vtkNotUsed(index), vtkIdList* ptIds, vtkPoints* pts)
611 {
612 double d1, d2;
613
614 pts->Reset();
615 ptIds->Reset();
616
617 // use minimum diagonal (Delaunay triangles) - assumed convex
618 d1 = vtkMath::Distance2BetweenPoints(this->Points->GetPoint(0), this->Points->GetPoint(2));
619 d2 = vtkMath::Distance2BetweenPoints(this->Points->GetPoint(1), this->Points->GetPoint(3));
620
621 if (d1 <= d2)
622 {
623 ptIds->InsertId(0, this->PointIds->GetId(0));
624 pts->InsertPoint(0, this->Points->GetPoint(0));
625 ptIds->InsertId(1, this->PointIds->GetId(1));
626 pts->InsertPoint(1, this->Points->GetPoint(1));
627 ptIds->InsertId(2, this->PointIds->GetId(2));
628 pts->InsertPoint(2, this->Points->GetPoint(2));
629
630 ptIds->InsertId(3, this->PointIds->GetId(0));
631 pts->InsertPoint(3, this->Points->GetPoint(0));
632 ptIds->InsertId(4, this->PointIds->GetId(2));
633 pts->InsertPoint(4, this->Points->GetPoint(2));
634 ptIds->InsertId(5, this->PointIds->GetId(3));
635 pts->InsertPoint(5, this->Points->GetPoint(3));
636 }
637 else
638 {
639 ptIds->InsertId(0, this->PointIds->GetId(0));
640 pts->InsertPoint(0, this->Points->GetPoint(0));
641 ptIds->InsertId(1, this->PointIds->GetId(1));
642 pts->InsertPoint(1, this->Points->GetPoint(1));
643 ptIds->InsertId(2, this->PointIds->GetId(3));
644 pts->InsertPoint(2, this->Points->GetPoint(3));
645
646 ptIds->InsertId(3, this->PointIds->GetId(1));
647 pts->InsertPoint(3, this->Points->GetPoint(1));
648 ptIds->InsertId(4, this->PointIds->GetId(2));
649 pts->InsertPoint(4, this->Points->GetPoint(2));
650 ptIds->InsertId(5, this->PointIds->GetId(3));
651 pts->InsertPoint(5, this->Points->GetPoint(3));
652 }
653
654 return 1;
655 }
656
657 //------------------------------------------------------------------------------
Derivatives(int vtkNotUsed (subId),const double pcoords[3],const double * values,int dim,double * derivs)658 void vtkQuad::Derivatives(
659 int vtkNotUsed(subId), const double pcoords[3], const double* values, int dim, double* derivs)
660 {
661 double v0[2], v1[2], v2[2], v3[2], v10[3], v20[3], lenX;
662 double x0[3], x1[3], x2[3], x3[3], n[3], vec20[3], vec30[3];
663 double *J[2], J0[2], J1[2];
664 double *JI[2], JI0[2], JI1[2];
665 double funcDerivs[8], sum[2], dBydx, dBydy;
666 int i, j;
667
668 // Project points of quad into 2D system
669 this->Points->GetPoint(0, x0);
670 this->Points->GetPoint(1, x1);
671 this->Points->GetPoint(2, x2);
672 ComputeNormal(this, x0, x1, x2, n);
673 this->Points->GetPoint(3, x3);
674
675 for (i = 0; i < 3; i++)
676 {
677 v10[i] = x1[i] - x0[i];
678 vec20[i] = x2[i] - x0[i];
679 vec30[i] = x3[i] - x0[i];
680 }
681
682 vtkMath::Cross(n, v10, v20); // creates local y' axis
683
684 if ((lenX = vtkMath::Normalize(v10)) <= 0.0 || vtkMath::Normalize(v20) <= 0.0) // degenerate
685 {
686 for (j = 0; j < dim; j++)
687 {
688 for (i = 0; i < 3; i++)
689 {
690 derivs[j * dim + i] = 0.0;
691 }
692 }
693 return;
694 }
695
696 v0[0] = v0[1] = 0.0; // convert points to 2D (i.e., local system)
697 v1[0] = lenX;
698 v1[1] = 0.0;
699 v2[0] = vtkMath::Dot(vec20, v10);
700 v2[1] = vtkMath::Dot(vec20, v20);
701 v3[0] = vtkMath::Dot(vec30, v10);
702 v3[1] = vtkMath::Dot(vec30, v20);
703
704 this->InterpolationDerivs(pcoords, funcDerivs);
705
706 // Compute Jacobian and inverse Jacobian
707 J[0] = J0;
708 J[1] = J1;
709 JI[0] = JI0;
710 JI[1] = JI1;
711
712 J[0][0] =
713 v0[0] * funcDerivs[0] + v1[0] * funcDerivs[1] + v2[0] * funcDerivs[2] + v3[0] * funcDerivs[3];
714 J[0][1] =
715 v0[1] * funcDerivs[0] + v1[1] * funcDerivs[1] + v2[1] * funcDerivs[2] + v3[1] * funcDerivs[3];
716 J[1][0] =
717 v0[0] * funcDerivs[4] + v1[0] * funcDerivs[5] + v2[0] * funcDerivs[6] + v3[0] * funcDerivs[7];
718 J[1][1] =
719 v0[1] * funcDerivs[4] + v1[1] * funcDerivs[5] + v2[1] * funcDerivs[6] + v3[1] * funcDerivs[7];
720
721 // Compute inverse Jacobian, return if Jacobian is singular
722 if (!vtkMath::InvertMatrix(J, JI, 2))
723 {
724 for (j = 0; j < dim; j++)
725 {
726 for (i = 0; i < 3; i++)
727 {
728 derivs[j * dim + i] = 0.0;
729 }
730 }
731 return;
732 }
733
734 // Loop over "dim" derivative values. For each set of values,
735 // compute derivatives
736 // in local system and then transform into modelling system.
737 // First compute derivatives in local x'-y' coordinate system
738 for (j = 0; j < dim; j++)
739 {
740 sum[0] = sum[1] = 0.0;
741 for (i = 0; i < 4; i++) // loop over interp. function derivatives
742 {
743 sum[0] += funcDerivs[i] * values[dim * i + j];
744 sum[1] += funcDerivs[4 + i] * values[dim * i + j];
745 }
746 dBydx = sum[0] * JI[0][0] + sum[1] * JI[0][1];
747 dBydy = sum[0] * JI[1][0] + sum[1] * JI[1][1];
748
749 // Transform into global system (dot product with global axes)
750 derivs[3 * j] = dBydx * v10[0] + dBydy * v20[0];
751 derivs[3 * j + 1] = dBydx * v10[1] + dBydy * v20[1];
752 derivs[3 * j + 2] = dBydx * v10[2] + dBydy * v20[2];
753 }
754 }
755
756 //------------------------------------------------------------------------------
757 // support quad clipping
758 typedef int QUAD_EDGE_LIST;
759 struct QUAD_CASES_t
760 {
761 QUAD_EDGE_LIST edges[14];
762 };
763 using QUAD_CASES = struct QUAD_CASES_t;
764
765 static QUAD_CASES quadCases[] = {
766 { { -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 } }, // 0
767 { { 3, 100, 0, 3, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 } }, // 1
768 { { 3, 101, 1, 0, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 } }, // 2
769 { { 4, 100, 101, 1, 3, -1, -1, -1, -1, -1, -1, -1, -1, -1 } }, // 3
770 { { 3, 102, 2, 1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 } }, // 4
771 { { 3, 100, 0, 3, 3, 102, 2, 1, 4, 0, 1, 2, 3, -1 } }, // 5
772 { { 4, 101, 102, 2, 0, -1, -1, -1, -1, -1, -1, -1, -1, -1 } }, // 6
773 { { 3, 100, 101, 3, 3, 101, 2, 3, 3, 101, 102, 2, -1, -1 } }, // 7
774 { { 3, 103, 3, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 } }, // 8
775 { { 4, 100, 0, 2, 103, -1, -1, -1, -1, -1, -1, -1, -1, -1 } }, // 9
776 { { 3, 101, 1, 0, 3, 103, 3, 2, 4, 0, 1, 2, 3, -1 } }, // 10
777 { { 3, 100, 101, 1, 3, 100, 1, 2, 3, 100, 2, 103, -1, -1 } }, // 11
778 { { 4, 102, 103, 3, 1, -1, -1, -1, -1, -1, -1, -1, -1, -1 } }, // 12
779 { { 3, 100, 0, 103, 3, 0, 1, 103, 3, 1, 102, 103, -1, -1 } }, // 13
780 { { 3, 0, 101, 102, 3, 0, 102, 3, 3, 102, 103, 3, -1, -1 } }, // 14
781 { { 4, 100, 101, 102, 103, -1, -1, -1, -1, -1, -1, -1, -1, -1 } }, // 15
782 };
783
784 static QUAD_CASES quadCasesComplement[] = {
785 { { -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 } }, // 0
786 { { 3, 100, 0, 3, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 } }, // 1
787 { { 3, 101, 1, 0, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 } }, // 2
788 { { 4, 100, 101, 1, 3, -1, -1, -1, -1, -1, -1, -1, -1, -1 } }, // 3
789 { { 3, 102, 2, 1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 } }, // 4
790 { { 3, 100, 0, 3, 3, 102, 2, 1, -1, -1, -1, -1, -1, -1 } }, // 5
791 { { 4, 101, 102, 2, 0, -1, -1, -1, -1, -1, -1, -1, -1, -1 } }, // 6
792 { { 3, 100, 101, 3, 3, 101, 2, 3, 3, 101, 102, 2, -1, -1 } }, // 7
793 { { 3, 103, 3, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 } }, // 8
794 { { 4, 100, 0, 2, 103, -1, -1, -1, -1, -1, -1, -1, -1, -1 } }, // 9
795 { { 3, 101, 1, 0, 3, 103, 3, 2, -1, -1, -1, -1, -1, -1 } }, // 10
796 { { 3, 100, 101, 1, 3, 100, 1, 2, 3, 100, 2, 103, -1, -1 } }, // 11
797 { { 4, 102, 103, 3, 1, -1, -1, -1, -1, -1, -1, -1, -1, -1 } }, // 12
798 { { 3, 100, 0, 103, 3, 0, 1, 103, 3, 1, 102, 103, -1, -1 } }, // 13
799 { { 3, 0, 101, 102, 3, 0, 102, 3, 3, 102, 103, 3, -1, -1 } }, // 14
800 { { 4, 100, 101, 102, 103, -1, -1, -1, -1, -1, -1, -1, -1, -1 } }, // 15
801 };
802
803 //------------------------------------------------------------------------------
804 // Clip this quad using scalar value provided. Like contouring, except
805 // that it cuts the quad to produce other quads and/or triangles.
Clip(double value,vtkDataArray * cellScalars,vtkIncrementalPointLocator * locator,vtkCellArray * polys,vtkPointData * inPd,vtkPointData * outPd,vtkCellData * inCd,vtkIdType cellId,vtkCellData * outCd,int insideOut)806 void vtkQuad::Clip(double value, vtkDataArray* cellScalars, vtkIncrementalPointLocator* locator,
807 vtkCellArray* polys, vtkPointData* inPd, vtkPointData* outPd, vtkCellData* inCd, vtkIdType cellId,
808 vtkCellData* outCd, int insideOut)
809 {
810 static const int CASE_MASK[4] = { 1, 2, 4, 8 };
811 QUAD_CASES* quadCase;
812 QUAD_EDGE_LIST* edge;
813 int i, j, index;
814 const vtkIdType* vert;
815 int e1, e2;
816 int newCellId;
817 vtkIdType pts[4];
818 int vertexId;
819 double t, x1[3], x2[3], x[3], deltaScalar;
820 double scalar0, scalar1, e1Scalar;
821
822 // Build the index into the case table
823 if (insideOut)
824 {
825 for (i = 0, index = 0; i < 4; i++)
826 {
827 if (cellScalars->GetComponent(i, 0) <= value)
828 {
829 index |= CASE_MASK[i];
830 }
831 }
832 // Select case based on the index and get the list of edges for this case
833 quadCase = quadCases + index;
834 }
835 else
836 {
837 for (i = 0, index = 0; i < 4; i++)
838 {
839 if (cellScalars->GetComponent(i, 0) > value)
840 {
841 index |= CASE_MASK[i];
842 }
843 }
844 // Select case based on the index and get the list of edges for this case
845 quadCase = quadCasesComplement + index;
846 }
847
848 edge = quadCase->edges;
849
850 // generate each quad
851 for (; edge[0] > -1; edge += edge[0] + 1)
852 {
853 for (i = 0; i < edge[0]; i++) // insert quad or triangle
854 {
855 // vertex exists, and need not be interpolated
856 if (edge[i + 1] >= 100)
857 {
858 vertexId = edge[i + 1] - 100;
859 this->Points->GetPoint(vertexId, x);
860 if (locator->InsertUniquePoint(x, pts[i]))
861 {
862 outPd->CopyData(inPd, this->PointIds->GetId(vertexId), pts[i]);
863 }
864 }
865
866 else // new vertex, interpolate
867 {
868 vert = edges[edge[i + 1]];
869
870 // calculate a preferred interpolation direction
871 scalar0 = cellScalars->GetComponent(vert[0], 0);
872 scalar1 = cellScalars->GetComponent(vert[1], 0);
873 deltaScalar = scalar1 - scalar0;
874
875 if (deltaScalar > 0)
876 {
877 e1 = vert[0];
878 e2 = vert[1];
879 e1Scalar = scalar0;
880 }
881 else
882 {
883 e1 = vert[1];
884 e2 = vert[0];
885 e1Scalar = scalar1;
886 deltaScalar = -deltaScalar;
887 }
888
889 // linear interpolation
890 if (deltaScalar == 0.0)
891 {
892 t = 0.0;
893 }
894 else
895 {
896 t = (value - e1Scalar) / deltaScalar;
897 }
898
899 this->Points->GetPoint(e1, x1);
900 this->Points->GetPoint(e2, x2);
901
902 for (j = 0; j < 3; j++)
903 {
904 x[j] = x1[j] + t * (x2[j] - x1[j]);
905 }
906
907 if (locator->InsertUniquePoint(x, pts[i]))
908 {
909 vtkIdType p1 = this->PointIds->GetId(e1);
910 vtkIdType p2 = this->PointIds->GetId(e2);
911 outPd->InterpolateEdge(inPd, pts[i], p1, p2, t);
912 }
913 }
914 }
915 // check for degenerate output
916 if (edge[0] == 3) // i.e., a triangle
917 {
918 if (pts[0] == pts[1] || pts[0] == pts[2] || pts[1] == pts[2])
919 {
920 continue;
921 }
922 }
923 else // a quad
924 {
925 if ((pts[0] == pts[3] && pts[1] == pts[2]) || (pts[0] == pts[1] && pts[3] == pts[2]))
926 {
927 continue;
928 }
929 }
930
931 newCellId = polys->InsertNextCell(edge[0], pts);
932 outCd->CopyData(inCd, cellId, newCellId);
933 }
934 }
935
936 //------------------------------------------------------------------------------
937 static double vtkQuadCellPCoords[12] = {
938 0.0, 0.0, 0.0, //
939 1.0, 0.0, 0.0, //
940 1.0, 1.0, 0.0, //
941 0.0, 1.0, 0.0 //
942 };
GetParametricCoords()943 double* vtkQuad::GetParametricCoords()
944 {
945 return vtkQuadCellPCoords;
946 }
947
948 //------------------------------------------------------------------------------
PrintSelf(ostream & os,vtkIndent indent)949 void vtkQuad::PrintSelf(ostream& os, vtkIndent indent)
950 {
951 this->Superclass::PrintSelf(os, indent);
952
953 os << indent << "Line:\n";
954 this->Line->PrintSelf(os, indent.GetNextIndent());
955 os << indent << "Triangle:\n";
956 this->Triangle->PrintSelf(os, indent.GetNextIndent());
957 }
958