vtkImageMarchingCubes.cxx 23.9 KB
Newer Older
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20
/*=========================================================================

  Program:   Visualization Toolkit
  Module:    vtkImageMarchingCubes.cxx

  Copyright (c) Ken Martin, Will Schroeder, Bill Lorensen
  All rights reserved.
  See Copyright.txt or http://www.kitware.com/Copyright.htm for details.

     This software is distributed WITHOUT ANY WARRANTY; without even
     the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
     PURPOSE.  See the above copyright notice for more information.

=========================================================================*/
#include "vtkImageMarchingCubes.h"

#include "vtkCellArray.h"
#include "vtkCommand.h"
#include "vtkFloatArray.h"
#include "vtkImageData.h"
21
#include "vtkImageTransform.h"
22
#include "vtkInformation.h"
23
#include "vtkInformationExecutivePortKey.h"
24
#include "vtkInformationVector.h"
25
#include "vtkMarchingCubesTriangleCases.h"
26 27 28 29
#include "vtkObjectFactory.h"
#include "vtkPointData.h"
#include "vtkPolyData.h"
#include "vtkStreamingDemandDrivenPipeline.h"
30
#include <cmath>
31 32 33 34 35 36 37 38 39 40 41 42 43 44

vtkStandardNewMacro(vtkImageMarchingCubes);

//----------------------------------------------------------------------------
// Description:
// Construct object with initial range (0,1) and single contour value
// of 0.0. ComputeNormal is on, ComputeGradients is off and ComputeScalars is on.
vtkImageMarchingCubes::vtkImageMarchingCubes()
{
  this->ContourValues = vtkContourValues::New();
  this->ComputeNormals = 1;
  this->ComputeGradients = 0;
  this->ComputeScalars = 1;

45
  this->LocatorPointIds = nullptr;
46
  this->InputMemoryLimit = 10240;  // 10 mega Bytes
47 48 49 50 51 52 53 54 55 56
}

vtkImageMarchingCubes::~vtkImageMarchingCubes()
{
  this->ContourValues->Delete();
}

// Description:
// Overload standard modified time function. If contour values are modified,
// then this object is modified as well.
57
vtkMTimeType vtkImageMarchingCubes::GetMTime()
58
{
59 60
  vtkMTimeType mTime=this->Superclass::GetMTime();
  vtkMTimeType contourValuesMTime=this->ContourValues->GetMTime();
61

62 63 64 65 66 67 68 69 70 71 72 73
  mTime = ( contourValuesMTime > mTime ? contourValuesMTime : mTime );

  return mTime;
}

//----------------------------------------------------------------------------
template <class T>
int vtkImageMarchingCubesGetTypeSize(T*)
{
  return sizeof(T);
}

74 75 76 77 78 79 80 81 82 83 84 85 86 87 88
//----------------------------------------------------------------------------
int vtkImageMarchingCubes::RequestUpdateExtent(
  vtkInformation* vtkNotUsed(request),
  vtkInformationVector** inputVector,
  vtkInformationVector* vtkNotUsed(outputVector))
{
  // start with an empty UPDATE_EXTENT to ensure proper streaming
  // (the UPDATE_EXTENT will be set properly in RequestData()).
  int extent[6] = { 0, -1, 0, -1, 0, -1 };
  vtkInformation* inInfo = inputVector[0]->GetInformationObject(0);
  inInfo->Set(vtkStreamingDemandDrivenPipeline::UPDATE_EXTENT(),
              extent, 6);
  return 1;
}

