1 /*=========================================================================
2
3 Program: Visualization Toolkit
4 Module: vtkLine.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 "vtkLine.h"
16
17 #include "vtkCellArray.h"
18 #include "vtkCellData.h"
19 #include "vtkMath.h"
20 #include "vtkObjectFactory.h"
21 #include "vtkPointData.h"
22 #include "vtkIncrementalPointLocator.h"
23 #include "vtkPoints.h"
24
25 vtkStandardNewMacro(vtkLine);
26
27 //----------------------------------------------------------------------------
28 // Construct the line with two points.
vtkLine()29 vtkLine::vtkLine()
30 {
31 this->Points->SetNumberOfPoints(2);
32 this->PointIds->SetNumberOfIds(2);
33 for (int i = 0; i < 2; i++)
34 {
35 this->Points->SetPoint(i, 0.0, 0.0, 0.0);
36 this->PointIds->SetId(i,0);
37 }
38 }
39
40 //----------------------------------------------------------------------------
41 static const int VTK_NO_INTERSECTION=0;
42 static const int VTK_YES_INTERSECTION=2;
43 static const int VTK_ON_LINE=3;
44
EvaluatePosition(const double x[3],double closestPoint[3],int & subId,double pcoords[3],double & dist2,double weights[])45 int vtkLine::EvaluatePosition(const double x[3], double closestPoint[3],
46 int& subId, double pcoords[3],
47 double& dist2, double weights[])
48 {
49 double a1[3], a2[3];
50
51 subId = 0;
52 pcoords[0] = pcoords[1] = pcoords[2] = 0.0;
53
54 this->Points->GetPoint(0, a1);
55 this->Points->GetPoint(1, a2);
56
57 // DistanceToLine sets pcoords[0] to a value t
58 dist2 = this->DistanceToLine(x,a1,a2,pcoords[0],closestPoint);
59
60 // pcoords[0] == t, need weights to be 1-t and t
61 weights[0] = 1.0 - pcoords[0];
62 weights[1] = pcoords[0];
63
64 if ( pcoords[0] < 0.0 || pcoords[0] > 1.0 )
65 {
66 return 0;
67 }
68 else
69 {
70 return 1;
71 }
72 }
73
74 //----------------------------------------------------------------------------
EvaluateLocation(int & vtkNotUsed (subId),const double pcoords[3],double x[3],double * weights)75 void vtkLine::EvaluateLocation(int& vtkNotUsed(subId), const double pcoords[3],
76 double x[3], double *weights)
77 {
78 int i;
79 double a1[3], a2[3];
80 this->Points->GetPoint(0, a1);
81 this->Points->GetPoint(1, a2);
82
83 for (i=0; i<3; i++)
84 {
85 x[i] = a1[i] + pcoords[0]*(a2[i] - a1[i]);
86 }
87
88 weights[0] = 1.0 - pcoords[0];
89 weights[1] = pcoords[0];
90 }
91
92 //----------------------------------------------------------------------------
93 // Performs intersection of the projection of two finite 3D lines onto a 2D
94 // plane. An intersection is found if the projection of the two lines onto
95 // the plane perpendicular to the cross product of the two lines intersect.
96 // The parameters (u,v) are the parametric coordinates of the lines at the
97 // position of closest approach.
Intersection(const double a1[3],const double a2[3],const double b1[3],const double b2[3],double & u,double & v)98 int vtkLine::Intersection (const double a1[3], const double a2[3],
99 const double b1[3], const double b2[3],
100 double& u, double& v)
101 {
102 double a21[3], b21[3], b1a1[3];
103 double c[2];
104 double *A[2], row1[2], row2[2];
105
106 // Initialize
107 u = v = 0.0;
108
109 // Determine line vectors.
110 a21[0] = a2[0] - a1[0]; b21[0] = b2[0] - b1[0]; b1a1[0] = b1[0] - a1[0];
111 a21[1] = a2[1] - a1[1]; b21[1] = b2[1] - b1[1]; b1a1[1] = b1[1] - a1[1];
112 a21[2] = a2[2] - a1[2]; b21[2] = b2[2] - b1[2]; b1a1[2] = b1[2] - a1[2];
113
114 // Compute the system (least squares) matrix.
115 A[0] = row1;
116 A[1] = row2;
117 row1[0] = vtkMath::Dot ( a21, a21 );
118 row1[1] = -vtkMath::Dot ( a21, b21 );
119 row2[0] = row1[1];
120 row2[1] = vtkMath::Dot ( b21, b21 );
121
122 // Compute the least squares system constant term.
123 c[0] = vtkMath::Dot ( a21, b1a1 );
124 c[1] = -vtkMath::Dot ( b21, b1a1 );
125
126
127 // Solve the system of equations
128 if ( vtkMath::SolveLinearSystem(A,c,2) == 0 )
129 {
130 // The lines are colinear. Therefore, one of the four endpoints is the
131 // point of closest approach
132 double minDist = VTK_DOUBLE_MAX;
133 const double* p[4] = {a1,a2,b1,b2};
134 const double* l1[4] = {b1,b1,a1,a1};
135 const double* l2[4] = {b2,b2,a2,a2};
136 double* uv1[4] = {&v,&v,&u,&u};
137 double* uv2[4] = {&u,&u,&v,&v};
138 double t=0;
139 for (unsigned i=0;i<4;i++)
140 {
141 double dist = vtkLine::DistanceToLine(p[i],l1[i],l2[i],t);
142 if (dist < minDist)
143 {
144 minDist = dist;
145 *(uv1[i]) = t;
146 *(uv2[i]) = static_cast<double>(i%2); // the corresponding extremum
147 }
148 }
149 return VTK_ON_LINE;
150 }
151 else
152 {
153 u = c[0];
154 v = c[1];
155 }
156
157 // Check parametric coordinates for intersection.
158 if ( (0.0 <= u) && (u <= 1.0) && (0.0 <= v) && (v <= 1.0) )
159 {
160 return VTK_YES_INTERSECTION;
161 }
162 else
163 {
164 return VTK_NO_INTERSECTION;
165 }
166 }
167
168 //----------------------------------------------------------------------------
Intersection3D(double a1[3],double a2[3],double b1[3],double b2[3],double & u,double & v)169 int vtkLine::Intersection3D(double a1[3], double a2[3],
170 double b1[3], double b2[3],
171 double& u, double& v)
172 {
173 // Description:
174 // Performs intersection of two finite 3D lines. An intersection is found if
175 // the projection of the two lines onto the plane perpendicular to the cross
176 // product of the two lines intersect, and if the distance between the
177 // closest points of approach are within a relative tolerance. The parameters
178 // (u,v) are the parametric coordinates of the lines at the position of
179 // closest approach.
180 int projectedIntersection = vtkLine::Intersection(a1,a2,b1,b2,u,v);
181
182 if (projectedIntersection == VTK_YES_INTERSECTION)
183 {
184 double a_i,b_i;
185 double lenA = 0.;
186 double lenB = 0.;
187 double dist = 0.;
188 for (unsigned i=0;i<3;i++)
189 {
190 a_i = a1[i] + (a2[i]-a1[i])*u;
191 b_i = b1[i] + (b2[i]-b1[i])*v;
192 lenA += (a2[i]-a1[i])*(a2[i]-a1[i]);
193 lenB += (b2[i]-b1[i])*(b2[i]-b1[i]);
194 dist += (a_i-b_i)*(a_i-b_i);
195 }
196 if (dist > 1.e-6*(lenA > lenB ? lenA : lenB))
197 return VTK_NO_INTERSECTION;
198 }
199
200 return projectedIntersection;
201 }
202
203 //----------------------------------------------------------------------------
CellBoundary(int vtkNotUsed (subId),const double pcoords[3],vtkIdList * pts)204 int vtkLine::CellBoundary(int vtkNotUsed(subId), const double pcoords[3],
205 vtkIdList *pts)
206 {
207 pts->SetNumberOfIds(1);
208
209 if ( pcoords[0] >= 0.5 )
210 {
211 pts->SetId(0,this->PointIds->GetId(1));
212 if ( pcoords[0] > 1.0 )
213 {
214 return 0;
215 }
216 else
217 {
218 return 1;
219 }
220 }
221 else
222 {
223 pts->SetId(0,this->PointIds->GetId(0));
224 if ( pcoords[0] < 0.0 )
225 {
226 return 0;
227 }
228 else
229 {
230 return 1;
231 }
232 }
233 }
234
235 //----------------------------------------------------------------------------
236 //
237 // marching lines case table
238 //
239 typedef int VERT_LIST;
240
241 typedef struct {
242 VERT_LIST verts[2];
243 } VERT_CASES;
244
245 static VERT_CASES vertCases[4]= {
246 {{-1,-1}},
247 {{1,0}},
248 {{0,1}},
249 {{-1,-1}}};
250
251 //----------------------------------------------------------------------------
Contour(double value,vtkDataArray * cellScalars,vtkIncrementalPointLocator * locator,vtkCellArray * verts,vtkCellArray * vtkNotUsed (lines),vtkCellArray * vtkNotUsed (polys),vtkPointData * inPd,vtkPointData * outPd,vtkCellData * inCd,vtkIdType cellId,vtkCellData * outCd)252 void vtkLine::Contour(double value, vtkDataArray *cellScalars,
253 vtkIncrementalPointLocator *locator, vtkCellArray *verts,
254 vtkCellArray *vtkNotUsed(lines),
255 vtkCellArray *vtkNotUsed(polys),
256 vtkPointData *inPd, vtkPointData *outPd,
257 vtkCellData *inCd, vtkIdType cellId, vtkCellData *outCd)
258 {
259 static const int CASE_MASK[2] = {1,2};
260 int index, i, newCellId;
261 VERT_CASES *vertCase;
262 VERT_LIST *vert;
263 double t, x[3], x1[3], x2[3];
264 vtkIdType pts[1];
265
266 //
267 // Build the case table
268 //
269 for ( i=0, index = 0; i < 2; i++)
270 {
271 if (cellScalars->GetComponent(i,0) >= value)
272 {
273 index |= CASE_MASK[i];
274 }
275 }
276
277 vertCase = vertCases + index;
278 vert = vertCase->verts;
279
280 if ( vert[0] > -1 )
281 {
282 t = (value - cellScalars->GetComponent(vert[0],0)) /
283 (cellScalars->GetComponent(vert[1],0) -
284 cellScalars->GetComponent(vert[0],0));
285 this->Points->GetPoint(vert[0], x1);
286 this->Points->GetPoint(vert[1], x2);
287 for (i=0; i<3; i++)
288 {
289 x[i] = x1[i] + t * (x2[i] - x1[i]);
290 }
291
292 if ( locator->InsertUniquePoint(x, pts[0]) )
293 {
294 if ( outPd )
295 {
296 vtkIdType p1 = this->PointIds->GetId(vert[0]);
297 vtkIdType p2 = this->PointIds->GetId(vert[1]);
298 outPd->InterpolateEdge(inPd,pts[0],p1,p2,t);
299 }
300 }
301 newCellId = verts->InsertNextCell(1,pts);
302 if (outCd)
303 {
304 outCd->CopyData(inCd, cellId, newCellId);
305 }
306 }
307 }
308
309 //----------------------------------------------------------------------------
DistanceBetweenLines(double l0[3],double l1[3],double m0[3],double m1[3],double closestPt1[3],double closestPt2[3],double & t1,double & t2)310 double vtkLine::DistanceBetweenLines(
311 double l0[3], double l1[3], // line 1
312 double m0[3], double m1[3], // line 2
313 double closestPt1[3], double closestPt2[3], // closest points
314 double &t1, double &t2 ) // parametric coords of the closest points
315 {
316 // Part of this function was adapted from "GeometryAlgorithms.com"
317 //
318 // Copyright 2001, softSurfer (www.softsurfer.com)
319 // This code may be freely used and modified for any purpose
320 // providing that this copyright notice is included with it.
321 // SoftSurfer makes no warranty for this code, and cannot be held
322 // liable for any real or imagined damage resulting from its use.
323 // Users of this code must verify correctness for their application.
324
325 const double u[3] = { l1[0]-l0[0], l1[1]-l0[1], l1[2]-l0[2] };
326 const double v[3] = { m1[0]-m0[0], m1[1]-m0[1], m1[2]-m0[2] };
327 const double w[3] = { l0[0]-m0[0], l0[1]-m0[1], l0[2]-m0[2] };
328 const double a = vtkMath::Dot(u,u);
329 const double b = vtkMath::Dot(u,v);
330 const double c = vtkMath::Dot(v,v); // always >= 0
331 const double d = vtkMath::Dot(u,w);
332 const double e = vtkMath::Dot(v,w);
333 const double D = a*c - b*b; // always >= 0
334
335 // compute the line parameters of the two closest points
336 if (D < 1e-6)
337 { // the lines are almost parallel
338 t1 = 0.0;
339 t2 = (b > c ? d/b : e/c); // use the largest denominator
340 }
341 else
342 {
343 t1 = (b*e - c*d) / D;
344 t2 = (a*e - b*d) / D;
345 }
346
347 for (unsigned int i = 0; i < 3; i++)
348 {
349 closestPt1[i] = l0[i] + t1 * u[i];
350 closestPt2[i] = m0[i] + t2 * v[i];
351 }
352
353 // Return the distance squared between the lines =
354 // the mag squared of the distance between the two closest points
355 // = L1(t1) - L2(t2)
356 return vtkMath::Distance2BetweenPoints( closestPt1, closestPt2 );
357 }
358
359 //----------------------------------------------------------------------------
DistanceBetweenLineSegments(double l0[3],double l1[3],double m0[3],double m1[3],double closestPt1[3],double closestPt2[3],double & t1,double & t2)360 double vtkLine::DistanceBetweenLineSegments(
361 double l0[3], double l1[3], // line segment 1
362 double m0[3], double m1[3], // line segment 2
363 double closestPt1[3], double closestPt2[3], // closest points
364 double &t1, double &t2 ) // parametric coords
365 // of the closest points
366 {
367 // Part of this function was adapted from "GeometryAlgorithms.com"
368 //
369 // Copyright 2001, softSurfer (www.softsurfer.com)
370 // This code may be freely used and modified for any purpose
371 // providing that this copyright notice is included with it.
372 // SoftSurfer makes no warranty for this code, and cannot be held
373 // liable for any real or imagined damage resulting from its use.
374 // Users of this code must verify correctness for their application.
375
376 const double u[3] = { l1[0]-l0[0], l1[1]-l0[1], l1[2]-l0[2] };
377 const double v[3] = { m1[0]-m0[0], m1[1]-m0[1], m1[2]-m0[2] };
378 const double w[3] = { l0[0]-m0[0], l0[1]-m0[1], l0[2]-m0[2] };
379 const double a = vtkMath::Dot(u,u);
380 const double b = vtkMath::Dot(u,v);
381 const double c = vtkMath::Dot(v,v); // always >= 0
382 const double d = vtkMath::Dot(u,w);
383 const double e = vtkMath::Dot(v,w);
384 const double D = a*c - b*b; // always >= 0
385 double sN, sD = D; // sc = sN / sD, default sD = D >= 0
386 double tN, tD = D; // tc = tN / tD, default tD = D >= 0
387
388 // compute the line parameters of the two closest points
389
390 if (D < 1e-6)
391 {
392 // The lines are colinear. Therefore, one of the four endpoints is the
393 // point of closest approach
394 double minDist = VTK_DOUBLE_MAX;
395 double* p[4] = {l0,l1,m0,m1};
396 double* a1[4] = {m0,m0,l0,l0};
397 double* a2[4] = {m1,m1,l1,l1};
398 double* uv1[4] = {&t2,&t2,&t1,&t1};
399 double* uv2[4] = {&t1,&t1,&t2,&t2};
400 double* pn1[4] = {closestPt2,closestPt2,closestPt1,closestPt1};
401 double* pn2[4] = {closestPt1,closestPt1,closestPt2,closestPt2};
402 double dist,pn_[3];
403 for (unsigned i=0;i<4;i++)
404 {
405 double t = 0;
406 dist = vtkLine::DistanceToLine(p[i],a1[i],a2[i],t,pn_);
407 if (dist < minDist)
408 {
409 minDist = dist;
410 *(uv1[i]) = (t < 0. ? 0. : (t > 1. ? 1. : t));
411 *(uv2[i]) = static_cast<double>(i%2); // the corresponding extremum
412 for (unsigned j=0;j<3;j++)
413 {
414 pn1[i][j] = pn_[j];
415 pn2[i][j] = p[i][j];
416 }
417 }
418 }
419 return minDist;
420 }
421
422 // The lines aren't parallel.
423
424 else
425 { // get the closest points on the infinite lines
426 sN = (b*e - c*d);
427 tN = (a*e - b*d);
428 if (sN < 0.0)
429 { // sc < 0 => the s=0 edge is visible
430 sN = 0.0;
431 tN = e;
432 tD = c;
433 }
434 else if (sN > sD)
435 { // sc > 1 => the s=1 edge is visible
436 sN = sD;
437 tN = e + b;
438 tD = c;
439 }
440 }
441
442 if (tN < 0.0)
443 { // tc < 0 => the t=0 edge is visible
444 tN = 0.0;
445
446 // recompute sc for this edge
447 if (-d < 0.0)
448 {
449 sN = 0.0;
450 }
451 else if (-d > a)
452 {
453 sN = sD;
454 }
455 else
456 {
457 sN = -d;
458 sD = a;
459 }
460 }
461 else if (tN > tD)
462 { // tc > 1 => the t=1 edge is visible
463 tN = tD;
464
465 // recompute sc for this edge
466 if ((-d + b) < 0.0)
467 {
468 sN = 0;
469 }
470 else if ((-d + b) > a)
471 {
472 sN = sD;
473 }
474 else
475 {
476 sN = (-d + b);
477 sD = a;
478 }
479 }
480
481 // finally do the division to get sc and tc
482 t1 = (fabs(sN) < 1e-6 ? 0.0 : sN / sD);
483 t2 = (fabs(tN) < 1e-6 ? 0.0 : tN / tD);
484
485 // Closest Point on segment1 = S1(t1) = l0 + t1*u
486 // Closest Point on segment2 = S1(t2) = m0 + t2*v
487
488 for (unsigned int i = 0; i < 3; i++)
489 {
490 closestPt1[i] = l0[i] + t1 * u[i];
491 closestPt2[i] = m0[i] + t2 * v[i];
492 }
493
494 // Return the distance squared between the lines =
495 // the mag squared of the distance between the two closest points
496 // = S1(t1) - S2(t2)
497 return vtkMath::Distance2BetweenPoints( closestPt1, closestPt2 );
498 }
499
500 //----------------------------------------------------------------------------
501 // Compute distance to finite line. Returns parametric coordinate t
502 // and point location on line.
DistanceToLine(const double x[3],const double p1[3],const double p2[3],double & t,double closestPoint[3])503 double vtkLine::DistanceToLine(const double x[3], const double p1[3], const double p2[3],
504 double &t, double closestPoint[3])
505 {
506 const double *closest = nullptr;
507 //
508 // Determine appropriate vectors
509 //
510 double p21[3] = { p2[0]- p1[0], p2[1]- p1[1], p2[2]- p1[2] };
511
512 //
513 // Get parametric location
514 //
515 double num = p21[0]*(x[0]-p1[0]) + p21[1]*(x[1]-p1[1]) + p21[2]*(x[2]-p1[2]);
516 if (num == 0.0)
517 {
518 t = 0;
519 closest = p1;
520 }
521 else
522 {
523 double denom = vtkMath::Dot(p21,p21);
524
525 // trying to avoid an expensive fabs
526 double tolerance = VTK_TOL*num;
527 if (tolerance < 0.0)
528 {
529 tolerance = -tolerance;
530 }
531 if ( denom < tolerance ) //numerically bad!
532 {
533 if (num > 0)
534 {
535 closest = p2;
536 t = VTK_DOUBLE_MAX;
537 }
538 else
539 {
540 closest = p1;
541 t = VTK_DOUBLE_MIN;
542 }
543 }
544 //
545 // If parametric coordinate is within 0<=p<=1, then the point is closest to
546 // the line. Otherwise, it's closest to a point at the end of the line.
547 //
548 else if ( (t=num/denom) < 0.0 )
549 {
550 closest = p1;
551 }
552 else if ( t > 1.0 )
553 {
554 closest = p2;
555 }
556 else
557 {
558 closest = p21;
559 p21[0] = p1[0] + t*p21[0];
560 p21[1] = p1[1] + t*p21[1];
561 p21[2] = p1[2] + t*p21[2];
562 }
563 }
564
565 if (closestPoint)
566 {
567 closestPoint[0] = closest[0];
568 closestPoint[1] = closest[1];
569 closestPoint[2] = closest[2];
570 }
571 return vtkMath::Distance2BetweenPoints(closest,x);
572 }
573
574 //----------------------------------------------------------------------------
575 //
576 // Determine the distance of the current vertex to the edge defined by
577 // the vertices provided. Returns distance squared. Note: line is assumed
578 // infinite in extent.
579 //
DistanceToLine(const double x[3],const double p1[3],const double p2[3])580 double vtkLine::DistanceToLine (const double x[3], const double p1[3], const double p2[3])
581 {
582 int i;
583 double np1[3], p1p2[3], proj, den;
584
585 for (i=0; i<3; i++)
586 {
587 np1[i] = x[i] - p1[i];
588 p1p2[i] = p1[i] - p2[i];
589 }
590
591 if ( (den=vtkMath::Norm(p1p2)) != 0.0 )
592 {
593 for (i=0; i<3; i++)
594 {
595 p1p2[i] /= den;
596 }
597 }
598 else
599 {
600 return vtkMath::Dot(np1,np1);
601 }
602
603 proj = vtkMath::Dot(np1,p1p2);
604
605 return (vtkMath::Dot(np1,np1) - proj*proj);
606 }
607
608 //----------------------------------------------------------------------------
609 // Line-line intersection. Intersection has to occur within [0,1] parametric
610 // coordinates and with specified tolerance.
IntersectWithLine(const double p1[3],const double p2[3],double tol,double & t,double x[3],double pcoords[3],int & subId)611 int vtkLine::IntersectWithLine(const double p1[3], const double p2[3], double tol, double& t,
612 double x[3], double pcoords[3], int& subId)
613 {
614 double a1[3], a2[3];
615 double projXYZ[3];
616 int i;
617
618 subId = 0;
619 pcoords[1] = pcoords[2] = 0.0;
620
621 this->Points->GetPoint(0, a1);
622 this->Points->GetPoint(1, a2);
623
624 if ( this->Intersection(p1, p2, a1, a2, t, pcoords[0]) ==
625 VTK_YES_INTERSECTION )
626 {
627 // make sure we are within tolerance
628 for (i=0; i<3; i++)
629 {
630 x[i] = a1[i] + pcoords[0]*(a2[i]-a1[i]);
631 projXYZ[i] = p1[i] + t*(p2[i]-p1[i]);
632 }
633 if ( vtkMath::Distance2BetweenPoints(x,projXYZ) <= tol*tol )
634 {
635 return 1;
636 }
637 else
638 {
639 return 0;
640 }
641 }
642
643 else //check to see if it lies within tolerance
644 {
645 // one of the parametric coords must be outside 0-1
646 if (t < 0.0)
647 {
648 t = 0.0;
649 if (vtkLine::DistanceToLine(p1,a1,a2,pcoords[0],x) <= tol*tol)
650 {
651 return 1;
652 }
653 else
654 {
655 return 0;
656 }
657 }
658 if (t > 1.0)
659 {
660 t = 1.0;
661 if (vtkLine::DistanceToLine(p2,a1,a2,pcoords[0],x) <= tol*tol)
662 {
663 return 1;
664 }
665 else
666 {
667 return 0;
668 }
669 }
670 if (pcoords[0] < 0.0)
671 {
672 pcoords[0] = 0.0;
673 if (vtkLine::DistanceToLine(a1,p1,p2,t,x) <= tol*tol)
674 {
675 return 1;
676 }
677 else
678 {
679 return 0;
680 }
681 }
682 if (pcoords[0] > 1.0)
683 {
684 pcoords[0] = 1.0;
685 if (vtkLine::DistanceToLine(a2,p1,p2,t,x) <= tol*tol)
686 {
687 return 1;
688 }
689 else
690 {
691 return 0;
692 }
693 }
694 }
695 return 0;
696 }
697
698 //----------------------------------------------------------------------------
Triangulate(int vtkNotUsed (index),vtkIdList * ptIds,vtkPoints * pts)699 int vtkLine::Triangulate(int vtkNotUsed(index), vtkIdList *ptIds,
700 vtkPoints *pts)
701 {
702 pts->Reset();
703 ptIds->Reset();
704
705 ptIds->InsertId(0,this->PointIds->GetId(0));
706 pts->InsertPoint(0,this->Points->GetPoint(0));
707
708 ptIds->InsertId(1,this->PointIds->GetId(1));
709 pts->InsertPoint(1,this->Points->GetPoint(1));
710
711 return 1;
712 }
713
714 //----------------------------------------------------------------------------
Derivatives(int vtkNotUsed (subId),const double vtkNotUsed (pcoords)[3],const double * values,int dim,double * derivs)715 void vtkLine::Derivatives(int vtkNotUsed(subId),
716 const double vtkNotUsed(pcoords)[3],
717 const double *values,
718 int dim, double *derivs)
719 {
720 double x0[3], x1[3], deltaX[3];
721 int i, j;
722
723 this->Points->GetPoint(0, x0);
724 this->Points->GetPoint(1, x1);
725
726 for (i=0; i<3; i++)
727 {
728 deltaX[i] = x1[i] - x0[i];
729 }
730 for (i=0; i<dim; i++)
731 {
732 for (j=0; j<3; j++)
733 {
734 if ( deltaX[j] != 0 )
735 {
736 derivs[3*i+j] = (values[i+dim] - values[i]) / deltaX[j];
737 }
738 else
739 {
740 derivs[3*i+j] =0;
741 }
742 }
743 }
744 }
745
746
747 //----------------------------------------------------------------------------
748 // support line clipping
749 namespace { //required so we don't violate ODR
750 typedef int LINE_LIST;
751 typedef struct {
752 LINE_LIST lines[2];
753 } LINE_CASES;
754
755 static LINE_CASES lineCases[] = {
756 {{ -1, -1}}, // 0
757 {{100, 1}}, // 1
758 {{ 0, 101}}, // 2
759 {{100, 101}}}; // 3
760 }
761
762 // Clip this line using scalar value provided. Like contouring, except
763 // that it cuts the line to produce other lines.
Clip(double value,vtkDataArray * cellScalars,vtkIncrementalPointLocator * locator,vtkCellArray * lines,vtkPointData * inPd,vtkPointData * outPd,vtkCellData * inCd,vtkIdType cellId,vtkCellData * outCd,int insideOut)764 void vtkLine::Clip(double value, vtkDataArray *cellScalars,
765 vtkIncrementalPointLocator *locator, vtkCellArray *lines,
766 vtkPointData *inPd, vtkPointData *outPd,
767 vtkCellData *inCd, vtkIdType cellId, vtkCellData *outCd,
768 int insideOut)
769 {
770 static const int CASE_MASK[3] = {1,2};
771 LINE_CASES *lineCase;
772 int i, j, index, *vert, newCellId;
773 vtkIdType pts[3];
774 int vertexId;
775 double t, x1[3], x2[3], x[3];
776
777 // Build the case table
778 if ( insideOut )
779 {
780 for ( i=0, index = 0; i < 2; i++)
781 {
782 if (cellScalars->GetComponent(i,0) <= value)
783 {
784 index |= CASE_MASK[i];
785 }
786 }
787 }
788 else
789 {
790 for ( i=0, index = 0; i < 2; i++)
791 {
792 if (cellScalars->GetComponent(i,0) > value)
793 {
794 index |= CASE_MASK[i];
795 }
796 }
797 }
798
799 // Select the case based on the index and get the list of lines for this case
800 lineCase = lineCases + index;
801 vert = lineCase->lines;
802
803 // generate clipped line
804 if ( vert[0] > -1 )
805 {
806 for (i=0; i<2; i++) // insert line
807 {
808 // vertex exists, and need not be interpolated
809 if (vert[i] >= 100)
810 {
811 vertexId = vert[i] - 100;
812 this->Points->GetPoint(vertexId, x);
813 if ( locator->InsertUniquePoint(x, pts[i]) )
814 {
815 outPd->CopyData(inPd,this->PointIds->GetId(vertexId),pts[i]);
816 }
817 }
818
819 else //new vertex, interpolate
820 {
821 t = (value - cellScalars->GetComponent(0,0)) /
822 (cellScalars->GetComponent(1,0) - cellScalars->GetComponent(0,0));
823
824 this->Points->GetPoint(0, x1);
825 this->Points->GetPoint(1, x2);
826 for (j=0; j<3; j++)
827 {
828 x[j] = x1[j] + t * (x2[j] - x1[j]);
829 }
830
831 if ( locator->InsertUniquePoint(x, pts[i]) )
832 {
833 vtkIdType p1 = this->PointIds->GetId(0);
834 vtkIdType p2 = this->PointIds->GetId(1);
835 outPd->InterpolateEdge(inPd,pts[i],p1,p2,t);
836 }
837 }
838 }
839 // check for degenerate lines
840 if ( pts[0] != pts[1] )
841 {
842 newCellId = lines->InsertNextCell(2,pts);
843 outCd->CopyData(inCd,cellId,newCellId);
844 }
845 }
846 }
847
848 //----------------------------------------------------------------------------
849 //
850 // Compute interpolation functions
851 //
InterpolationFunctions(const double pcoords[3],double weights[2])852 void vtkLine::InterpolationFunctions(const double pcoords[3], double weights[2])
853 {
854 weights[0] = 1.0 - pcoords[0];
855 weights[1] = pcoords[0];
856 }
857
858 //----------------------------------------------------------------------------
InterpolationDerivs(const double vtkNotUsed (pcoords)[3],double derivs[2])859 void vtkLine::InterpolationDerivs(const double vtkNotUsed(pcoords)[3], double derivs[2])
860 {
861 derivs[0] = -1.0;
862 derivs[1] = 1.0;
863 }
864
865 //----------------------------------------------------------------------------
866 static double vtkLineCellPCoords[6] = {0.0,0.0,0.0, 1.0,0.0,0.0};
GetParametricCoords()867 double *vtkLine::GetParametricCoords()
868 {
869 return vtkLineCellPCoords;
870 }
871
872 //----------------------------------------------------------------------------
PrintSelf(ostream & os,vtkIndent indent)873 void vtkLine::PrintSelf(ostream& os, vtkIndent indent)
874 {
875 this->Superclass::PrintSelf(os,indent);
876 }
877