text: /*=========================================================================
  
    Program:   Visualization Toolkit
    Module:    $RCSfile$
  
    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 "vtkDiscreteMarchingCubes.h"
  
  #include "vtkCellArray.h"
  #include "vtkCharArray.h"
  #include "vtkDoubleArray.h"
  #include "vtkFloatArray.h"
  #include "vtkInformation.h"
  #include "vtkInformationVector.h"
  #include "vtkIntArray.h"
  #include "vtkLongArray.h"
  #include "vtkMarchingCubesCases.h"
  #include "vtkMath.h"
  #include "vtkMergePoints.h"
  #include "vtkObjectFactory.h"
  #include "vtkPointData.h"
  #include "vtkPolyData.h"
  #include "vtkCellData.h"
  #include "vtkShortArray.h"
  #include "vtkStructuredPoints.h"
  #include "vtkUnsignedCharArray.h"
  #include "vtkUnsignedIntArray.h"
  #include "vtkUnsignedLongArray.h"
  #include "vtkUnsignedShortArray.h"
  #include "vtkStreamingDemandDrivenPipeline.h"
  
  #ifndef vtkFloatingPointType
  #define vtkFloatingPointType float
  #endif
  
  vtkCxxRevisionMacro(vtkDiscreteMarchingCubes, "$Revision$");
  vtkStandardNewMacro(vtkDiscreteMarchingCubes);
  
  // Description:
  // Construct object with initial range (0,1) and single contour value
  // of 0.0. ComputeNormals is off, ComputeGradients is off and ComputeScalars is on.
  vtkDiscreteMarchingCubes::vtkDiscreteMarchingCubes()
  {
    this->ComputeNormals = 0;
    this->ComputeGradients = 0;
    this->ComputeScalars = 1;
  }
  
  vtkDiscreteMarchingCubes::~vtkDiscreteMarchingCubes()
  {
  }
  
  
  //
  // Contouring filter specialized for volumes and "short int" data values.  
  //
  template <class T>
  void vtkDiscreteMarchingCubesComputeGradient(
    vtkDiscreteMarchingCubes *self,T *scalars, int dims[3], 
    vtkFloatingPointType origin[3], vtkFloatingPointType Spacing[3],
    vtkIncrementalPointLocator *locator, 
    vtkDataArray *newCellScalars, 
    vtkCellArray *newPolys, vtkFloatingPointType *values, 
    int numValues)
  {
    vtkFloatingPointType s[8], value;
    int i, j, k, sliceSize;
    static int CASE_MASK[8] = {1,2,4,8,16,32,64,128};
    vtkMarchingCubesTriangleCases *triCase, *triCases;
    EDGE_LIST  *edge;
    int contNum, jOffset, kOffset, idx, ii, index, *vert;
    vtkIdType ptIds[3];
    int extent[6];
    int ComputeScalars = newCellScalars != NULL;
    vtkFloatingPointType t, *x1, *x2, x[3], min, max;
    vtkFloatingPointType pts[8][3], xp, yp, zp;
    static int edges[12][2] = { {0,1}, {1,2}, {3,2}, {0,3},
                                {4,5}, {5,6}, {7,6}, {4,7},
                                {0,4}, {1,5}, {3,7}, {2,6}};
  
    vtkInformation *inInfo = self->GetExecutive()->GetInputInformation(0, 0);
    inInfo->Get(vtkStreamingDemandDrivenPipeline::WHOLE_EXTENT(),extent);
  
    triCases =  vtkMarchingCubesTriangleCases::GetCases();
  
  //
  // Get min/max contour values
  //
    if ( numValues < 1 )
      {
      return;
      }
    for ( min=max=values[0], i=1; i < numValues; i++)
      {
      if ( values[i] < min )
        {
        min = values[i];
        }
      if ( values[i] > max )
        {
        max = values[i];
        }
      }
  //
  // Traverse all voxel cells, generating triangles
  // using marching cubes algorithm.
  //  
    sliceSize = dims[0] * dims[1];
    for ( k=0; k < (dims[2]-1); k++)
      {
      self->UpdateProgress (static_cast<vtkFloatingPointType>(k)/(dims[2] - 1));
      if (self->GetAbortExecute())
        {
        break;
        }
      kOffset = k*sliceSize;
      pts[0][2] = origin[2] + (k+extent[4])*Spacing[2];
      zp = pts[0][2] + Spacing[2];
      for ( j=0; j < (dims[1]-1); j++)
        {
        jOffset = j*dims[0];
        pts[0][1] = origin[1] + (j+extent[2])*Spacing[1];
        yp = pts[0][1] + Spacing[1];
        for ( i=0; i < (dims[0]-1); i++)
          {
          //get scalar values
          idx = i + jOffset + kOffset;
          s[0] = scalars[idx];
          s[1] = scalars[idx+1];
          s[2] = scalars[idx+1 + dims[0]];
          s[3] = scalars[idx + dims[0]];
          s[4] = scalars[idx + sliceSize];
          s[5] = scalars[idx+1 + sliceSize];
          s[6] = scalars[idx+1 + dims[0] + sliceSize];
          s[7] = scalars[idx + dims[0] + sliceSize];
  
          if ( (s[0] < min && s[1] < min && s[2] < min && s[3] < min &&
                s[4] < min && s[5] < min && s[6] < min && s[7] < min) ||
               (s[0] > max && s[1] > max && s[2] > max && s[3] > max &&
                s[4] > max && s[5] > max && s[6] > max && s[7] > max) )
            {
            continue; // no contours possible
            }
  
          //create voxel points
          pts[0][0] = origin[0] + (i+extent[0])*Spacing[0];
          xp = pts[0][0] + Spacing[0];
  
          pts[1][0] = xp;
          pts[1][1] = pts[0][1];
          pts[1][2] = pts[0][2];
  
          pts[2][0] = xp;
          pts[2][1] = yp;
          pts[2][2] = pts[0][2];
  
          pts[3][0] = pts[0][0];
          pts[3][1] = yp;
          pts[3][2] = pts[0][2];
  
          pts[4][0] = pts[0][0];
          pts[4][1] = pts[0][1];
          pts[4][2] = zp;
  
          pts[5][0] = xp;
          pts[5][1] = pts[0][1];
          pts[5][2] = zp;
  
          pts[6][0] = xp;
          pts[6][1] = yp;
          pts[6][2] = zp;
  
          pts[7][0] = pts[0][0];
          pts[7][1] = yp;
          pts[7][2] = zp;
  
          for (contNum=0; contNum < numValues; contNum++)
            {
            value = values[contNum];
            // Build the case table
            for ( ii=0, index = 0; ii < 8; ii++)
              {
              // for discrete marching cubes, we are looking for an
              // exact match of a scalar at a vertex to a value
              if ( s[ii] == value )
                {
                index |= CASE_MASK[ii];
                }
              }
            if ( index == 0 || index == 255 ) //no surface
              {
              continue;
              }
  
            triCase = triCases+ index;
            edge = triCase->edges;
  
            for ( ; edge[0] > -1; edge += 3 )
              {
              for (ii=0; ii<3; ii++) //insert triangle
                {
                vert = edges[edge[ii]];
                // for discrete marching cubes, the interpolation point
                // is always 0.5.
                t = 0.5;
                x1 = pts[vert[0]];
                x2 = pts[vert[1]];
                x[0] = x1[0] + t * (x2[0] - x1[0]);
                x[1] = x1[1] + t * (x2[1] - x1[1]);
                x[2] = x1[2] + t * (x2[2] - x1[2]);
  
                // add point
                locator->InsertUniquePoint(x, ptIds[ii]);
                }
              // check for degenerate triangle
              if ( ptIds[0] != ptIds[1] &&
                   ptIds[0] != ptIds[2] &&
                   ptIds[1] != ptIds[2] )
                  {
                  newPolys->InsertNextCell(3,ptIds);
                  // Note that DiscreteMarchingCubes stores the scalar
                  // data in the cells. It does not use the point data
                  // since cells from different labeled segments may use
                  // the same point.
                  if (ComputeScalars)
                    {
                    newCellScalars->InsertNextTuple(&value);
                    }
                  }
              }//for each triangle
            }//for all contours
          }//for i
        }//for j
      }//for k
  }
  
  //
  // Contouring filter specialized for volumes and "short int" data values.  
  //
  int vtkDiscreteMarchingCubes::RequestData(
    vtkInformation *vtkNotUsed(request),
    vtkInformationVector **inputVector,
    vtkInformationVector *outputVector)
  {
    vtkInformation *inInfo = inputVector[0]->GetInformationObject(0);
    vtkInformation *outInfo = outputVector->GetInformationObject(0);
  
    vtkPoints *newPts;
    vtkCellArray *newPolys;
    vtkFloatArray *newCellScalars;
    vtkImageData *input = vtkImageData::SafeDownCast(
      inInfo->Get(vtkDataObject::DATA_OBJECT()));
    vtkPointData *pd;
    vtkDataArray *inScalars;
    int dims[3], extent[6];
    int estimatedSize;
    vtkFloatingPointType Spacing[3], origin[3];
    vtkFloatingPointType bounds[6];
    vtkPolyData *output = vtkPolyData::SafeDownCast(
      outInfo->Get(vtkDataObject::DATA_OBJECT()));
    int numContours=this->ContourValues->GetNumberOfContours();
    vtkFloatingPointType *values=this->ContourValues->GetValues();
    
    vtkDebugMacro(<< "Executing marching cubes");
  
  //
  // Initialize and check input
  //
    pd=input->GetPointData();
    if (pd ==NULL)
      {
      vtkErrorMacro(<<"PointData is NULL");
      return 1;
      }
    inScalars=pd->GetScalars();
    if ( inScalars == NULL )
      {
      vtkErrorMacro(<<"Scalars must be defined for contouring");
      return 1;
      }
  
    if ( input->GetDataDimension() != 3 )
      {
      vtkErrorMacro(<<"Cannot contour data of dimension != 3");
      return 1;
      }
    input->GetDimensions(dims);
    input->GetOrigin(origin);
    input->GetSpacing(Spacing);
  
    inInfo->Get(vtkStreamingDemandDrivenPipeline::WHOLE_EXTENT(), extent);
  
    // estimate the number of points from the volume dimensions
    estimatedSize = static_cast<int>(
      pow (static_cast<vtkFloatingPointType>(dims[0] * dims[1] * dims[2]), .75));
    estimatedSize = estimatedSize / 1024 * 1024; //multiple of 1024
    if (estimatedSize < 1024)
      {
      estimatedSize = 1024;
      }
    vtkDebugMacro(<< "Estimated allocation size is " << estimatedSize);
    newPts = vtkPoints::New(); newPts->Allocate(estimatedSize,estimatedSize/2);
    // compute bounds for merging points
    for ( int i=0; i<3; i++)
      {
      bounds[2*i] = origin[i] + extent[2*i] * Spacing[i];
      bounds[2*i+1] = origin[i] + extent[2*i+1] * Spacing[i];
      }
    if ( this->Locator == NULL )
      {
      this->CreateDefaultLocator();
      }
    this->Locator->InitPointInsertion (newPts, bounds, estimatedSize);
  
    newPolys = vtkCellArray::New();
    newPolys->Allocate(newPolys->EstimateSize(estimatedSize,3));
  
    if (this->ComputeScalars)
      {
      newCellScalars = vtkFloatArray::New();
      newCellScalars->Allocate(estimatedSize,3);
      }
    else
      {
      newCellScalars = NULL;
      }
  
    if (inScalars->GetNumberOfComponents() == 1 )
      {
      void* scalars = inScalars->GetVoidPointer(0);
      switch (inScalars->GetDataType())
        {
        vtkTemplateMacro(
          vtkDiscreteMarchingCubesComputeGradient(this,
                                                  static_cast<VTK_TT*>(scalars),
                                                  dims, origin, Spacing,
                                                  this->Locator, newCellScalars,
                                                  newPolys, values, numContours)
          );
        } //switch
      }
  
    else //multiple components - have to convert
      {
      int dataSize = dims[0] * dims[1] * dims[2];
      vtkDoubleArray *image=vtkDoubleArray::New(); 
      image->SetNumberOfComponents(inScalars->GetNumberOfComponents());
      image->SetNumberOfTuples(image->GetNumberOfComponents()*dataSize);
      inScalars->GetTuples(0,dataSize,image);
  
      double *scalars = image->GetPointer(0);
      vtkDiscreteMarchingCubesComputeGradient(this,
                                              scalars,
                                              dims,
                                              origin,
                                              Spacing,
                                              this->Locator,
                                              newCellScalars,
                                              newPolys,
                                              values,
                                              numContours);
      image->Delete();
      }
    
    vtkDebugMacro(<<"Created: " 
                  << newPts->GetNumberOfPoints() << " points, " 
                  << newPolys->GetNumberOfCells() << " triangles");
    //
    // Update ourselves.  Because we don't know up front how many triangles
    // we've created, take care to reclaim memory. 
    //
    output->SetPoints(newPts);
    newPts->Delete();
  
    output->SetPolys(newPolys);
    newPolys->Delete();
  
    if (newCellScalars)
      {
      output->GetCellData()->SetScalars(newCellScalars);
      newCellScalars->Delete();
      }
    output->Squeeze();
    if (this->Locator)
      {
      this->Locator->Initialize(); //free storage
      }
  
    return 1;
  }
  