89 90 91 92 93 94 95 96 97 98
//----------------------------------------------------------------------------
int vtkImageMarchingCubes::RequestData(
  vtkInformation *vtkNotUsed(request),
  vtkInformationVector **inputVector,
  vtkInformationVector *outputVector)
{
  // get the info objects
  vtkInformation *inInfo = inputVector[0]->GetInformationObject(0);
  vtkInformation *outInfo = outputVector->GetInformationObject(0);

99
  // get the input and output
100 101 102 103 104
  vtkImageData *inData = vtkImageData::SafeDownCast(
    inInfo->Get(vtkDataObject::DATA_OBJECT()));
  vtkPolyData *output = vtkPolyData::SafeDownCast(
    outInfo->Get(vtkDataObject::DATA_OBJECT()));

105 106 107
  vtkDemandDrivenPipeline* inputExec =
    vtkDemandDrivenPipeline::SafeDownCast(
      vtkExecutive::PRODUCER()->GetExecutive(inInfo));
108

109 110
  int numContours=this->ContourValues->GetNumberOfContours();
  double *values=this->ContourValues->GetValues();
111

112
  vtkDebugMacro("Starting Execute Method");
113

114 115
  // Gradients must be computed (but not saved) if Compute normals is on.
  this->NeedGradients = this->ComputeGradients || this->ComputeNormals;
116

117
  // Determine the number of slices per request from input memory limit.
118
  int minSlicesPerChunk, chunkOverlap;
119
  if (this->NeedGradients)
120
  {
121 122
    minSlicesPerChunk = 4;
    chunkOverlap = 3;
123
  }
124
  else
125
  {
126 127
    minSlicesPerChunk = 2;
    chunkOverlap = 1;
128
  }
129
  inputExec->UpdateInformation();
130
  // Each data type requires a different amount of memory.
131
  vtkIdType temp;
132
  switch (inData->GetScalarType())
133
  {
134
    vtkTemplateMacro(
135
      temp = vtkImageMarchingCubesGetTypeSize(static_cast<VTK_TT*>(nullptr))
136 137 138 139
      );
    default:
      vtkErrorMacro(<< "Could not determine input scalar type.");
      return 1;
140
  }
141

142
  // Get the WHOLE_EXTENT and divide it into chunks
143
  int extent[6];
144 145 146 147 148
  inInfo->Get(vtkStreamingDemandDrivenPipeline::WHOLE_EXTENT(), extent);
  // multiply by the area of each slice
  temp *= extent[1] - extent[0] + 1;
  temp *= extent[3] - extent[2] + 1;
  // temp holds memory per image. (+1 to avoid dividing by zero)
149 150
  this->NumberOfSlicesPerChunk =
    static_cast<int>(this->InputMemoryLimit * 1024 / (temp + 1));
151
  if (this->NumberOfSlicesPerChunk < minSlicesPerChunk)
152
  {
153 154
    vtkWarningMacro("Execute: Need "
      <<  minSlicesPerChunk*(temp/1024) << " KB to load "
155 156
      << minSlicesPerChunk << " slices.\n");
    this->NumberOfSlicesPerChunk = minSlicesPerChunk;
157
  }
158
  vtkDebugMacro("Execute: NumberOfSlicesPerChunk = "
159 160
                << this->NumberOfSlicesPerChunk);
  this->NumberOfSlicesPerChunk -= chunkOverlap;
161

162 163
  // Create the points, scalars, normals and Cell arrays for the output.
  // estimate the number of points from the volume dimensions
164 165 166 167 168
  vtkIdType estimatedSize = static_cast<vtkIdType>(extent[1] - extent[0] + 1);
  estimatedSize *= static_cast<vtkIdType>(extent[3] - extent[2] + 1);
  estimatedSize *= static_cast<vtkIdType>(extent[5] - extent[4] + 1);
  estimatedSize = static_cast<vtkIdType>(pow(1.0*estimatedSize, 0.75));
  estimatedSize = (estimatedSize / 1024) * 1024; //multiple of 1024
169
  if (estimatedSize < 1024)
170
  {
171
    estimatedSize = 1024;
172
  }
173 174 175 176 177 178
  vtkDebugMacro(<< "Estimated number of points/triangles: " << estimatedSize);
  this->Points = vtkPoints::New();
  this->Points->Allocate(estimatedSize,estimatedSize/2);
  this->Triangles = vtkCellArray::New();
  this->Triangles->Allocate(estimatedSize,estimatedSize/2);
  if (this->ComputeScalars)
179
  {
180 181
    this->Scalars = vtkFloatArray::New();
    this->Scalars->Allocate(estimatedSize,estimatedSize/2);
182
  }
183
  if (this->ComputeNormals)
184
  {
185 186 187
    this->Normals = vtkFloatArray::New();
    this->Normals->SetNumberOfComponents(3);
    this->Normals->Allocate(3*estimatedSize,3*estimatedSize/2);
188
  }
189
  if (this->ComputeGradients)
190
  {
191 192 193
    this->Gradients = vtkFloatArray::New();
    this->Gradients->SetNumberOfComponents(3);
    this->Gradients->Allocate(3*estimatedSize,3*estimatedSize/2);
194
  }
195

luz.paz's avatar
luz.paz committed
196
  // Initialize the internal point locator (edge table for one image of cubes).
197
  this->InitializeLocator(extent[0], extent[1], extent[2], extent[3]);
198

199
  // Loop through the chunks running marching cubes on each one
200 201 202
  int zMin = extent[4];
  int zMax = extent[5];
  for(int chunkMin = zMin, chunkMax; chunkMin < zMax; chunkMin = chunkMax)
203
  {
204 205 206
    // Get the chunk from the input
    chunkMax = chunkMin + this->NumberOfSlicesPerChunk;
    if (chunkMax > zMax)
207
    {
208
      chunkMax = zMax;
209
    }
210 211 212 213
    extent[4] = chunkMin;
    extent[5] = chunkMax;
    // Expand if computing gradients with central differences
    if (this->NeedGradients)
214
    {
215 216
      --extent[4];
      ++extent[5];
217
    }
218 219
    // Don't go over boundary of data.
    if (extent[4] < zMin)
220
    {
221
      extent[4] = zMin;
222
    }
223
    if (extent[5] > zMax)
224
    {
225
      extent[5] = zMax;
226
    }
227 228 229
    // Get the chunk from the input
    inInfo->Set(vtkStreamingDemandDrivenPipeline::UPDATE_EXTENT(),
                extent, 6);
230
    inputExec->Update();
231

232 233
    this->March(inData, chunkMin, chunkMax, numContours, values);
    if ( !this->AbortExecute )
234
    {
235
      this->UpdateProgress(static_cast<double>(chunkMax-zMin)/(zMax-zMin));
236
    }
237

238 239
    if (vtkDataObject::GetGlobalReleaseDataFlag() ||
        inInfo->Has(vtkStreamingDemandDrivenPipeline::RELEASE_DATA()))
240
    {
241 242
      inData->ReleaseData();
    }
243
  }
244

245
  // Put results in our output
246 247
  vtkDebugMacro(<<"Created: "
               << this->Points->GetNumberOfPoints() << " points, "
248 249 250
               << this->Triangles->GetNumberOfCells() << " triangles");
  output->SetPoints(this->Points);
  this->Points->Delete();
251
  this->Points = nullptr;
252 253
  output->SetPolys(this->Triangles);
  this->Triangles->Delete();
254
  this->Triangles = nullptr;
255
  if (this->ComputeScalars)
256
  {
257 258 259
    int idx = output->GetPointData()->AddArray(this->Scalars);
    output->GetPointData()->SetActiveAttribute(idx, vtkDataSetAttributes::SCALARS);
    this->Scalars->Delete();
260
    this->Scalars = nullptr;
261
  }
262
  if (this->ComputeNormals)
263
  {
264 265
    output->GetPointData()->SetNormals(this->Normals);
    this->Normals->Delete();
266
    this->Normals = nullptr;
267
  }
268

269 270
  // Recover extra space.
  output->Squeeze();
271

272 273
  vtkImageTransform::TransformPointSet(inData,output);

274 275 276 277 278 279 280 281 282 283
  // release the locators memory
  this->DeleteLocator();

  return 1;
}

