1 /*=========================================================================
2
3 Program: Visualization Toolkit
4 Module: vtkSliceCubes.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 "vtkSliceCubes.h"
16
17 #include "vtkByteSwap.h"
18 #include "vtkCharArray.h"
19 #include "vtkDoubleArray.h"
20 #include "vtkFloatArray.h"
21 #include "vtkImageData.h"
22 #include "vtkIntArray.h"
23 #include "vtkLongArray.h"
24 #include "vtkMarchingCubesTriangleCases.h"
25 #include "vtkMath.h"
26 #include "vtkObjectFactory.h"
27 #include "vtkPointData.h"
28 #include "vtkShortArray.h"
29 #include "vtkUnsignedCharArray.h"
30 #include "vtkUnsignedIntArray.h"
31 #include "vtkUnsignedLongArray.h"
32 #include "vtkUnsignedShortArray.h"
33 #include "vtkVolumeReader.h"
34
35 vtkStandardNewMacro(vtkSliceCubes);
36
37 vtkCxxSetObjectMacro(vtkSliceCubes,Reader,vtkVolumeReader);
38
39 // Description:
40 // Construct with NULL reader, output FileName specification, and limits
41 // FileName.
vtkSliceCubes()42 vtkSliceCubes::vtkSliceCubes()
43 {
44 this->Reader = NULL;
45 this->FileName = NULL;
46 this->LimitsFileName = NULL;
47 this->Value = 0.0;
48 }
49
~vtkSliceCubes()50 vtkSliceCubes::~vtkSliceCubes()
51 {
52 delete [] this->FileName;
53 delete [] this->LimitsFileName;
54 this->SetReader(NULL);
55 }
56
57 // Description:
58 // Method causes object to read slices and generate isosurface.
Update()59 void vtkSliceCubes::Update()
60 {
61 this->Execute();
62 }
63
64 // Calculate the gradient using central difference.
65 // NOTE: We calculate the negative of the gradient for efficiency
66 template <class T>
ComputePointGradient(int i,int j,int k,int dims[3],double Spacing[3],double n[3],T * s0,T * s1,T * s2)67 void ComputePointGradient(int i, int j, int k, int dims[3],
68 double Spacing[3], double n[3], T *s0, T *s1, T *s2)
69 {
70 double sp, sm;
71
72 // x-direction
73 if ( i == 0 )
74 {
75 sp = s1[i+1 + j*dims[0]];
76 sm = s1[i + j*dims[0]];
77 n[0] = (sm - sp) / Spacing[0];
78 }
79 else if ( i == (dims[0]-1) )
80 {
81 sp = s1[i + j*dims[0]];
82 sm = s1[i-1 + j*dims[0]];
83 n[0] = (sm - sp) / Spacing[0];
84 }
85 else
86 {
87 sp = s1[i+1 + j*dims[0]];
88 sm = s1[i-1 + j*dims[0]];
89 n[0] = 0.5 * (sm - sp) / Spacing[0];
90 }
91
92 // y-direction
93 if ( j == 0 )
94 {
95 sp = s1[i + (j+1)*dims[0]];
96 sm = s1[i + j*dims[0]];
97 n[1] = (sm - sp) / Spacing[1];
98 }
99 else if ( j == (dims[1]-1) )
100 {
101 sp = s1[i + j*dims[0]];
102 sm = s1[i + (j-1)*dims[0]];
103 n[1] = (sm - sp) / Spacing[1];
104 }
105 else
106
107 {
108 sp = s1[i + (j+1)*dims[0]];
109 sm = s1[i + (j-1)*dims[0]];
110 n[1] = 0.5 * (sm - sp) / Spacing[1];
111 }
112
113 // z-direction
114 // z-direction
115 if ( k == 0 )
116 {
117 sp = s2[i + j*dims[0]];
118 sm = s1[i + j*dims[0]];
119 n[2] = (sm - sp) / Spacing[2];
120 }
121 else if ( k == (dims[2]-1) )
122 {
123 sp = s1[i + j*dims[0]];
124 sm = s0[i + j*dims[0]];
125 n[2] = (sm - sp) / Spacing[2];
126 }
127 else
128 {
129 sp = s2[i + j*dims[0]];
130 sm = s0[i + j*dims[0]];
131 n[2] = 0.5 * (sm - sp) / Spacing[2];
132 }
133 }
134
135 template <class T, class S>
vtkSliceCubesContour(T * slice,S * scalars,int imageRange[2],int dims[3],double origin[3],double Spacing[3],double value,double xmin[3],double xmax[3],FILE * outFP,vtkVolumeReader * reader,unsigned char debug)136 int vtkSliceCubesContour(T *slice, S *scalars, int imageRange[2], int dims[3],
137 double origin[3], double Spacing[3], double value,
138 double xmin[3], double xmax[3], FILE *outFP,
139 vtkVolumeReader *reader, unsigned char debug)
140 {
141 S *slice0scalars=NULL, *slice1scalars=NULL;
142 S *slice2scalars, *slice3scalars;
143 T *slice0, *slice1, *slice2, *slice3;
144 vtkImageData *sp;
145 vtkDoubleArray *doubleScalars=NULL;
146 int numTriangles=0, numComp = 0;
147 double s[8];
148 int i, j, k, idx, jOffset, ii, index, *vert, jj, sliceSize=0;
149 static int CASE_MASK[8] = {1,2,4,8,16,32,64,128};
150 vtkMarchingCubesTriangleCases *triCase, *triCases;
151 EDGE_LIST *edge;
152 double pts[8][3], grad[8][3];
153 double t, *x1, *x2, *n1, *n2;
154 double xp, yp, zp;
155 float point[6];
156 static int edges[12][2] = { {0,1}, {1,2}, {3,2}, {0,3},
157 {4,5}, {5,6}, {7,6}, {4,7},
158 {0,4}, {1,5}, {3,7}, {2,6}};
159
160 triCases = vtkMarchingCubesTriangleCases::GetCases();
161
162 if ( slice == NULL ) //have to do conversion to double slice-by-slice
163 {
164 sliceSize = dims[0] * dims[1];
165 doubleScalars = vtkDoubleArray::New();
166 doubleScalars->Allocate(sliceSize);
167 }
168
169 slice2scalars = scalars;
170 slice2scalars->Register(NULL);
171
172 if (debug) vtkGenericWarningMacro(<< " Slice# " << imageRange[0]);
173
174 if ( slice2scalars == NULL )
175 {
176 return 0;
177 }
178 if ( slice != NULL )
179 {
180 slice1 = slice2 = slice2scalars->GetPointer(0);
181 }
182 else
183 {//get as double
184 numComp = scalars->GetNumberOfComponents();
185 slice2scalars->GetData(0,sliceSize-1,0,numComp-1,doubleScalars);
186 slice1 = slice2 = (T *) doubleScalars->GetPointer(0);
187 }
188
189 sp = reader->GetImage(imageRange[0]+1);
190 slice3scalars = (S *) sp->GetPointData()->GetScalars();
191 slice3scalars->Register(NULL);
192 sp->Delete();
193
194 if (debug) vtkGenericWarningMacro(<< " Slice# " << imageRange[0]+1 );
195
196 if ( slice != NULL )
197 {
198 slice3 = slice3scalars->GetPointer(0);
199 }
200 else
201 {//get as double: cast is ok because this code is only executed for double type
202 slice3scalars->GetData(0,sliceSize-1,0,numComp-1,doubleScalars);
203 slice3 = (T *) doubleScalars->GetPointer(0);
204 }
205
206 if ( !slice2 || !slice3 )
207 {
208 vtkGenericWarningMacro(<< "Cannot allocate data!");
209 return 0;
210 }
211
212 // Generate triangles and normals from slices
213 for (k=0; k < (dims[2]-1); k++)
214 {
215
216 // swap things around
217 if ( slice0scalars != NULL )
218 {
219 slice0scalars->Delete();
220 }
221 slice0scalars = slice1scalars;
222 slice0 = slice1;
223 slice1scalars = slice2scalars;
224 slice1 = slice2;
225 slice2scalars = slice3scalars;
226 slice2 = slice3;
227 if ( k < (dims[2]-2) )
228 {
229 if (debug) vtkGenericWarningMacro(<< " Slice# " << imageRange[0]+k+2);
230 sp = reader->GetImage(imageRange[0]+k+2);
231 slice3scalars = (S *) sp->GetPointData()->GetScalars();
232 if ( slice3scalars == NULL )
233 {
234 vtkGenericWarningMacro(<< "Can't read all the requested slices");
235 goto PREMATURE_TERMINATION;
236 }
237 slice3scalars->Register(NULL);
238 sp->Delete();
239 if ( slice != NULL )
240 {
241 slice3 = slice3scalars->GetPointer(0);
242 }
243 else
244 {//get as double
245 slice3scalars->GetData(0,sliceSize-1,0,numComp-1,doubleScalars);
246 slice3 = (T *) doubleScalars->GetPointer(0);
247 }
248 }
249
250 pts[0][2] = origin[2] + k*Spacing[2];
251 zp = origin[2] + (k+1)*Spacing[2];
252 for ( j=0; j < (dims[1]-1); j++)
253 {
254 jOffset = j*dims[0];
255 pts[0][1] = origin[1] + j*Spacing[1];
256 yp = origin[1] + (j+1)*Spacing[1];
257 for ( i=0; i < (dims[0]-1); i++)
258 {
259 //get scalar values
260 idx = i + jOffset;
261 s[0] = slice1[idx];
262 s[1] = slice1[idx+1];
263 s[2] = slice1[idx+1 + dims[0]];
264 s[3] = slice1[idx + dims[0]];
265 s[4] = slice2[idx];
266 s[5] = slice2[idx+1];
267 s[6] = slice2[idx+1 + dims[0]];
268 s[7] = slice2[idx + dims[0]];
269
270 // Build the case table
271 for ( ii=0, index = 0; ii < 8; ii++)
272 {
273 if ( s[ii] >= value )
274 {
275 index |= CASE_MASK[ii];
276 }
277 }
278
279 if ( index == 0 || index == 255 ) // no surface
280 {
281 continue;
282 }
283 //create voxel points
284 pts[0][0] = origin[0] + i*Spacing[0];
285 xp = origin[0] + (i+1)*Spacing[0];
286
287 pts[1][0] = xp;
288 pts[1][1] = pts[0][1];
289 pts[1][2] = pts[0][2];
290
291 pts[2][0] = xp;
292 pts[2][1] = yp;
293 pts[2][2] = pts[0][2];
294
295 pts[3][0] = pts[0][0];
296 pts[3][1] = yp;
297 pts[3][2] = pts[0][2];
298
299 pts[4][0] = pts[0][0];
300 pts[4][1] = pts[0][1];
301 pts[4][2] = zp;
302
303 pts[5][0] = xp;
304 pts[5][1] = pts[0][1];
305 pts[5][2] = zp;
306
307 pts[6][0] = xp;
308 pts[6][1] = yp;
309 pts[6][2] = zp;
310
311 pts[7][0] = pts[0][0];
312 pts[7][1] = yp;
313 pts[7][2] = zp;
314
315 //create gradients
316 ComputePointGradient(i,j, k, dims, Spacing, grad[0],
317 slice0, slice1, slice2);
318 ComputePointGradient(i+1,j, k, dims, Spacing, grad[1],
319 slice0, slice1, slice2);
320 ComputePointGradient(i+1,j+1, k, dims, Spacing, grad[2],
321 slice0, slice1, slice2);
322 ComputePointGradient(i,j+1, k, dims, Spacing, grad[3],
323 slice0, slice1, slice2);
324 ComputePointGradient(i,j, k+1, dims, Spacing, grad[4],
325 slice1, slice2, slice3);
326 ComputePointGradient(i+1,j, k+1, dims, Spacing, grad[5],
327 slice1, slice2, slice3);
328 ComputePointGradient(i+1,j+1, k+1, dims, Spacing, grad[6],
329 slice1, slice2, slice3);
330 ComputePointGradient(i,j+1, k+1, dims, Spacing, grad[7],
331 slice1, slice2, slice3);
332
333 triCase = triCases + index;
334 edge = triCase->edges;
335
336 for ( ; edge[0] > -1; edge += 3 )
337 {
338 for (ii=0; ii<3; ii++) //insert triangle
339 {
340 vert = edges[edge[ii]];
341 t = (value - s[vert[0]]) / (s[vert[1]] - s[vert[0]]);
342 x1 = pts[vert[0]];
343 x2 = pts[vert[1]];
344 n1 = grad[vert[0]];
345 n2 = grad[vert[1]];
346 for (jj=0; jj<3; jj++)
347 {
348 point[jj] = x1[jj] + t * (x2[jj] - x1[jj]);
349 point[jj+3] = n1[jj] + t * (n2[jj] - n1[jj]);
350 if (point[jj] < xmin[jj] )
351 {
352 xmin[jj] = point[jj];
353 }
354 if (point[jj] > xmax[jj] )
355 {
356 xmax[jj] = point[jj];
357 }
358 }
359 vtkMath::Normalize(point+3);
360 // swap bytes if necessary
361 bool status=vtkByteSwap::SwapWrite4BERange(point,6,outFP);
362 if(!status)
363 {
364 vtkGenericWarningMacro(<< "SwapWrite failed!");
365 }
366 }
367 numTriangles++;
368 }//for each triangle
369 }//for i
370 }//for j
371 }//for k
372
373 // Close things down
374 PREMATURE_TERMINATION:
375
376 fclose(outFP);
377 if ( slice == NULL )
378 {
379 doubleScalars->Delete();
380 }
381 if (slice0scalars && slice0scalars != slice1scalars)
382 {
383 slice0scalars->Delete();
384 }
385 if (slice3scalars && slice3scalars != slice2scalars)
386 {
387 slice3scalars->Delete();
388 }
389 if (slice1scalars)
390 {
391 slice1scalars->Delete();
392 }
393 slice2scalars->Delete();
394
395 return numTriangles;
396 }
397
Execute()398 void vtkSliceCubes::Execute()
399 {
400 FILE *outFP;
401 vtkImageData *tempStructPts;
402 vtkDataArray *inScalars;
403 int dims[3], imageRange[2];
404 double xmin[3], xmax[3];
405 double origin[3], Spacing[3];
406
407 // check input/initalize
408 vtkDebugMacro(<< "Executing slice cubes");
409 if ( this->Reader == NULL )
410 {
411 vtkErrorMacro(<<"No reader specified...can't generate isosurface");
412 return;
413 }
414
415 if ( this->FileName == NULL )
416 {
417 vtkErrorMacro(<<"No FileName specified...can't output isosurface");
418 return;
419 }
420
421 if ( (outFP = fopen(this->FileName, "wb")) == NULL )
422 {
423 vtkErrorMacro(<<"Cannot open specified output file...");
424 return;
425 }
426
427 // get image dimensions from the readers first slice
428 this->Reader->GetImageRange(imageRange);
429 tempStructPts = this->Reader->GetImage(imageRange[0]);
430 tempStructPts->GetDimensions(dims);
431 tempStructPts->GetOrigin(origin);
432 tempStructPts->GetSpacing(Spacing);
433
434 dims[2] = (imageRange[1] - imageRange[0] + 1);
435
436 if ( (dims[0]*dims[1]*dims[2]) <= 1 || dims[2] < 2 )
437 {
438 vtkErrorMacro(<<"Bad dimensions...slice must be 3D volume");
439 fclose(outFP);
440 return;
441 }
442
443 xmin[0]=xmin[1]=xmin[2] = VTK_DOUBLE_MAX;
444 xmax[0]=xmax[1]=xmax[2] = -VTK_DOUBLE_MAX;
445
446 inScalars = tempStructPts->GetPointData()->GetScalars();
447 if ( inScalars == NULL )
448 {
449 vtkErrorMacro(<<"Must have scalars to generate isosurface");
450 tempStructPts->Delete();
451 fclose(outFP);
452 return;
453 }
454 inScalars->Register(this);
455 tempStructPts->Delete();
456
457 if (inScalars->GetNumberOfComponents() == 1 )
458 {
459 switch (inScalars->GetDataType())
460 {
461 case VTK_CHAR:
462 {
463 vtkCharArray *scalars = (vtkCharArray *)inScalars;
464 char *s = scalars->GetPointer(0);
465 vtkSliceCubesContour(s,scalars,imageRange,dims,origin,
466 Spacing,this->Value,
467 xmin,xmax,outFP,this->Reader,this->Debug);
468 }
469 break;
470 case VTK_UNSIGNED_CHAR:
471 {
472 vtkUnsignedCharArray *scalars = (vtkUnsignedCharArray *)inScalars;
473 unsigned char *s = scalars->GetPointer(0);
474 vtkSliceCubesContour(s,scalars,imageRange,dims,origin,
475 Spacing,this->Value,
476 xmin,xmax,outFP,this->Reader,this->Debug);
477 }
478 break;
479 case VTK_SHORT:
480 {
481 vtkShortArray *scalars = (vtkShortArray *)inScalars;
482 short *s = scalars->GetPointer(0);
483 vtkSliceCubesContour(s,scalars,imageRange,dims,origin,
484 Spacing,this->Value,
485 xmin,xmax,outFP,this->Reader,this->Debug);
486 }
487 break;
488 case VTK_UNSIGNED_SHORT:
489 {
490 vtkUnsignedShortArray *scalars = (vtkUnsignedShortArray *)inScalars;
491 unsigned short *s = scalars->GetPointer(0);
492 vtkSliceCubesContour(s,scalars,imageRange,dims,origin,
493 Spacing,this->Value,
494 xmin,xmax,outFP,this->Reader,this->Debug);
495 }
496 break;
497 case VTK_INT:
498 {
499 vtkIntArray *scalars = (vtkIntArray *)inScalars;
500 int *s = scalars->GetPointer(0);
501 vtkSliceCubesContour(s,scalars,imageRange,dims,origin,
502 Spacing,this->Value,
503 xmin,xmax,outFP,this->Reader,this->Debug);
504 }
505 break;
506 case VTK_UNSIGNED_INT:
507 {
508 vtkUnsignedIntArray *scalars = (vtkUnsignedIntArray *)inScalars;
509 unsigned int *s = scalars->GetPointer(0);
510 vtkSliceCubesContour(s,scalars,imageRange,dims,origin,
511 Spacing,this->Value,
512 xmin,xmax,outFP,this->Reader,this->Debug);
513 }
514 break;
515 case VTK_LONG:
516 {
517 vtkLongArray *scalars = (vtkLongArray *)inScalars;
518 long *s = scalars->GetPointer(0);
519 vtkSliceCubesContour(s,scalars,imageRange,dims,origin,
520 Spacing,this->Value,
521 xmin,xmax,outFP,this->Reader,this->Debug);
522 }
523 break;
524 case VTK_UNSIGNED_LONG:
525 {
526 vtkUnsignedLongArray *scalars = (vtkUnsignedLongArray *)inScalars;
527 unsigned long *s = scalars->GetPointer(0);
528 vtkSliceCubesContour(s,scalars,imageRange,dims,origin,
529 Spacing,this->Value,
530 xmin,xmax,outFP,this->Reader,this->Debug);
531 }
532 break;
533 case VTK_FLOAT:
534 {
535 vtkFloatArray *scalars = (vtkFloatArray *)inScalars;
536 float *s = scalars->GetPointer(0);
537 vtkSliceCubesContour(s,scalars,imageRange,dims,origin,
538 Spacing,this->Value,
539 xmin,xmax,outFP,this->Reader,this->Debug);
540 }
541 break;
542 case VTK_DOUBLE:
543 {
544 vtkDoubleArray *scalars = (vtkDoubleArray *)inScalars;
545 double *s = scalars->GetPointer(0);
546 vtkSliceCubesContour(s,scalars,imageRange,dims,origin,
547 Spacing,this->Value,
548 xmin,xmax,outFP,this->Reader,this->Debug);
549 }
550 break;
551 }//switch
552 }
553
554 else //multiple components have to convert
555 {
556 vtkDoubleArray *scalars = (vtkDoubleArray *)inScalars;
557 double *s = NULL; //clue to convert data to double
558 vtkSliceCubesContour(s,scalars,imageRange,dims,origin,
559 Spacing,this->Value,
560 xmin,xmax,outFP,this->Reader,this->Debug);
561 }
562
563 inScalars->UnRegister(this);
564
565 if ( this->LimitsFileName )
566 {
567 int i;
568 float t;
569
570 if ( (outFP = fopen(this->LimitsFileName, "wb")) == NULL )
571 {
572 vtkWarningMacro(<<"Sorry, couldn't write limits file...");
573 }
574 else
575 {
576 bool status=true;
577 float forigin[3];
578 for (i=0; i<3 && status; i++)
579 {
580 t = origin[i] + (dims[i] - 1)*Spacing[i];
581 forigin[i] = (float)origin[i];
582 status=vtkByteSwap::SwapWrite4BERange(forigin+i,1,outFP);
583 if(!status)
584 {
585 vtkWarningMacro(<<"SwapWrite failed.");
586 }
587 // swap if necessary
588 if(status)
589 {
590 status=vtkByteSwap::SwapWrite4BERange(&t,1,outFP);
591 if(!status)
592 {
593 vtkWarningMacro(<<"SwapWrite failed.");
594 }
595 }
596 }
597 float ftmp;
598 for (i=0; i<3 && status; i++)
599 {
600 ftmp = (float)xmin[i];
601 status=vtkByteSwap::SwapWrite4BERange(&ftmp,1,outFP);
602 if(!status)
603 {
604 vtkWarningMacro(<<"SwapWrite failed.");
605 }
606 ftmp = (float)xmax[i];
607 if(status)
608 {
609 status=vtkByteSwap::SwapWrite4BERange(&ftmp,1,outFP);
610 if(!status)
611 {
612 vtkWarningMacro(<<"SwapWrite failed.");
613 }
614 }
615 }
616 }
617 fclose(outFP);
618 }
619 }
620
PrintSelf(ostream & os,vtkIndent indent)621 void vtkSliceCubes::PrintSelf(ostream& os, vtkIndent indent)
622 {
623 this->Superclass::PrintSelf(os,indent);
624
625 os << indent << "Iso Value: " << this->Value << "\n";
626
627 if ( this->Reader )
628 {
629 os << indent << "Reader:\n";
630 this->Reader->PrintSelf(os,indent.GetNextIndent());
631 }
632 else
633 {
634 os << indent << "Reader: (none)\n";
635 }
636
637 os << indent << "File Name: "
638 << (this->FileName ? this->FileName : "(none)") << "\n";
639 os << indent << "Limits File Name: "
640 << (this->LimitsFileName ? this->LimitsFileName : "(none)") << "\n";
641 }
642