1 /*=========================================================================
2
3 Program: Visualization Toolkit
4 Module: vtkVoxel.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 "vtkVoxel.h"
16
17 #include "vtkCellArray.h"
18 #include "vtkCellData.h"
19 #include "vtkLine.h"
20 #include "vtkMath.h"
21 #include "vtkObjectFactory.h"
22 #include "vtkPixel.h"
23 #include "vtkPointData.h"
24 #include "vtkIncrementalPointLocator.h"
25 #include "vtkPoints.h"
26 #include "vtkBox.h"
27 #include "vtkMarchingCubesTriangleCases.h"
28
29 vtkStandardNewMacro(vtkVoxel);
30
31 //----------------------------------------------------------------------------
32 // Construct the voxel with eight points.
vtkVoxel()33 vtkVoxel::vtkVoxel()
34 {
35 int i;
36
37 this->Points->SetNumberOfPoints(8);
38 this->PointIds->SetNumberOfIds(8);
39 for (i = 0; i < 8; i++)
40 {
41 this->Points->SetPoint(i, 0.0, 0.0, 0.0);
42 }
43 for (i = 0; i < 8; i++)
44 {
45 this->PointIds->SetId(i,0);
46 }
47 this->Line = 0;
48 this->Pixel = 0;
49 }
50
51 //----------------------------------------------------------------------------
~vtkVoxel()52 vtkVoxel::~vtkVoxel()
53 {
54 if (this->Line)
55 {
56 this->Line->Delete();
57 }
58 if (this->Pixel)
59 {
60 this->Pixel->Delete();
61 }
62 }
63
64 //----------------------------------------------------------------------------
EvaluatePosition(double x[3],double * closestPoint,int & subId,double pcoords[3],double & dist2,double * weights)65 int vtkVoxel::EvaluatePosition(double x[3], double* closestPoint,
66 int& subId, double pcoords[3],
67 double& dist2, double *weights)
68 {
69 double pt1[3], pt2[3], pt3[3], pt4[3];
70 int i;
71
72 subId = 0;
73 //
74 // Get coordinate system
75 //
76 this->Points->GetPoint(0, pt1);
77 this->Points->GetPoint(1, pt2);
78 this->Points->GetPoint(2, pt3);
79 this->Points->GetPoint(4, pt4);
80 //
81 // Develop parametric coordinates
82 //
83 pcoords[0] = (x[0] - pt1[0]) / (pt2[0] - pt1[0]);
84 pcoords[1] = (x[1] - pt1[1]) / (pt3[1] - pt1[1]);
85 pcoords[2] = (x[2] - pt1[2]) / (pt4[2] - pt1[2]);
86
87 if ( pcoords[0] >= 0.0 && pcoords[0] <= 1.0 &&
88 pcoords[1] >= 0.0 && pcoords[1] <= 1.0 &&
89 pcoords[2] >= 0.0 && pcoords[2] <= 1.0 )
90 {
91 if (closestPoint)
92 {
93 closestPoint[0] = x[0]; closestPoint[1] = x[1]; closestPoint[2] = x[2];
94 }
95 dist2 = 0.0; // inside voxel
96 this->InterpolationFunctions(pcoords,weights);
97 return 1;
98 }
99 else
100 {
101 double pc[3], w[8];
102 if (closestPoint)
103 {
104 for (i=0; i<3; i++)
105 {
106 if (pcoords[i] < 0.0)
107 {
108 pc[i] = 0.0;
109 }
110 else if (pcoords[i] > 1.0)
111 {
112 pc[i] = 1.0;
113 }
114 else
115 {
116 pc[i] = pcoords[i];
117 }
118 }
119 this->EvaluateLocation(subId, pc, closestPoint,
120 static_cast<double *>(w));
121 dist2 = vtkMath::Distance2BetweenPoints(closestPoint,x);
122 }
123 return 0;
124 }
125 }
126
127 //----------------------------------------------------------------------------
EvaluateLocation(int & vtkNotUsed (subId),double pcoords[3],double x[3],double * weights)128 void vtkVoxel::EvaluateLocation(int& vtkNotUsed(subId), double pcoords[3],
129 double x[3], double *weights)
130 {
131 double pt1[3], pt2[3], pt3[3], pt4[3];
132 int i;
133
134 this->Points->GetPoint(0, pt1);
135 this->Points->GetPoint(1, pt2);
136 this->Points->GetPoint(2, pt3);
137 this->Points->GetPoint(4, pt4);
138
139 for (i=0; i<3; i++)
140 {
141 x[i] = pt1[i] + pcoords[0]*(pt2[i] - pt1[i]) +
142 pcoords[1]*(pt3[i] - pt1[i]) +
143 pcoords[2]*(pt4[i] - pt1[i]);
144 }
145
146 this->InterpolationFunctions(pcoords,weights);
147 }
148
149 //----------------------------------------------------------------------------
150 //
151 // Compute Interpolation functions
152 //
InterpolationFunctions(double pcoords[3],double sf[8])153 void vtkVoxel::InterpolationFunctions(double pcoords[3], double sf[8])
154 {
155 double rm, sm, tm;
156
157 double r = pcoords[0], s = pcoords[1], t = pcoords[2];
158
159 rm = 1. - r;
160 sm = 1. - s;
161 tm = 1. - t;
162
163 sf[0] = rm * sm * tm;
164 sf[1] = r * sm * tm;
165 sf[2] = rm * s * tm;
166 sf[3] = r * s * tm;
167 sf[4] = rm * sm * t;
168 sf[5] = r * sm * t;
169 sf[6] = rm * s * t;
170 sf[7] = r * s * t;
171 }
172
173 //----------------------------------------------------------------------------
InterpolationDerivs(double pcoords[3],double derivs[24])174 void vtkVoxel::InterpolationDerivs(double pcoords[3], double derivs[24])
175 {
176 double rm, sm, tm;
177
178 rm = 1. - pcoords[0];
179 sm = 1. - pcoords[1];
180 tm = 1. - pcoords[2];
181
182 // r derivatives
183 derivs[0] = -sm*tm;
184 derivs[1] = sm*tm;
185 derivs[2] = -pcoords[1]*tm;
186 derivs[3] = pcoords[1]*tm;
187 derivs[4] = -sm*pcoords[2];
188 derivs[5] = sm*pcoords[2];
189 derivs[6] = -pcoords[1]*pcoords[2];
190 derivs[7] = pcoords[1]*pcoords[2];
191
192 // s derivatives
193 derivs[8] = -rm*tm;
194 derivs[9] = -pcoords[0]*tm;
195 derivs[10] = rm*tm;
196 derivs[11] = pcoords[0]*tm;
197 derivs[12] = -rm*pcoords[2];
198 derivs[13] = -pcoords[0]*pcoords[2];
199 derivs[14] = rm*pcoords[2];
200 derivs[15] = pcoords[0]*pcoords[2];
201
202 // t derivatives
203 derivs[16] = -rm*sm;
204 derivs[17] = -pcoords[0]*sm;
205 derivs[18] = -rm*pcoords[1];
206 derivs[19] = -pcoords[0]*pcoords[1];
207 derivs[20] = rm*sm;
208 derivs[21] = pcoords[0]*sm;
209 derivs[22] = rm*pcoords[1];
210 derivs[23] = pcoords[0]*pcoords[1];
211 }
212
213 //----------------------------------------------------------------------------
CellBoundary(int vtkNotUsed (subId),double pcoords[3],vtkIdList * pts)214 int vtkVoxel::CellBoundary(int vtkNotUsed(subId), double pcoords[3],
215 vtkIdList *pts)
216 {
217 double t1=pcoords[0]-pcoords[1];
218 double t2=1.0-pcoords[0]-pcoords[1];
219 double t3=pcoords[1]-pcoords[2];
220 double t4=1.0-pcoords[1]-pcoords[2];
221 double t5=pcoords[2]-pcoords[0];
222 double t6=1.0-pcoords[2]-pcoords[0];
223
224 pts->SetNumberOfIds(4);
225
226 // compare against six planes in parametric space that divide element
227 // into six pieces.
228 if ( t3 >= 0.0 && t4 >= 0.0 && t5 < 0.0 && t6 >= 0.0 )
229 {
230 pts->SetId(0,this->PointIds->GetId(0));
231 pts->SetId(1,this->PointIds->GetId(1));
232 pts->SetId(2,this->PointIds->GetId(3));
233 pts->SetId(3,this->PointIds->GetId(2));
234 }
235
236 else if ( t1 >= 0.0 && t2 < 0.0 && t5 < 0.0 && t6 < 0.0 )
237 {
238 pts->SetId(0,this->PointIds->GetId(1));
239 pts->SetId(1,this->PointIds->GetId(3));
240 pts->SetId(2,this->PointIds->GetId(7));
241 pts->SetId(3,this->PointIds->GetId(5));
242 }
243
244 else if ( t1 >= 0.0 && t2 >= 0.0 && t3 < 0.0 && t4 >= 0.0 )
245 {
246 pts->SetId(0,this->PointIds->GetId(0));
247 pts->SetId(1,this->PointIds->GetId(1));
248 pts->SetId(2,this->PointIds->GetId(5));
249 pts->SetId(3,this->PointIds->GetId(4));
250 }
251
252 else if ( t3 < 0.0 && t4 < 0.0 && t5 >= 0.0 && t6 < 0.0 )
253 {
254 pts->SetId(0,this->PointIds->GetId(4));
255 pts->SetId(1,this->PointIds->GetId(5));
256 pts->SetId(2,this->PointIds->GetId(7));
257 pts->SetId(3,this->PointIds->GetId(6));
258 }
259
260 else if ( t1 < 0.0 && t2 >= 0.0 && t5 >= 0.0 && t6 >= 0.0 )
261 {
262 pts->SetId(0,this->PointIds->GetId(0));
263 pts->SetId(1,this->PointIds->GetId(4));
264 pts->SetId(2,this->PointIds->GetId(6));
265 pts->SetId(3,this->PointIds->GetId(2));
266 }
267
268 else // if ( t1 < 0.0 && t2 < 0.0 && t3 >= 0.0 && t6 < 0.0 )
269 {
270 pts->SetId(0,this->PointIds->GetId(3));
271 pts->SetId(1,this->PointIds->GetId(2));
272 pts->SetId(2,this->PointIds->GetId(6));
273 pts->SetId(3,this->PointIds->GetId(7));
274 }
275
276 if ( pcoords[0] < 0.0 || pcoords[0] > 1.0 ||
277 pcoords[1] < 0.0 || pcoords[1] > 1.0 ||
278 pcoords[2] < 0.0 || pcoords[2] > 1.0 )
279 {
280 return 0;
281 }
282 else
283 {
284 return 1;
285 }
286 }
287
288 //----------------------------------------------------------------------------
289 static int edges[12][2] = { {0,1}, {1,3}, {2,3}, {0,2},
290 {4,5}, {5,7}, {6,7}, {4,6},
291 {0,4}, {1,5}, {2,6}, {3,7}};
292 // define in terms vtkPixel understands
293 static int faces[6][4] = { {2,0,6,4}, {1,3,5,7},
294 {0,1,4,5}, {3,2,7,6},
295 {1,0,3,2}, {4,5,6,7} };
296
297 //----------------------------------------------------------------------------
298 //
299 // Marching cubes case table
300 //
301 #include "vtkMarchingCubesCases.h"
302
Contour(double value,vtkDataArray * cellScalars,vtkIncrementalPointLocator * locator,vtkCellArray * verts,vtkCellArray * lines,vtkCellArray * polys,vtkPointData * inPd,vtkPointData * outPd,vtkCellData * inCd,vtkIdType cellId,vtkCellData * outCd)303 void vtkVoxel::Contour(double value, vtkDataArray *cellScalars,
304 vtkIncrementalPointLocator *locator,
305 vtkCellArray *verts,
306 vtkCellArray *lines,
307 vtkCellArray *polys,
308 vtkPointData *inPd, vtkPointData *outPd,
309 vtkCellData *inCd, vtkIdType cellId, vtkCellData *outCd)
310 {
311 static int CASE_MASK[8] = {1,2,4,8,16,32,64,128};
312 vtkMarchingCubesTriangleCases *triCase;
313 EDGE_LIST *edge;
314 int i, j, index, *vert;
315 static int vertMap[8] = { 0, 1, 3, 2, 4, 5, 7, 6 };
316 int newCellId;
317 vtkIdType pts[3];
318 double t, x1[3], x2[3], x[3];
319 vtkIdType offset = verts->GetNumberOfCells() + lines->GetNumberOfCells();
320
321 // Build the case table
322 for ( i=0, index = 0; i < 8; i++)
323 {
324 if (cellScalars->GetComponent(vertMap[i],0) >= value)
325 {
326 index |= CASE_MASK[i];
327 }
328 }
329
330 triCase = vtkMarchingCubesTriangleCases::GetCases() + index;
331 edge = triCase->edges;
332
333 for ( ; edge[0] > -1; edge += 3 )
334 {
335 for (i=0; i<3; i++) // insert triangle
336 {
337 vert = edges[edge[i]];
338 t = (value - cellScalars->GetComponent(vert[0],0)) /
339 (cellScalars->GetComponent(vert[1],0)
340 - cellScalars->GetComponent(vert[0],0));
341 this->Points->GetPoint(vert[0], x1);
342 this->Points->GetPoint(vert[1], x2);
343 for (j=0; j<3; j++)
344 {
345 x[j] = x1[j] + t * (x2[j] - x1[j]);
346 }
347 if ( locator->InsertUniquePoint(x, pts[i]) )
348 {
349 if ( outPd )
350 {
351 int p1 = this->PointIds->GetId(vert[0]);
352 int p2 = this->PointIds->GetId(vert[1]);
353 outPd->InterpolateEdge(inPd,pts[i],p1,p2,t);
354 }
355 }
356 }
357 // check for degenerate triangle
358 if ( pts[0] != pts[1] &&
359 pts[0] != pts[2] &&
360 pts[1] != pts[2] )
361 {
362 newCellId = offset + polys->InsertNextCell(3,pts);
363 outCd->CopyData(inCd,cellId,newCellId);
364 }
365 }
366 }
367
368
369 //----------------------------------------------------------------------------
GetEdgeArray(int edgeId)370 int *vtkVoxel::GetEdgeArray(int edgeId)
371 {
372 return edges[edgeId];
373 }
374
375 //----------------------------------------------------------------------------
GetEdge(int edgeId)376 vtkCell *vtkVoxel::GetEdge(int edgeId)
377 {
378 if (!this->Line)
379 {
380 this->Line = vtkLine::New();
381 }
382
383 int *verts;
384
385 verts = edges[edgeId];
386
387 // load point id's
388 this->Line->PointIds->SetId(0,this->PointIds->GetId(verts[0]));
389 this->Line->PointIds->SetId(1,this->PointIds->GetId(verts[1]));
390
391 // load coordinates
392 this->Line->Points->SetPoint(0,this->Points->GetPoint(verts[0]));
393 this->Line->Points->SetPoint(1,this->Points->GetPoint(verts[1]));
394
395 return this->Line;
396 }
397
398 //----------------------------------------------------------------------------
GetFaceArray(int faceId)399 int *vtkVoxel::GetFaceArray(int faceId)
400 {
401 return faces[faceId];
402 }
403
404 //----------------------------------------------------------------------------
GetFace(int faceId)405 vtkCell *vtkVoxel::GetFace(int faceId)
406 {
407 if (!this->Pixel)
408 {
409 this->Pixel = vtkPixel::New();
410 }
411
412 int *verts, i;
413
414 verts = faces[faceId];
415
416 for (i=0; i<4; i++)
417 {
418 this->Pixel->PointIds->SetId(i,this->PointIds->GetId(verts[i]));
419 this->Pixel->Points->SetPoint(i,this->Points->GetPoint(verts[i]));
420 }
421
422 return this->Pixel;
423 }
424
425 //----------------------------------------------------------------------------
426 //
427 // Intersect voxel with line using "bounding box" intersection.
428 //
IntersectWithLine(double p1[3],double p2[3],double vtkNotUsed (tol),double & t,double x[3],double pcoords[3],int & subId)429 int vtkVoxel::IntersectWithLine(double p1[3], double p2[3],
430 double vtkNotUsed(tol),
431 double& t, double x[3],
432 double pcoords[3], int& subId)
433 {
434 double minPt[3], maxPt[3];
435 double bounds[6];
436 double p21[3];
437 int i;
438
439 subId = 0;
440
441 this->Points->GetPoint(0, minPt);
442 this->Points->GetPoint(7, maxPt);
443
444 for (i=0; i<3; i++)
445 {
446 p21[i] = p2[i] - p1[i];
447 bounds[2*i] = minPt[i];
448 bounds[2*i+1] = maxPt[i];
449 }
450
451 if ( ! vtkBox::IntersectBox(bounds, p1, p21, x, t) )
452 {
453 return 0;
454 }
455
456 //
457 // Evaluate intersection
458 //
459 for (i=0; i<3; i++)
460 {
461 pcoords[i] = (x[i] - minPt[i]) / (maxPt[i] - minPt[i]);
462 }
463
464 return 1;
465 }
466
467 //----------------------------------------------------------------------------
Triangulate(int index,vtkIdList * ptIds,vtkPoints * pts)468 int vtkVoxel::Triangulate(int index, vtkIdList *ptIds, vtkPoints *pts)
469 {
470 int p[4], i;
471
472 ptIds->Reset();
473 pts->Reset();
474 //
475 // Create five tetrahedron. Triangulation varies depending upon index. This
476 // is necessary to insure compatible voxel triangulations.
477 //
478 if ( (index % 2) )
479 {
480 p[0] = 0; p[1] = 1; p[2] = 2; p[3] = 4;
481 for ( i=0; i < 4; i++ )
482 {
483 ptIds->InsertNextId(this->PointIds->GetId(p[i]));
484 pts->InsertNextPoint(this->Points->GetPoint(p[i]));
485 }
486
487 p[0] = 1; p[1] = 4; p[2] = 5; p[3] = 7;
488 for ( i=0; i < 4; i++ )
489 {
490 ptIds->InsertNextId(this->PointIds->GetId(p[i]));
491 pts->InsertNextPoint(this->Points->GetPoint(p[i]));
492 }
493
494 p[0] = 1; p[1] = 4; p[2] = 7; p[3] = 2;
495 for ( i=0; i < 4; i++ )
496 {
497 ptIds->InsertNextId(this->PointIds->GetId(p[i]));
498 pts->InsertNextPoint(this->Points->GetPoint(p[i]));
499 }
500
501 p[0] = 1; p[1] = 2; p[2] = 7; p[3] = 3;
502 for ( i=0; i < 4; i++ )
503 {
504 ptIds->InsertNextId(this->PointIds->GetId(p[i]));
505 pts->InsertNextPoint(this->Points->GetPoint(p[i]));
506 }
507
508 p[0] = 2; p[1] = 7; p[2] = 6; p[3] = 4;
509 for ( i=0; i < 4; i++ )
510 {
511 ptIds->InsertNextId(this->PointIds->GetId(p[i]));
512 pts->InsertNextPoint(this->Points->GetPoint(p[i]));
513 }
514 }
515 else
516 {
517 p[0] = 3; p[1] = 1; p[2] = 5; p[3] = 0;
518 for ( i=0; i < 4; i++ )
519 {
520 ptIds->InsertNextId(this->PointIds->GetId(p[i]));
521 pts->InsertNextPoint(this->Points->GetPoint(p[i]));
522 }
523
524 p[0] = 0; p[1] = 3; p[2] = 2; p[3] = 6;
525 for ( i=0; i < 4; i++ )
526 {
527 ptIds->InsertNextId(this->PointIds->GetId(p[i]));
528 pts->InsertNextPoint(this->Points->GetPoint(p[i]));
529 }
530
531 p[0] = 3; p[1] = 5; p[2] = 7; p[3] = 6;
532 for ( i=0; i < 4; i++ )
533 {
534 ptIds->InsertNextId(this->PointIds->GetId(p[i]));
535 pts->InsertNextPoint(this->Points->GetPoint(p[i]));
536 }
537
538 p[0] = 0; p[1] = 6; p[2] = 4; p[3] = 5;
539 for ( i=0; i < 4; i++ )
540 {
541 ptIds->InsertNextId(this->PointIds->GetId(p[i]));
542 pts->InsertNextPoint(this->Points->GetPoint(p[i]));
543 }
544
545 p[0] = 0; p[1] = 3; p[2] = 6; p[3] = 5;
546 for ( i=0; i < 4; i++ )
547 {
548 ptIds->InsertNextId(this->PointIds->GetId(p[i]));
549 pts->InsertNextPoint(this->Points->GetPoint(p[i]));
550 }
551 }
552
553 return 1;
554 }
555
556 //----------------------------------------------------------------------------
Derivatives(int vtkNotUsed (subId),double pcoords[3],double * values,int dim,double * derivs)557 void vtkVoxel::Derivatives(int vtkNotUsed(subId), double pcoords[3],
558 double *values, int dim, double *derivs)
559 {
560 double functionDerivs[24], sum;
561 int i, j, k;
562 double x0[3], x1[3], x2[3], x4[3], spacing[3];
563
564 this->Points->GetPoint(0, x0);
565 this->Points->GetPoint(1, x1);
566 spacing[0] = x1[0] - x0[0];
567
568 this->Points->GetPoint(2, x2);
569 spacing[1] = x2[1] - x0[1];
570
571 this->Points->GetPoint(4, x4);
572 spacing[2] = x4[2] - x0[2];
573
574 // get derivatives in r-s-t directions
575 this->InterpolationDerivs(pcoords, functionDerivs);
576
577 // since the x-y-z axes are aligned with r-s-t axes, only need to scale
578 // the derivative values by the data spacing.
579 for (k=0; k < dim; k++) //loop over values per vertex
580 {
581 for (j=0; j < 3; j++) //loop over derivative directions
582 {
583 for (sum=0.0, i=0; i < 8; i++) //loop over interp. function derivatives
584 {
585 sum += functionDerivs[8*j + i] * values[dim*i + k];
586 }
587 derivs[3*k + j] = sum / spacing[j];
588 }
589 }
590 }
591
592 //----------------------------------------------------------------------------
GetEdgePoints(int edgeId,int * & pts)593 void vtkVoxel::GetEdgePoints(int edgeId, int* &pts)
594 {
595 pts = this->GetEdgeArray(edgeId);
596 }
597
598 //----------------------------------------------------------------------------
GetFacePoints(int faceId,int * & pts)599 void vtkVoxel::GetFacePoints(int faceId, int* &pts)
600 {
601 pts = this->GetFaceArray(faceId);
602 }
603
604 static double vtkVoxelCellPCoords[24] = {0.0,0.0,0.0, 1.0,0.0,0.0,
605 0.0,1.0,0.0, 1.0,1.0,0.0,
606 0.0,0.0,1.0, 1.0,0.0,1.0,
607 0.0,1.0,1.0, 1.0,1.0,1.0};
608
609 //----------------------------------------------------------------------------
GetParametricCoords()610 double *vtkVoxel::GetParametricCoords()
611 {
612 return vtkVoxelCellPCoords;
613 }
614
615 //----------------------------------------------------------------------------
PrintSelf(ostream & os,vtkIndent indent)616 void vtkVoxel::PrintSelf(ostream& os, vtkIndent indent)
617 {
618 this->Superclass::PrintSelf(os,indent);
619
620 os << indent << "Line:\n";
621 if (this->Line)
622 {
623 this->Line->PrintSelf(os,indent.GetNextIndent());
624 }
625 else
626 {
627 os << "None\n";
628 }
629 os << indent << "Pixel:\n";
630 if (this->Pixel)
631 {
632 this->Pixel->PrintSelf(os,indent.GetNextIndent());
633 }
634 else
635 {
636 os << "None\n";
637 }
638 }
639