//----------------------------------------------------------------------------
// This method uses central differences to compute the gradient
// of a point. Note: This method assumes that max > min for all 3 axes!
// It does not consider the dataset spacing.
284
// b0 (b1, b2) indicates the boundary conditions for the three axes.:
285 286 287 288 289 290 291 292 293
// b0 = -1 => pixel is on x axis minimum of region.
// b0 = 0 => no boundary conditions
// b0 = +1 => pixel is on x axis maximum of region.
template <class T>
void vtkImageMarchingCubesComputePointGradient(T *ptr, double *g,
                                               int inc0, int inc1, int inc2,
                                               short b0, short b1, short b2)
{
  if (b0 < 0)
294
  {
295
    g[0] = (double)(ptr[inc0]) - (double)(*ptr);
296
  }
297
  else if (b0 > 0)
298
  {
299
    g[0] = (double)(*ptr) - (double)(ptr[-inc0]);
300
  }
301
  else
302
  {
303
    g[0] = (double)(ptr[inc0]) - (double)(ptr[-inc0]);
304
  }
305 306

  if (b1 < 0)
307
  {
308
    g[1] = (double)(ptr[inc1]) - (double)(*ptr);
309
  }
310
  else if (b1 > 0)
311
  {
312
    g[1] = (double)(*ptr) - (double)(ptr[-inc1]);
313
  }
314
  else
315
  {
316
    g[1] = (double)(ptr[inc1]) - (double)(ptr[-inc1]);
317
  }
318 319

  if (b2 < 0)
320
  {
321
    g[2] = (double)(ptr[inc2]) - (double)(*ptr);
322
  }
323
  else if (b2 > 0)
324
  {
325
    g[2] = (double)(*ptr) - (double)(ptr[-inc2]);
326
  }
327
  else
328
  {
329
    g[2] = (double)(ptr[inc2]) - (double)(ptr[-inc2]);
330
  }
331 332 333 334
}


