1 /*=========================================================================
2
3 Program: Visualization Toolkit
4 Module: vtkPixel.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 "vtkPixel.h"
16
17 #include "vtkObjectFactory.h"
18 #include "vtkQuad.h"
19 #include "vtkTriangle.h"
20 #include "vtkPlane.h"
21 #include "vtkMath.h"
22 #include "vtkCellArray.h"
23 #include "vtkLine.h"
24 #include "vtkIncrementalPointLocator.h"
25 #include "vtkPointData.h"
26 #include "vtkCellData.h"
27 #include "vtkPoints.h"
28 #include "vtkMarchingSquaresLineCases.h"
29
30 vtkStandardNewMacro(vtkPixel);
31
32 //----------------------------------------------------------------------------
33 // Construct the pixel with four points.
vtkPixel()34 vtkPixel::vtkPixel()
35 {
36 int i;
37
38 this->Points->SetNumberOfPoints(4);
39 this->PointIds->SetNumberOfIds(4);
40 for (i = 0; i < 4; i++)
41 {
42 this->Points->SetPoint(i, 0.0, 0.0, 0.0);
43 }
44 for (i = 0; i < 4; i++)
45 {
46 this->PointIds->SetId(i,0);
47 }
48 this->Line = vtkLine::New();
49 }
50
51 //----------------------------------------------------------------------------
~vtkPixel()52 vtkPixel::~vtkPixel()
53 {
54 this->Line->Delete();
55 }
56
57 //----------------------------------------------------------------------------
EvaluatePosition(double x[3],double * closestPoint,int & subId,double pcoords[3],double & dist2,double * weights)58 int vtkPixel::EvaluatePosition(double x[3], double* closestPoint,
59 int& subId, double pcoords[3],
60 double& dist2, double *weights)
61 {
62 double pt1[3], pt2[3], pt3[3];
63 int i;
64 double p[3], p21[3], p31[3], cp[3];
65 double l21, l31, n[3];
66
67 subId = 0;
68 pcoords[2] = 0.0;
69
70 // Get normal for pixel
71 //
72 this->Points->GetPoint(0, pt1);
73 this->Points->GetPoint(1, pt2);
74 this->Points->GetPoint(2, pt3);
75
76 vtkTriangle::ComputeNormal (pt1, pt2, pt3, n);
77
78 // Project point to plane
79 //
80 vtkPlane::ProjectPoint(x,pt1,n,cp);
81
82 for (i=0; i<3; i++)
83 {
84 p21[i] = pt2[i] - pt1[i];
85 p31[i] = pt3[i] - pt1[i];
86 p[i] = x[i] - pt1[i];
87 }
88
89 if ( (l21=vtkMath::Norm(p21)) == 0.0 )
90 {
91 l21 = 1.0;
92 }
93 if ( (l31=vtkMath::Norm(p31)) == 0.0 )
94 {
95 l31 = 1.0;
96 }
97
98 pcoords[0] = vtkMath::Dot(p21,p) / (l21*l21);
99 pcoords[1] = vtkMath::Dot(p31,p) / (l31*l31);
100
101 this->InterpolationFunctions(pcoords, weights);
102
103 if ( pcoords[0] >= 0.0 && pcoords[0] <= 1.0 &&
104 pcoords[1] >= 0.0 && pcoords[1] <= 1.0 )
105 {
106 if (closestPoint)
107 {
108 closestPoint[0] = cp[0];
109 closestPoint[1] = cp[1];
110 closestPoint[2] = cp[2];
111 dist2 =
112 vtkMath::Distance2BetweenPoints(closestPoint,x); //projection distance
113 }
114 return 1;
115 }
116 else
117 {
118 double pc[3], w[4];
119 if (closestPoint)
120 {
121 for (i=0; i<2; i++)
122 {
123 if (pcoords[i] < 0.0)
124 {
125 pc[i] = 0.0;
126 }
127 else if (pcoords[i] > 1.0)
128 {
129 pc[i] = 1.0;
130 }
131 else
132 {
133 pc[i] = pcoords[i];
134 }
135 }
136 this->EvaluateLocation(subId, pc, closestPoint,
137 static_cast<double *>(w));
138 dist2 = vtkMath::Distance2BetweenPoints(closestPoint,x);
139 }
140 return 0;
141 }
142 }
143
144 //----------------------------------------------------------------------------
EvaluateLocation(int & subId,double pcoords[3],double x[3],double * weights)145 void vtkPixel::EvaluateLocation(int& subId, double pcoords[3], double x[3],
146 double *weights)
147 {
148 double pt1[3], pt2[3], pt3[3];
149 int i;
150
151 subId = 0;
152
153 this->Points->GetPoint(0, pt1);
154 this->Points->GetPoint(1, pt2);
155 this->Points->GetPoint(2, pt3);
156
157 for (i=0; i<3; i++)
158 {
159 x[i] = pt1[i] + pcoords[0]*(pt2[i] - pt1[i]) +
160 pcoords[1]*(pt3[i] - pt1[i]);
161 }
162
163 this->InterpolationFunctions(pcoords, weights);
164 }
165
166 //----------------------------------------------------------------------------
CellBoundary(int vtkNotUsed (subId),double pcoords[3],vtkIdList * pts)167 int vtkPixel::CellBoundary(int vtkNotUsed(subId), double pcoords[3], vtkIdList *pts)
168 {
169 double t1=pcoords[0]-pcoords[1];
170 double t2=1.0-pcoords[0]-pcoords[1];
171
172 pts->SetNumberOfIds(2);
173
174 // compare against two lines in parametric space that divide element
175 // into four pieces.
176 if ( t1 >= 0.0 && t2 >= 0.0 )
177 {
178 pts->SetId(0,this->PointIds->GetId(0));
179 pts->SetId(1,this->PointIds->GetId(1));
180 }
181
182 else if ( t1 >= 0.0 && t2 < 0.0 )
183 {
184 pts->SetId(0,this->PointIds->GetId(1));
185 pts->SetId(1,this->PointIds->GetId(3));
186 }
187
188 else if ( t1 < 0.0 && t2 < 0.0 )
189 {
190 pts->SetId(0,this->PointIds->GetId(3));
191 pts->SetId(1,this->PointIds->GetId(2));
192 }
193
194 else //( t1 < 0.0 && t2 >= 0.0 )
195 {
196 pts->SetId(0,this->PointIds->GetId(2));
197 pts->SetId(1,this->PointIds->GetId(0));
198 }
199
200 if ( pcoords[0] < 0.0 || pcoords[0] > 1.0 ||
201 pcoords[1] < 0.0 || pcoords[1] > 1.0 )
202 {
203 return 0;
204 }
205 else
206 {
207 return 1;
208 }
209 }
210
211 //----------------------------------------------------------------------------
212 //
213 // Marching squares
214 //
215 #include "vtkMarchingSquaresCases.h"
216
217 static int edges[4][2] = { {0,1}, {1,3}, {2,3}, {0,2} };
218
Contour(double value,vtkDataArray * cellScalars,vtkIncrementalPointLocator * locator,vtkCellArray * vtkNotUsed (verts),vtkCellArray * lines,vtkCellArray * vtkNotUsed (polys),vtkPointData * inPd,vtkPointData * outPd,vtkCellData * inCd,vtkIdType cellId,vtkCellData * outCd)219 void vtkPixel::Contour(double value, vtkDataArray *cellScalars,
220 vtkIncrementalPointLocator *locator,
221 vtkCellArray *vtkNotUsed(verts),
222 vtkCellArray *lines,
223 vtkCellArray *vtkNotUsed(polys),
224 vtkPointData *inPd, vtkPointData *outPd,
225 vtkCellData *inCd, vtkIdType cellId, vtkCellData *outCd)
226 {
227 static int CASE_MASK[4] = {1,2,8,4}; //note differenceom quad!
228 vtkMarchingSquaresLineCases *lineCase;
229 EDGE_LIST *edge;
230 int i, j, index, *vert;
231 int newCellId;
232 vtkIdType pts[2];
233 double t, x1[3], x2[3], x[3];
234
235 // Build the case table
236 for ( i=0, index = 0; i < 4; i++)
237 {
238 if (cellScalars->GetComponent(i,0) >= value)
239 {
240 index |= CASE_MASK[i];
241 }
242 }
243
244 lineCase = vtkMarchingSquaresLineCases::GetCases() + index;
245 edge = lineCase->edges;
246
247 for ( ; edge[0] > -1; edge += 2 )
248 {
249 for (i=0; i<2; i++) // insert line
250 {
251 vert = edges[edge[i]];
252 t = (value - cellScalars->GetComponent(vert[0],0)) /
253 (cellScalars->GetComponent(vert[1],0) -
254 cellScalars->GetComponent(vert[0],0));
255 this->Points->GetPoint(vert[0], x1);
256 this->Points->GetPoint(vert[1], x2);
257 for (j=0; j<3; j++)
258 {
259 x[j] = x1[j] + t * (x2[j] - x1[j]);
260 }
261 if ( locator->InsertUniquePoint(x, pts[i]) )
262 {
263 if ( outPd )
264 {
265 int p1 = this->PointIds->GetId(vert[0]);
266 int p2 = this->PointIds->GetId(vert[1]);
267 outPd->InterpolateEdge(inPd,pts[i],p1,p2,t);
268 }
269 }
270 }
271 // check for degenerate line
272 if ( pts[0] != pts[1] )
273 {
274 newCellId = lines->InsertNextCell(2,pts);
275 outCd->CopyData(inCd,cellId,newCellId);
276 }
277 }
278 }
279
280 //----------------------------------------------------------------------------
GetEdge(int edgeId)281 vtkCell *vtkPixel::GetEdge(int edgeId)
282 {
283 int *verts;
284
285 verts = edges[edgeId];
286
287 // load point id's
288 this->Line->PointIds->SetId(0,this->PointIds->GetId(verts[0]));
289 this->Line->PointIds->SetId(1,this->PointIds->GetId(verts[1]));
290
291 // load coordinates
292 this->Line->Points->SetPoint(0,this->Points->GetPoint(verts[0]));
293 this->Line->Points->SetPoint(1,this->Points->GetPoint(verts[1]));
294
295 return this->Line;
296 }
297
298 //----------------------------------------------------------------------------
299 //
300 // Compute interpolation functions (similar but different than Quad interpolation
301 // functions)
302 //
InterpolationFunctions(double pcoords[3],double sf[4])303 void vtkPixel::InterpolationFunctions(double pcoords[3], double sf[4])
304 {
305 double rm, sm;
306
307 rm = 1. - pcoords[0];
308 sm = 1. - pcoords[1];
309
310 sf[0] = rm * sm;
311 sf[1] = pcoords[0] * sm;
312 sf[2] = rm * pcoords[1];
313 sf[3] = pcoords[0] * pcoords[1];
314 }
315 //----------------------------------------------------------------------------
316 //
317 // Compute derivatives of interpolation functions.
318 //
InterpolationDerivs(double pcoords[3],double derivs[8])319 void vtkPixel::InterpolationDerivs(double pcoords[3], double derivs[8])
320 {
321 double rm, sm;
322
323 rm = 1. - pcoords[0];
324 sm = 1. - pcoords[1];
325
326 // r derivatives
327 derivs[0] = -sm;
328 derivs[1] = sm;
329 derivs[2] = -pcoords[1];
330 derivs[3] = pcoords[1];
331
332 // s derivatives
333 derivs[4] = -rm;
334 derivs[5] = -pcoords[0];
335 derivs[6] = rm;
336 derivs[7] = pcoords[0];
337 }
338
339 //----------------------------------------------------------------------------
340 //
341 // Intersect plane; see whether point is inside.
342 //
IntersectWithLine(double p1[3],double p2[3],double tol,double & t,double x[3],double pcoords[3],int & subId)343 int vtkPixel::IntersectWithLine(double p1[3], double p2[3], double tol, double& t,
344 double x[3], double pcoords[3], int& subId)
345 {
346 double pt1[3], pt4[3], n[3];
347 double tol2 = tol*tol;
348 double closestPoint[3];
349 double dist2, weights[4];
350 int i;
351
352 subId = 0;
353 pcoords[0] = pcoords[1] = pcoords[2] = 0.0;
354 //
355 // Get normal for triangle
356 //
357 this->Points->GetPoint(0, pt1);
358 this->Points->GetPoint(3, pt4);
359
360 n[0] = n[1] = n[2] = 0.0;
361 for (i=0; i<3; i++)
362 {
363 if ( (pt4[i] - pt1[i]) <= 0.0 )
364 {
365 n[i] = 1.0;
366 break;
367 }
368 }
369 //
370 // Intersect plane of pixel with line
371 //
372 if ( ! vtkPlane::IntersectWithLine(p1,p2,n,pt1,t,x) )
373 {
374 return 0;
375 }
376 //
377 // Use evaluate position
378 //
379 if (this->EvaluatePosition(x, closestPoint, subId, pcoords, dist2, weights))
380 {
381 if ( dist2 <= tol2 )
382 {
383 return 1;
384 }
385 }
386
387 return 0;
388 }
389
390 //----------------------------------------------------------------------------
Triangulate(int index,vtkIdList * ptIds,vtkPoints * pts)391 int vtkPixel::Triangulate(int index, vtkIdList *ptIds, vtkPoints *pts)
392 {
393 pts->Reset();
394 ptIds->Reset();
395
396 if ( (index % 2) )
397 {
398 ptIds->InsertId(0,this->PointIds->GetId(0));
399 pts->InsertPoint(0,this->Points->GetPoint(0));
400 ptIds->InsertId(1,this->PointIds->GetId(1));
401 pts->InsertPoint(1,this->Points->GetPoint(1));
402 ptIds->InsertId(2,this->PointIds->GetId(2));
403 pts->InsertPoint(2,this->Points->GetPoint(2));
404
405 ptIds->InsertId(3,this->PointIds->GetId(1));
406 pts->InsertPoint(3,this->Points->GetPoint(1));
407 ptIds->InsertId(4,this->PointIds->GetId(3));
408 pts->InsertPoint(4,this->Points->GetPoint(3));
409 ptIds->InsertId(5,this->PointIds->GetId(2));
410 pts->InsertPoint(5,this->Points->GetPoint(2));
411 }
412 else
413 {
414 ptIds->InsertId(0,this->PointIds->GetId(0));
415 pts->InsertPoint(0,this->Points->GetPoint(0));
416 ptIds->InsertId(1,this->PointIds->GetId(1));
417 pts->InsertPoint(1,this->Points->GetPoint(1));
418 ptIds->InsertId(2,this->PointIds->GetId(3));
419 pts->InsertPoint(2,this->Points->GetPoint(3));
420
421 ptIds->InsertId(3,this->PointIds->GetId(0));
422 pts->InsertPoint(3,this->Points->GetPoint(0));
423 ptIds->InsertId(4,this->PointIds->GetId(3));
424 pts->InsertPoint(4,this->Points->GetPoint(3));
425 ptIds->InsertId(5,this->PointIds->GetId(2));
426 pts->InsertPoint(5,this->Points->GetPoint(2));
427 }
428
429 return 1;
430 }
431
432 //----------------------------------------------------------------------------
Derivatives(int vtkNotUsed (subId),double pcoords[3],double * values,int dim,double * derivs)433 void vtkPixel::Derivatives(int vtkNotUsed(subId),
434 double pcoords[3],
435 double *values,
436 int dim, double *derivs)
437 {
438 double functionDerivs[8], sum;
439 int i, j, k, plane, idx[2], jj;
440 double x0[3], x1[3], x2[3], x3[3], spacing[3];
441
442 this->Points->GetPoint(0, x0);
443 this->Points->GetPoint(1, x1);
444 this->Points->GetPoint(2, x2);
445 this->Points->GetPoint(3, x3);
446
447 //figure which plane this pixel is in
448 for (i=0; i < 3; i++)
449 {
450 spacing[i] = x3[i] - x0[i];
451 }
452
453 if ( spacing[0] > spacing[2] && spacing[1] > spacing[2] ) // z-plane
454 {
455 plane = 2;
456 idx[0] = 0; idx[1] = 1;
457 }
458 else if ( spacing[0] > spacing[1] && spacing[2] > spacing[1] ) // y-plane
459 {
460 plane = 1;
461 idx[0] = 0; idx[1] = 2;
462 }
463 else // x-plane
464 {
465 plane = 0;
466 idx[0] = 1; idx[1] = 2;
467 }
468
469 spacing[0] = x1[idx[0]] - x0[idx[0]];
470 spacing[1] = x2[idx[1]] - x0[idx[1]];
471
472 // get derivatives in r-s directions
473 this->InterpolationDerivs(pcoords, functionDerivs);
474
475 // since two of the x-y-z axes are aligned with r-s axes, only need to scale
476 // the derivative values by the data spacing.
477 for (k=0; k < dim; k++) //loop over values per vertex
478 {
479 for (jj=j=0; j < 3; j++) //loop over derivative directions
480 {
481 if ( j == plane ) // 0-derivate values in this direction
482 {
483 sum = 0.0;
484 }
485 else //compute derivatives
486 {
487 for (sum=0.0, i=0; i < 4; i++) //loop over interp. function derivatives
488 {
489 sum += functionDerivs[4*jj + i] * values[dim*i + k];
490 }
491 sum /= spacing[idx[jj++]];
492 }
493 derivs[3*k + j] = sum;
494 }
495 }
496 }
497
498 //----------------------------------------------------------------------------
499 // support pixel clipping
500 typedef int PIXEL_EDGE_LIST;
501 typedef struct {
502 PIXEL_EDGE_LIST edges[14];
503 } PIXEL_CASES;
504
505 static PIXEL_CASES pixelCases[] = {
506 {{ -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}}, // 0
507 {{ 3, 100, 0, 3, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}}, // 1
508 {{ 3, 101, 1, 0, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}}, // 2
509 {{ 4, 100, 101, 1, 3, -1, -1, -1, -1, -1, -1, -1, -1, -1}}, // 3
510 {{ 3, 103, 2, 1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}}, // 4
511 {{ 3, 100, 0, 3, 3, 103, 2, 1, 4, 0, 1, 2, 3, -1}}, // 5
512 {{ 4, 101, 103, 2, 0, -1, -1, -1, -1, -1, -1, -1, -1, -1}}, // 6
513 {{ 3, 100, 101, 3, 3, 101, 2, 3, 3, 101, 103, 2, -1, -1}}, // 7
514 {{ 3, 102, 3, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}}, // 8
515 {{ 4, 100, 0, 2, 102, -1, -1, -1, -1, -1, -1, -1, -1, -1}}, // 9
516 {{ 3, 101, 1, 0, 3, 102, 3, 2, 4, 0, 1, 2, 3, -1}}, // 10
517 {{ 3, 100, 101, 1, 3, 100, 1, 2, 3, 100, 2, 102, -1, -1}}, // 11
518 {{ 4, 103, 102, 3, 1, -1, -1, -1, -1, -1, -1, -1, -1, -1}}, // 12
519 {{ 3, 100, 0, 102, 3, 0, 1, 102, 3, 1, 103, 102, -1, -1}}, // 13
520 {{ 3, 0, 101, 103, 3, 0, 103, 3, 3, 103, 102, 3, -1, -1}}, // 14
521 {{ 4, 100, 101, 103, 102, -1, -1, -1, -1, -1, -1, -1, -1, -1}}, // 15
522 };
523
524 static PIXEL_CASES pixelCasesComplement[] = {
525 {{ -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}}, // 0
526 {{ 3, 100, 0, 3, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}}, // 1
527 {{ 3, 101, 1, 0, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}}, // 2
528 {{ 4, 100, 101, 1, 3, -1, -1, -1, -1, -1, -1, -1, -1, -1}}, // 3
529 {{ 3, 103, 2, 1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}}, // 4
530 {{ 3, 100, 0, 3, 3, 103, 2, 1, -1, -1, -1, -1, -1, -1}}, // 5
531 {{ 4, 101, 103, 2, 0, -1, -1, -1, -1, -1, -1, -1, -1, -1}}, // 6
532 {{ 3, 100, 101, 3, 3, 101, 2, 3, 3, 101, 103, 2, -1, -1}}, // 7
533 {{ 3, 102, 3, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}}, // 8
534 {{ 4, 100, 0, 2, 102, -1, -1, -1, -1, -1, -1, -1, -1, -1}}, // 9
535 {{ 3, 101, 1, 0, 3, 102, 3, 2, -1, -1, -1, -1, -1, -1}}, // 10
536 {{ 3, 100, 101, 1, 3, 100, 1, 2, 3, 100, 2, 102, -1, -1}}, // 11
537 {{ 4, 103, 102, 3, 1, -1, -1, -1, -1, -1, -1, -1, -1, -1}}, // 12
538 {{ 3, 100, 0, 102, 3, 0, 1, 102, 3, 1, 103, 102, -1, -1}}, // 13
539 {{ 3, 0, 101, 103, 3, 0, 103, 3, 3, 103, 102, 3, -1, -1}}, // 14
540 {{ 4, 100, 101, 103, 102, -1, -1, -1, -1, -1, -1, -1, -1, -1}}, // 15
541 };
542
543
544 //----------------------------------------------------------------------------
545 // Clip this pixel using scalar value provided. Like contouring, except
546 // that it cuts the pixel to produce 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)547 void vtkPixel::Clip(double value, vtkDataArray *cellScalars,
548 vtkIncrementalPointLocator *locator, vtkCellArray *polys,
549 vtkPointData *inPd, vtkPointData *outPd,
550 vtkCellData *inCd, vtkIdType cellId, vtkCellData *outCd,
551 int insideOut)
552 {
553 static int CASE_MASK[4] = {1,2,8,4}; //note difference from quad!
554 PIXEL_CASES *pixelCase;
555 PIXEL_EDGE_LIST *edge;
556 int i, j, index, *vert;
557 int e1, e2;
558 int newCellId;
559 vtkIdType pts[4];
560 int vertexId;
561 double t, x1[3], x2[3], x[3], deltaScalar;
562 double scalar0, scalar1, e1Scalar;
563
564 // Build the index into the case table
565 if ( insideOut )
566 {
567 for ( i=0, index = 0; i < 4; i++)
568 {
569 if (cellScalars->GetComponent(i,0) <= value)
570 {
571 index |= CASE_MASK[i];
572 }
573 }
574 // Select case based on the index and get the list of edges for this case
575 pixelCase = pixelCases + index;
576 }
577 else
578 {
579 for ( i=0, index = 0; i < 4; i++)
580 {
581 if (cellScalars->GetComponent(i,0) > value)
582 {
583 index |= CASE_MASK[i];
584 }
585 }
586 // Select case based on the index and get the list of edges for this case
587 pixelCase = pixelCasesComplement + index;
588 }
589
590 edge = pixelCase->edges;
591
592 // generate each pixel
593 for ( ; edge[0] > -1; edge += edge[0]+1 )
594 {
595 for (i=0; i < edge[0]; i++) // insert pixel or triangle
596 {
597 // vertex exists, and need not be interpolated
598 if (edge[i+1] >= 100)
599 {
600 vertexId = edge[i+1] - 100;
601 this->Points->GetPoint(vertexId, x);
602 if ( locator->InsertUniquePoint(x, pts[i]) )
603 {
604 outPd->CopyData(inPd,this->PointIds->GetId(vertexId),pts[i]);
605 }
606 }
607
608 else //new vertex, interpolate
609 {
610 vert = edges[edge[i+1]];
611
612 // calculate a preferred interpolation direction
613 scalar0 = cellScalars->GetComponent(vert[0],0);
614 scalar1 = cellScalars->GetComponent(vert[1],0);
615 deltaScalar = scalar1 - scalar0;
616
617 if (deltaScalar > 0)
618 {
619 e1 = vert[0]; e2 = vert[1];
620 e1Scalar = scalar0;
621 }
622 else
623 {
624 e1 = vert[1]; e2 = vert[0];
625 e1Scalar = scalar1;
626 deltaScalar = -deltaScalar;
627 }
628
629 // linear interpolation
630 if (deltaScalar == 0.0)
631 {
632 t = 0.0;
633 }
634 else
635 {
636 t = (value - e1Scalar) / deltaScalar;
637 }
638
639 this->Points->GetPoint(e1, x1);
640 this->Points->GetPoint(e2, x2);
641
642 for (j=0; j<3; j++)
643 {
644 x[j] = x1[j] + t * (x2[j] - x1[j]);
645 }
646
647 if ( locator->InsertUniquePoint(x, pts[i]) )
648 {
649 int p1 = this->PointIds->GetId(e1);
650 int p2 = this->PointIds->GetId(e2);
651 outPd->InterpolateEdge(inPd,pts[i],p1,p2,t);
652 }
653 }
654 }
655 // check for degenerate output
656 if ( edge[0] == 3 ) //i.e., a triangle
657 {
658 if (pts[0] == pts[1] || pts[0] == pts[2] || pts[1] == pts[2] )
659 {
660 continue;
661 }
662 }
663 else // a pixel
664 {
665 if ((pts[0] == pts[3] && pts[1] == pts[2]) ||
666 (pts[0] == pts[1] && pts[3] == pts[2]) )
667 {
668 continue;
669 }
670 }
671
672 newCellId = polys->InsertNextCell(edge[0],pts);
673 outCd->CopyData(inCd,cellId,newCellId);
674 }
675 }
676
677 //----------------------------------------------------------------------------
678 static double vtkPixelCellPCoords[12] = {0.0,0.0,0.0, 1.0,0.0,0.0,
679 0.0,1.0,0.0, 1.0,1.0,0.0};
680
GetParametricCoords()681 double *vtkPixel::GetParametricCoords()
682 {
683 return vtkPixelCellPCoords;
684 }
685
686 //----------------------------------------------------------------------------
PrintSelf(ostream & os,vtkIndent indent)687 void vtkPixel::PrintSelf(ostream& os, vtkIndent indent)
688 {
689 this->Superclass::PrintSelf(os,indent);
690
691 os << indent << "Line:\n";
692 this->Line->PrintSelf(os,indent.GetNextIndent());
693 }
694
695