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