//----------------------------------------------------------------------------
335
// This method interpolates vertices to make a new point.
336 337 338 339
template <class T>
int vtkImageMarchingCubesMakeNewPoint(vtkImageMarchingCubes *self,
                                      int idx0, int idx1, int idx2,
                                      int inc0, int inc1, int inc2,
340
                                      T *ptr, int edge,
341 342 343 344
                                      int *imageExtent,
                                      double value)
{
  int edgeAxis = 0;
345
  T *ptrB = nullptr;
346 347 348 349
  double temp, pt[3];

  // decode the edge into starting point and axis direction
  switch (edge)
350
  {
351 352 353 354 355 356 357 358 359 360 361 362 363 364 365 366 367 368 369 370 371 372 373 374 375 376 377 378 379 380 381 382 383 384 385 386 387 388 389 390 391 392 393 394 395 396 397 398 399 400 401 402 403 404 405 406 407 408 409 410 411 412 413 414 415 416
    case 0:  // 0,1
      ptrB = ptr + inc0;
      edgeAxis = 0;
      break;
    case 1:  // 1,2
      ++idx0;
      ptr += inc0;
      ptrB = ptr + inc1;
      edgeAxis = 1;
      break;
    case 2:  // 3,2
      ++idx1;
      ptr += inc1;
      ptrB = ptr + inc0;
      edgeAxis = 0;
      break;
    case 3:  // 0,3
      ptrB = ptr + inc1;
      edgeAxis = 1;
      break;
    case 4:  // 4,5
      ++idx2;
      ptr += inc2;
      ptrB = ptr + inc0;
      edgeAxis = 0;
      break;
    case 5:  // 5,6
      ++idx0; ++idx2;
      ptr += inc0 + inc2;
      ptrB = ptr + inc1;
      edgeAxis = 1;
      break;
    case 6:  // 7,6
      ++idx1; ++idx2;
      ptr += inc1 + inc2;
      ptrB = ptr + inc0;
      edgeAxis = 0;
      break;
    case 7: // 4,7
      ++idx2;
      ptr += inc2;
      ptrB = ptr + inc1;
      edgeAxis = 1;
      break;
    case 8: // 0,4
      ptrB = ptr + inc2;
      edgeAxis = 2;
      break;
    case 9: // 1,5
      ++idx0;
      ptr += inc0;
      ptrB = ptr + inc2;
      edgeAxis = 2;
      break;
    case 10: // 3,7
      ++idx1;
      ptr += inc1;
      ptrB = ptr + inc2;
      edgeAxis = 2;
      break;
    case 11: // 2,6
      ++idx0; ++idx1;
      ptr += inc0 + inc1;
      ptrB = ptr + inc2;
      edgeAxis = 2;
      break;
417
  }
418 419 420

  // interpolation factor
  temp = (value - *ptr) / (*ptrB - *ptr);
421

422 423
  // interpolate the point position
  switch (edgeAxis)
424
  {
425
    case 0:
426 427 428
      pt[0] = (double)idx0 + temp;
      pt[1] = (double)idx1;
      pt[2] = (double)idx2;
429 430
      break;
    case 1:
431 432 433
      pt[0] = (double)idx0;
      pt[1] = (double)idx1 + temp;
      pt[2] = (double)idx2;
434 435
      break;
    case 2:
436 437 438
      pt[0] = (double)idx0;
      pt[1] = (double)idx1;
      pt[2] = (double)idx2 + temp;
439
      break;
440
  }
441

442 443
  // Save the scale if we are generating scalars
  if (self->ComputeScalars)
444
  {
445
    self->Scalars->InsertNextValue(value);
446
  }
447

448 449
  // Interpolate to find normal from vectors.
  if (self->NeedGradients)
450
  {
451 452 453 454 455
    short b0, b1, b2;
    double g[3], gB[3];
    // Find boundary conditions and compute gradient (first point)
    b0 = (idx0 == imageExtent[1]);
    if (idx0 == imageExtent[0])
456
    {
457
      b0 = -1;
458
    }
459 460
    b1 = (idx1 == imageExtent[3]);
    if (idx1 == imageExtent[2])
461
    {
462
      b1 = -1;
463
    }
464 465
    b2 = (idx2 == imageExtent[5]);
    if (idx2 == imageExtent[4])
466
    {
467
      b2 = -1;
468
    }
469
    vtkImageMarchingCubesComputePointGradient(ptr, g, inc0, inc1, inc2,
470 471 472
                                             b0, b1, b2);
    // Find boundary conditions and compute gradient (second point)
    switch (edgeAxis)
473
    {
474
      case 0:
475 476 477
        ++idx0;
        b0 = (idx0 == imageExtent[1]);
        break;
478
      case 1:
479 480 481
        ++idx1;
        b1 = (idx1 == imageExtent[3]);
        break;
482
      case 2:
483 484 485
        ++idx2;
        b2 = (idx2 == imageExtent[5]);
        break;
486
    }
487
    vtkImageMarchingCubesComputePointGradient(ptrB, gB, inc0, inc1, inc2,
488 489
                                             b0, b1, b2);
    // Interpolate Gradient
490 491 492
    g[0] = g[0] + temp * (gB[0] - g[0]);
    g[1] = g[1] + temp * (gB[1] - g[1]);
    g[2] = g[2] + temp * (gB[2] - g[2]);
493
    if (self->ComputeGradients)
494
    {
495
      self->Gradients->InsertNextTuple(g);
496
    }
497
    if (self->ComputeNormals)
498
    {
499 500 501 502 503 504
      temp = -1.0 / sqrt(g[0]*g[0] + g[1]*g[1] + g[2]*g[2]);
      g[0] *= temp;
      g[1] *= temp;
      g[2] *= temp;
      self->Normals->InsertNextTuple(g);
    }
505
  }
506

507 508 509 510 511 512 513 514 515 516 517 518 519 520 521 522 523 524 525 526 527 528 529 530
  return self->Points->InsertNextPoint(pt);
}

