Filters/Core/vtkMarchingCubes.cxx
//
// Contouring filter specialized for volumes and "short int" data values.
//
Contouring filter specialized for volumes and "short int" data values.
Contouring filter specialized for volumes and " short int " data values .
int vtkMarchingCubes::RequestData(
  vtkInformation *vtkNotUsed(request),
  vtkInformationVector **inputVector,
  vtkInformationVector *outputVector)
{
  vtkInformation *inInfo = inputVector[0]->GetInformationObject(0);
  vtkInformation *outInfo = outputVector->GetInformationObject(0);
  vtkImageData *input = vtkImageData::SafeDownCast(
    inInfo->Get(vtkDataObject::DATA_OBJECT()));
  vtkPolyData *output = vtkPolyData::SafeDownCast(
    outInfo->Get(vtkDataObject::DATA_OBJECT()));
  vtkPoints *newPts;
  vtkCellArray *newPolys;
  vtkFloatArray *newScalars;
  vtkFloatArray *newNormals;
  vtkFloatArray *newGradients;
  vtkPointData *pd;
  vtkDataArray *inScalars;
  int dims[3], extent[6];
  vtkIdType estimatedSize;
  double bounds[6];
  vtkIdType numContours=this->ContourValues->GetNumberOfContours();
  double *values=this->ContourValues->GetValues();
  vtkDebugMacro(<< "Executing marching cubes");
 Initialize and check input
  pd=input->GetPointData();
  if (pd ==nullptr)
  {
    vtkErrorMacro(<<"PointData is nullptr");
    return 1;
  }
  vtkInformationVector *inArrayVec = this->Information->Get(INPUT_ARRAYS_TO_PROCESS());
  if(inArrayVec)
  { we have been passed an input array
    inScalars=this->GetInputArrayToProcess(0,inputVector);
  }
  else
  {
    inScalars = pd->GetScalars();
  }
  if ( inScalars == nullptr )
  {
    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);
  inInfo->Get(vtkStreamingDemandDrivenPipeline::WHOLE_EXTENT(), extent);
   estimate the number of points from the volume dimensions
  estimatedSize = static_cast<vtkIdType>(
    pow(1.0*dims[0]*dims[1]*dims[2], 0.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] = extent[2*i];
    bounds[2*i+1] = extent[2*i+1];
  }
  if ( this->Locator == nullptr )
  {
    this->CreateDefaultLocator();
  }
  this->Locator->InitPointInsertion (newPts, bounds, estimatedSize);
  if (this->ComputeNormals)
  {
    newNormals = vtkFloatArray::New();
    newNormals->SetNumberOfComponents(3);
    newNormals->Allocate(3*estimatedSize,3*estimatedSize/2);
  }
  else
  {
    newNormals = nullptr;
  }
  if (this->ComputeGradients)
  {
    newGradients = vtkFloatArray::New();
    newGradients->SetNumberOfComponents(3);
    newGradients->Allocate(3*estimatedSize,3*estimatedSize/2);
  }
  else
  {
    newGradients = nullptr;
  }
  newPolys = vtkCellArray::New();
  newPolys->Allocate(newPolys->EstimateSize(estimatedSize,3));
  if (this->ComputeScalars)
  {
    newScalars = vtkFloatArray::New();
    newScalars->Allocate(estimatedSize,estimatedSize/2);
  }
  else
  {
    newScalars = nullptr;
  }
  if (inScalars->GetNumberOfComponents() == 1 )
  {
    void* scalars = inScalars->GetVoidPointer(0);
    switch (inScalars->GetDataType())
    {
      vtkTemplateMacro(
        vtkMarchingCubesComputeGradient(this, static_cast<VTK_TT*>(scalars),
                                        dims,this->Locator,
                                        newScalars,newGradients,
                                        newNormals,newPolys,values,
                                        numContours)
        );
    } switch
  }
  else multiple components - have to convert
  {
    vtkIdType dataSize = static_cast<vtkIdType>(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);
    vtkMarchingCubesComputeGradient(this,scalars,dims,this->Locator,
                  newScalars,newGradients,
                  newNormals,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 (newScalars)
  {
    int idx = output->GetPointData()->AddArray(newScalars);
    output->GetPointData()->SetActiveAttribute(idx, vtkDataSetAttributes::SCALARS);
    newScalars->Delete();
  }
  if (newGradients)
  {
    output->GetPointData()->SetVectors(newGradients);
    newGradients->Delete();
  }
  if (newNormals)
  {
    output->GetPointData()->SetNormals(newNormals);
    newNormals->Delete();
  }
  output->Squeeze();
  if (this->Locator)
  {
    this->Locator->Initialize(); free storage
  }
  vtkImageTransform::TransformPointSet(input, output);
  return 1;
}
int vtkMarchingCubes :: RequestData ( vtkInformation * vtkNotUsed ( request ) , vtkInformationVector * * inputVector , vtkInformationVector * outputVector ) { vtkInformation * inInfo = inputVector [ 0 ] -> GetInformationObject ( 0 ) ; vtkInformation * outInfo = outputVector -> GetInformationObject ( 0 ) ; vtkImageData * input = vtkImageData :: SafeDownCast ( inInfo -> Get ( vtkDataObject :: DATA_OBJECT ( ) ) ) ; vtkPolyData * output = vtkPolyData :: SafeDownCast ( outInfo -> Get ( vtkDataObject :: DATA_OBJECT ( ) ) ) ; vtkPoints * newPts ; vtkCellArray * newPolys ; vtkFloatArray * newScalars ; vtkFloatArray * newNormals ; vtkFloatArray * newGradients ; vtkPointData * pd ; vtkDataArray * inScalars ; int dims [ 3 ] , extent [ 6 ] ; vtkIdType estimatedSize ; double bounds [ 6 ] ; vtkIdType numContours = this -> ContourValues -> GetNumberOfContours ( ) ; double * values = this -> ContourValues -> GetValues ( ) ; vtkDebugMacro ( << " " ) ; pd = input -> GetPointData ( ) ; if ( pd == nullptr ) { vtkErrorMacro ( << " " ) ; return 1 ; } vtkInformationVector * inArrayVec = this -> Information -> Get ( INPUT_ARRAYS_TO_PROCESS ( ) ) ; if ( inArrayVec ) { inScalars = this -> GetInputArrayToProcess ( 0 , inputVector ) ; } else { inScalars = pd -> GetScalars ( ) ; } if ( inScalars == nullptr ) { vtkErrorMacro ( << " " ) ; return 1 ; } if ( input -> GetDataDimension ( ) != 3 ) { vtkErrorMacro ( << " " ) ; return 1 ; } input -> GetDimensions ( dims ) ; inInfo -> Get ( vtkStreamingDemandDrivenPipeline :: WHOLE_EXTENT ( ) , extent ) ; estimatedSize = static_cast < vtkIdType > ( pow ( 1.0 * dims [ 0 ] * dims [ 1 ] * dims [ 2 ] , 0.75 ) ) ; estimatedSize = estimatedSize / 1024 * 1024 ; if ( estimatedSize < 1024 ) { estimatedSize = 1024 ; } vtkDebugMacro ( << " " << estimatedSize ) ; newPts = vtkPoints :: New ( ) ; newPts -> Allocate ( estimatedSize , estimatedSize / 2 ) ; for ( int i = 0 ; i < 3 ; i ++ ) { bounds [ 2 * i ] = extent [ 2 * i ] ; bounds [ 2 * i + 1 ] = extent [ 2 * i + 1 ] ; } if ( this -> Locator == nullptr ) { this -> CreateDefaultLocator ( ) ; } this -> Locator -> InitPointInsertion ( newPts , bounds , estimatedSize ) ; if ( this -> ComputeNormals ) { newNormals = vtkFloatArray :: New ( ) ; newNormals -> SetNumberOfComponents ( 3 ) ; newNormals -> Allocate ( 3 * estimatedSize , 3 * estimatedSize / 2 ) ; } else { newNormals = nullptr ; } if ( this -> ComputeGradients ) { newGradients = vtkFloatArray :: New ( ) ; newGradients -> SetNumberOfComponents ( 3 ) ; newGradients -> Allocate ( 3 * estimatedSize , 3 * estimatedSize / 2 ) ; } else { newGradients = nullptr ; } newPolys = vtkCellArray :: New ( ) ; newPolys -> Allocate ( newPolys -> EstimateSize ( estimatedSize , 3 ) ) ; if ( this -> ComputeScalars ) { newScalars = vtkFloatArray :: New ( ) ; newScalars -> Allocate ( estimatedSize , estimatedSize / 2 ) ; } else { newScalars = nullptr ; } if ( inScalars -> GetNumberOfComponents ( ) == 1 ) { void * scalars = inScalars -> GetVoidPointer ( 0 ) ; switch ( inScalars -> GetDataType ( ) ) { vtkTemplateMacro ( vtkMarchingCubesComputeGradient ( this , static_cast < VTK_TT * > ( scalars ) , dims , this -> Locator , newScalars , newGradients , newNormals , newPolys , values , numContours ) ) ; } } else { vtkIdType dataSize = static_cast < vtkIdType > ( 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 ) ; vtkMarchingCubesComputeGradient ( this , scalars , dims , this -> Locator , newScalars , newGradients , newNormals , newPolys , values , numContours ) ; image -> Delete ( ) ; } vtkDebugMacro ( << " " << newPts -> GetNumberOfPoints ( ) << " " << newPolys -> GetNumberOfCells ( ) << " " ) ; output -> SetPoints ( newPts ) ; newPts -> Delete ( ) ; output -> SetPolys ( newPolys ) ; newPolys -> Delete ( ) ; if ( newScalars ) { int idx = output -> GetPointData ( ) -> AddArray ( newScalars ) ; output -> GetPointData ( ) -> SetActiveAttribute ( idx , vtkDataSetAttributes :: SCALARS ) ; newScalars -> Delete ( ) ; } if ( newGradients ) { output -> GetPointData ( ) -> SetVectors ( newGradients ) ; newGradients -> Delete ( ) ; } if ( newNormals ) { output -> GetPointData ( ) -> SetNormals ( newNormals ) ; newNormals -> Delete ( ) ; } output -> Squeeze ( ) ; if ( this -> Locator ) { this -> Locator -> Initialize ( ) ; } vtkImageTransform :: TransformPointSet ( input , output ) ; return 1 ; }
Contouring filter specialized for volumes and "short int" data values.
Contouring filter specialized for volumes and " short int " data values .
// get the info objects // get the input and output // // Initialize and check input // //we have been passed an input array // estimate the number of points from the volume dimensions //multiple of 1024 // compute bounds for merging points //switch //multiple components - have to convert // // Update ourselves.  Because we don't know up front how many triangles // we've created, take care to reclaim memory. // //free storage