//----------------------------------------------------------------------------
// This method runs marching cubes on one cube.
template <class T>
void vtkImageMarchingCubesHandleCube(vtkImageMarchingCubes *self,
                                     int cellX, int cellY, int cellZ,
                                     vtkImageData *inData,
                                     T *ptr, int numContours, double *values)
{
  vtkIdType inc0, inc1, inc2;
  int valueIdx;
  double value;
  int cubeIndex, ii;
  vtkIdType pointIds[3];
  vtkMarchingCubesTriangleCases *triCase, *triCases;

  vtkInformation *inInfo = self->GetExecutive()->GetInputInformation(0, 0);

  triCases =  vtkMarchingCubesTriangleCases::GetCases();

  inData->GetIncrements(inc0, inc1, inc2);
  for (valueIdx = 0; valueIdx < numContours; ++valueIdx)
531
  {
532 533 534 535
    value = values[valueIdx];
    // compute the case index
    cubeIndex = 0;
    if ((double)(ptr[0]) > value)
536
    {
537
      cubeIndex += 1;
538
    }
539
    if ((double)(ptr[inc0]) > value)
540
    {
541
      cubeIndex += 2;
542
    }
543
    if ((double)(ptr[inc0 + inc1]) > value)
544
    {
545
      cubeIndex += 4;
546
    }
547
    if ((double)(ptr[inc1]) > value)
548
    {
549
      cubeIndex += 8;
550
    }
551
    if ((double)(ptr[inc2]) > value)
552
    {
553
      cubeIndex += 16;
554
    }
555
    if ((double)(ptr[inc0 + inc2]) > value)
556
    {
557
      cubeIndex += 32;
558
    }
559
    if ((double)(ptr[inc0 + inc1 + inc2]) > value)
560
    {
561
      cubeIndex += 64;
562
    }
563
    if ((double)(ptr[inc1 + inc2]) > value)
564
    {
565
      cubeIndex += 128;
566
    }
567 568
    // Make sure we have trianlges
    if (cubeIndex != 0 && cubeIndex != 255)
569
    {
570 571
      // Get edges.
      triCase = triCases + cubeIndex;
572
      EDGE_LIST *edge = triCase->edges;
573
      // loop over triangles
574
      while(*edge > -1)
575
      {
576
        for (ii=0; ii<3; ++ii, ++edge) //insert triangle
577
        {
578 579 580 581
          // Get the index of the point
          pointIds[ii] = self->GetLocatorPoint(cellX, cellY, *edge);
          // If the point has not been created yet
          if (pointIds[ii] == -1)
582
          {
583 584
            int *extent =
              inInfo->Get(vtkStreamingDemandDrivenPipeline::WHOLE_EXTENT());
585

586 587 588 589
            pointIds[ii] = vtkImageMarchingCubesMakeNewPoint(self,
                                              cellX, cellY, cellZ,
                                              inc0, inc1, inc2,
                                              ptr, *edge, extent,
590
                                              value);
591 592
            self->AddLocatorPoint(cellX, cellY, *edge, pointIds[ii]);
          }
593
        }
594
        self->Triangles->InsertNextCell(3,pointIds);
595
      }//for each triangle
596
    }
597
  }
598 599 600 601 602 603 604 605 606 607 608 609 610 611
}

//----------------------------------------------------------------------------
template <class T>
void vtkImageMarchingCubesMarch(vtkImageMarchingCubes *self,
                                vtkImageData *inData, T *ptr,
                                int chunkMin, int chunkMax,
                                int numContours, double *values)
{
  int idx0, idx1, idx2;
  int min0, max0, min1, max1, min2, max2;
  vtkIdType inc0, inc1, inc2;
  T *ptr0, *ptr1, *ptr2;
  unsigned long target, count;
612

613
  // avoid warnings
614
  (void)ptr;
615

616 617 618 619 620
  // Get information to loop through images.
  inData->GetExtent(min0, max0, min1, max1, min2, max2);
  ptr2 = (T *)(inData->GetScalarPointer(min0, min1, chunkMin));
  inData->GetIncrements(inc0, inc1, inc2);

621
  // Setup the abort interval
622 623 624
  target = (unsigned long)((max0-min0+1) * (max1-min1+1) / 50.0);
  ++target;
  count = 0;
625

626 627
  // Loop over all the cubes
  for (idx2 = chunkMin; idx2 < chunkMax; ++idx2)
628
  {
629 630
    ptr1 = ptr2;
    for (idx1 = min1; idx1 < max1; ++idx1)
631
    {
632
      if (!(count%target))
633
      {
634
        if (self->GetAbortExecute())
635
        {
636 637
          return;
        }
638
      }
639 640 641 642
      count++;
      // continue with last loop
      ptr0 = ptr1;
      for (idx0 = min0; idx0 < max0; ++idx0)
643
      {
644 645 646 647 648 649
        // put magnitudes into the cube structure.
        vtkImageMarchingCubesHandleCube(self, idx0, idx1, idx2, inData, ptr0,
                                       numContours, values);

        ptr0 += inc0;
      }
650 651
      ptr1 += inc1;
    }
652 653
    ptr2 += inc2;
    self->IncrementLocatorZ();
654
  }
655 656 657 658 659 660
}



//----------------------------------------------------------------------------
// This method calls the proper templade function.
661
void vtkImageMarchingCubes::March(vtkImageData *inData,
662 663 664 665
                                 int chunkMin, int chunkMax,
                                 int numContours, double *values)
{
  void *ptr = inData->GetScalarPointer();
666

667
  switch (inData->GetScalarType())
668
  {
669 670 671 672 673 674 675
    vtkTemplateMacro(
      vtkImageMarchingCubesMarch(this, inData, static_cast<VTK_TT*>(ptr),
                                 chunkMin, chunkMax, numContours, values)
      );
    default:
      vtkErrorMacro(<< "Unknown output ScalarType");
      return;
676
  }
677 678 679 680
}


//============================================================================
681
// These method act as the point locator so vertices will be shared.
682 683 684 685 686 687 688 689 690 691
// One 2d array of cubes is stored. (z dimension is ignored).
// Points are indexed by their cube and edge.
// Shared edges are only represented once.  Cubes are responsible for
// edges on their min faces.  Their is an extra row and column of cubes
// to store the max edges of the last row/column of cubes,


//----------------------------------------------------------------------------
// This method allocates and initializes the point array.
// One 2d array of cubes is stored. (z dimension is ignored).
692
void vtkImageMarchingCubes::InitializeLocator(int min0, int max0,
693
                                              int min1, int max1)
694 695
{
  // Free old memory
696 697
  delete [] this->LocatorPointIds;

698 699 700 701 702 703
  // Extra row and column
  this->LocatorDimX = (max0 - min0 + 2);
  this->LocatorDimY = (max1 - min1 + 2);
  this->LocatorMinX = min0;
  this->LocatorMinY = min1;
  // 5 non shared edges.
704 705 706 707
  vtkIdType size = 5;
  size *= static_cast<vtkIdType>(this->LocatorDimX);
  size *= static_cast<vtkIdType>(this->LocatorDimY);
  this->LocatorPointIds = new vtkIdType[size];
708
  // Initialize the array
709
  for (vtkIdType idx = 0; idx < size; ++idx)
710
  {
711
    this->LocatorPointIds[idx] = -1;
712
  }
713 714 715 716 717 718 719
}

//----------------------------------------------------------------------------
// This method frees the locators memory.
void vtkImageMarchingCubes::DeleteLocator()
{
  // Free old memory
720
  delete [] this->LocatorPointIds;
721
  this->LocatorPointIds = nullptr;
722 723 724 725 726 727
}

//----------------------------------------------------------------------------
// This method moves the Z index of the locator up one slice.
void vtkImageMarchingCubes::IncrementLocatorZ()
{
728 729
  vtkIdType *ptr = this->LocatorPointIds;
  for (int y = 0; y < this->LocatorDimY; ++y)
730
  {
731
    for (int x = 0; x < this->LocatorDimX; ++x)
732
    {
733 734 735 736 737
      ptr[0] = ptr[4];
      ptr[3] = ptr[1];
      ptr[1] = ptr[2] = ptr[4] = -1;
      ptr += 5;
    }
738
  }
739 740 741 742 743 744 745 746 747 748 749
}

//----------------------------------------------------------------------------
// This method adds a point to the array.  Cube is the X/Y cube,
// segment is the index of the segment (same as marching cubes).(XYZ)
// (0,0,0)->(1,0,0): 0,  (1,0,0)->(1,1,0): 1,
// (1,1,0)->(0,1,0): 2,  (0,1,0)->(0,0,0): 3,
// (0,0,1)->(1,0,1): 4,  (1,0,1)->(1,1,1): 5,
// (1,1,1)->(0,1,1): 6,  (0,1,1)->(0,0,1): 7,
// (0,0,0)->(0,0,1): 8,  (1,0,0)->(1,0,1): 9,
// (0,1,0)->(0,1,1): 10, (1,1,0)->(1,1,1): 11.
luz.paz's avatar
luz.paz committed
750
// Shared edges are computed internally. (no error checking)
751
void vtkImageMarchingCubes::AddLocatorPoint(int cellX, int cellY, int edge,
752
                                            vtkIdType ptId)
753 754
{
  // Get the correct position in the array.
755
  vtkIdType *ptr = this->GetLocatorPointer(cellX, cellY, edge);
756 757 758 759 760
  *ptr = ptId;
}

//----------------------------------------------------------------------------
// This method gets a point from the locator.
761
vtkIdType vtkImageMarchingCubes::GetLocatorPoint(int cellX, int cellY, int edge)
762 763
{
  // Get the correct position in the array.
764
  vtkIdType *ptr = this->GetLocatorPointer(cellX, cellY, edge);
765 766 767 768 769
  return *ptr;
}

//----------------------------------------------------------------------------
// This method returns a pointer to an ID from a cube and an edge.
770
vtkIdType *vtkImageMarchingCubes::GetLocatorPointer(int cellX,int cellY,int edge)
771 772 773 774
{
  // Remove redundant edges (shared by more than one cube).
  // Take care of shared edges
  switch (edge)
775
  {
776 777 778 779 780 781 782
    case 9:  ++cellX;          edge = 8; break;
    case 10: ++cellY;          edge = 8; break;
    case 11: ++cellX; ++cellY; edge = 8; break;
    case 5:  ++cellX;          edge = 7; break;
    case 6:  ++cellY;          edge = 4; break;
    case 1:  ++cellX;          edge = 3; break;
    case 2:  ++cellY;          edge = 0; break;
783
  }
784

785 786 787 788
  // relative to min and max.
  cellX -= this->LocatorMinX;
  cellY -= this->LocatorMinY;

789
  // compute new indexes for edges (0 to 4)
790
  // must be compatible with LocatorIncrementZ.
791
  if (edge == 7)
792
  {
793
    edge = 1;
794
  }
795
  if (edge == 8)
796
  {
797
    edge = 2;
798
  }
799 800

  // return correct pointer
801
  return this->LocatorPointIds + edge
802
    + (cellX + cellY * static_cast<vtkIdType>(this->LocatorDimX)) * 5;
803 804 805 806 807 808 809 810 811 812 813 814 815 816 817 818 819 820 821 822 823 824
}

//----------------------------------------------------------------------------
int vtkImageMarchingCubes::FillInputPortInformation(int, vtkInformation *info)
{
  info->Set(vtkAlgorithm::INPUT_REQUIRED_DATA_TYPE(), "vtkImageData");
  return 1;
}

//----------------------------------------------------------------------------
void vtkImageMarchingCubes::PrintSelf(ostream& os, vtkIndent indent)
{
  this->Superclass::PrintSelf(os,indent);

  this->ContourValues->PrintSelf(os,indent.GetNextIndent());

  os << indent << "ComputeScalars: " << this->ComputeScalars << "\n";
  os << indent << "ComputeNormals: " << this->ComputeNormals << "\n";
  os << indent << "ComputeGradients: " << this->ComputeGradients << "\n";

  os << indent << "InputMemoryLimit: " << this->InputMemoryLimit <<"K bytes\n";
}