You can not select more than 25 topics
Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.
695 lines
21 KiB
695 lines
21 KiB
/*=========================================================================
|
|
|
|
Program: Visualization Toolkit
|
|
Module: $RCSfile: vtkVoxelContoursToSurfaceFilter.cxx,v $
|
|
|
|
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 "vtkVoxelContoursToSurfaceFilter.h"
|
|
|
|
#include "vtkAppendPolyData.h"
|
|
#include "vtkCellArray.h"
|
|
#include "vtkContourFilter.h"
|
|
#include "vtkInformation.h"
|
|
#include "vtkInformationVector.h"
|
|
#include "vtkObjectFactory.h"
|
|
#include "vtkPointData.h"
|
|
#include "vtkPolyData.h"
|
|
#include "vtkStructuredPoints.h"
|
|
|
|
vtkCxxRevisionMacro(vtkVoxelContoursToSurfaceFilter, "$Revision: 1.27 $");
|
|
vtkStandardNewMacro(vtkVoxelContoursToSurfaceFilter);
|
|
|
|
vtkVoxelContoursToSurfaceFilter::vtkVoxelContoursToSurfaceFilter()
|
|
{
|
|
this->MemoryLimitInBytes = 10000000;
|
|
this->Spacing[0] = 1.0;
|
|
this->Spacing[1] = 1.0;
|
|
this->Spacing[2] = 1.0;
|
|
this->LineList = new double[4*1000];
|
|
this->LineListLength = 0;
|
|
this->LineListSize = 1000;
|
|
this->SortedXList = NULL;
|
|
this->SortedYList = NULL;
|
|
this->WorkingList = NULL;
|
|
this->IntersectionList = NULL;
|
|
this->SortedListSize = 0;
|
|
}
|
|
|
|
vtkVoxelContoursToSurfaceFilter::~vtkVoxelContoursToSurfaceFilter()
|
|
{
|
|
if ( this->LineList )
|
|
{
|
|
delete [] this->LineList;
|
|
}
|
|
if ( this->SortedXList )
|
|
{
|
|
delete [] this->SortedXList;
|
|
}
|
|
if ( this->SortedYList )
|
|
{
|
|
delete [] this->SortedYList;
|
|
}
|
|
if ( this->WorkingList )
|
|
{
|
|
delete [] this->WorkingList;
|
|
}
|
|
if ( this->IntersectionList )
|
|
{
|
|
delete [] this->IntersectionList;
|
|
}
|
|
}
|
|
|
|
void vtkVoxelContoursToSurfaceFilter::AddLineToLineList( double x1, double y1,
|
|
double x2, double y2 )
|
|
{
|
|
// Do we need to increase the size of our list?
|
|
if ( this->LineListLength >= this->LineListSize )
|
|
{
|
|
// Double the space we had before
|
|
double *newList = new double[this->LineListSize*4*2];
|
|
memcpy( newList, this->LineList, 4*this->LineListSize*sizeof(double) );
|
|
delete [] this->LineList;
|
|
this->LineList = newList;
|
|
this->LineListSize *= 2;
|
|
}
|
|
|
|
// Now we are sure we have space - add the line
|
|
this->LineList[4*this->LineListLength + 0] = x1;
|
|
this->LineList[4*this->LineListLength + 1] = y1;
|
|
this->LineList[4*this->LineListLength + 2] = x2;
|
|
this->LineList[4*this->LineListLength + 3] = y2;
|
|
this->LineListLength++;
|
|
}
|
|
|
|
void vtkVoxelContoursToSurfaceFilter::SortLineList()
|
|
{
|
|
int i, j;
|
|
double tmp[4];
|
|
double tmpval;
|
|
|
|
|
|
// Make sure we have enough space in our sorted list
|
|
if ( this->SortedListSize < this->LineListLength )
|
|
{
|
|
if ( this->SortedXList )
|
|
{
|
|
delete [] this->SortedXList;
|
|
}
|
|
if ( this->SortedYList )
|
|
{
|
|
delete [] this->SortedYList;
|
|
}
|
|
if ( this->WorkingList )
|
|
{
|
|
delete [] this->WorkingList;
|
|
}
|
|
if ( this->IntersectionList )
|
|
{
|
|
delete [] this->IntersectionList;
|
|
}
|
|
|
|
this->SortedXList = new double[4*this->LineListLength];
|
|
this->SortedYList = new double[4*this->LineListLength];
|
|
this->SortedListSize = this->LineListLength;
|
|
|
|
// Create the space we'll need for our working list of indices
|
|
// The is the list of lines that we are currently considering
|
|
// for intersections. Lines move in then out of the list as
|
|
// we pass the first then the second endpoint. This will be
|
|
// used during the CastXLines and CastYLines methods, and is
|
|
// the same size as the number of lines.
|
|
this->WorkingList = new int[this->LineListLength];
|
|
|
|
// Create the space we'll need for the intersection list
|
|
// There can't be more intersections than there are possible
|
|
// lines. Although it is highly doubtful we'll actually use
|
|
// all this space, it isn't much and it makes the code simpler
|
|
// not to have to worry about exceeding the bounds. This will be
|
|
// used during the CastXLines and CastYLines methods, and is
|
|
// the same size as the number of lines.
|
|
this->IntersectionList = new double[this->LineListLength];
|
|
}
|
|
|
|
// Copy the lines into the lists
|
|
memcpy( this->SortedXList, this->LineList,
|
|
4*this->LineListLength*sizeof(double) );
|
|
memcpy( this->SortedYList, this->LineList,
|
|
4*this->LineListLength*sizeof(double) );
|
|
|
|
// Now sort on x and y
|
|
// Use a simple bubble sort - will improve if necessary
|
|
for ( i = 0; i < this->LineListLength; i++ )
|
|
{
|
|
// swap x entry if necessary to keep min x the first endpoint
|
|
if ( this->SortedXList[4*i + 0] > this->SortedXList[4*i + 2] )
|
|
{
|
|
tmpval = this->SortedXList[4*i];
|
|
this->SortedXList[4*i] = this->SortedXList[4*i + 2];
|
|
this->SortedXList[4*i + 2] = tmpval;
|
|
tmpval = this->SortedXList[4*i + 1];
|
|
this->SortedXList[4*i + 1] = this->SortedXList[4*i + 3];
|
|
this->SortedXList[4*i + 3] = tmpval;
|
|
}
|
|
|
|
// swap y entry if necessary to keep min y the first endpoint
|
|
if ( this->SortedYList[4*i + 1] > this->SortedYList[4*i + 3] )
|
|
{
|
|
tmpval = this->SortedYList[4*i];
|
|
this->SortedYList[4*i] = this->SortedYList[4*i + 2];
|
|
this->SortedYList[4*i + 2] = tmpval;
|
|
tmpval = this->SortedYList[4*i + 1];
|
|
this->SortedYList[4*i + 1] = this->SortedYList[4*i + 3];
|
|
this->SortedYList[4*i + 3] = tmpval;
|
|
}
|
|
|
|
// Sort x list
|
|
for ( j = i; j > 0; j-- )
|
|
{
|
|
if ( this->SortedXList[j*4] < this->SortedXList[(j-1)*4] )
|
|
{
|
|
memcpy( tmp, this->SortedXList + j*4, 4*sizeof(double) );
|
|
memcpy( this->SortedXList + j*4,
|
|
this->SortedXList + (j-1)*4, 4*sizeof(double) );
|
|
memcpy( this->SortedXList + (j-1)*4, tmp, 4*sizeof(double) );
|
|
}
|
|
else
|
|
{
|
|
break;
|
|
}
|
|
}
|
|
|
|
// Sort y list
|
|
for ( j = i; j > 0; j-- )
|
|
{
|
|
if ( this->SortedYList[j*4+1] < this->SortedYList[(j-1)*4+1] )
|
|
{
|
|
memcpy( tmp, this->SortedYList + j*4, 4*sizeof(double) );
|
|
memcpy( this->SortedYList + j*4, this->SortedYList + (j-1)*4,
|
|
4*sizeof(double) );
|
|
memcpy( this->SortedYList + (j-1)*4, tmp, 4*sizeof(double) );
|
|
}
|
|
else
|
|
{
|
|
break;
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
|
|
void vtkVoxelContoursToSurfaceFilter::CastLines( float *slicePtr,
|
|
double gridOrigin[3],
|
|
int gridSize[3],
|
|
int type )
|
|
{
|
|
double axis1, axis2;
|
|
double d1, d2;
|
|
int index;
|
|
int i, j;
|
|
double tmp;
|
|
double *line;
|
|
float *currSlicePtr;
|
|
int currSlice;
|
|
int currentIntersection;
|
|
double sign;
|
|
double *sortedList;
|
|
double low1, low2, high1, high2;
|
|
int increment1, increment2;
|
|
int offset1, offset2, offset3, offset4;
|
|
|
|
// this is the x direction
|
|
if ( type == 0 )
|
|
{
|
|
low1 = gridOrigin[0];
|
|
high1 = gridOrigin[0] + (double)gridSize[0];
|
|
low2 = gridOrigin[1];
|
|
high2 = gridOrigin[1] + (double)gridSize[1];
|
|
increment1 = gridSize[0];
|
|
increment2 = 1;
|
|
sortedList = this->SortedXList;
|
|
offset1 = 0;
|
|
offset2 = 2;
|
|
offset3 = 1;
|
|
offset4 = 3;
|
|
}
|
|
// This is the y direction
|
|
else
|
|
{
|
|
low1 = gridOrigin[1];
|
|
high1 = gridOrigin[1] + (double)gridSize[1];
|
|
low2 = gridOrigin[0];
|
|
high2 = gridOrigin[0] + (double)gridSize[0];
|
|
increment1 = 1;
|
|
increment2 = gridSize[0];
|
|
sortedList = this->SortedYList;
|
|
offset1 = 1;
|
|
offset2 = 3;
|
|
offset3 = 0;
|
|
offset4 = 2;
|
|
}
|
|
|
|
// Initialize the working list to nothing. We will start
|
|
// looking at index = 0 for the next line to add to the
|
|
// working list
|
|
this->WorkingListLength = 0;
|
|
index = 0;
|
|
|
|
// Loop through the x or y lines
|
|
for ( axis1 = low1, currSlice = 0; axis1 < high1; axis1 += 1.0, currSlice++ )
|
|
{
|
|
// Initialize the intersection list to nothing
|
|
this->IntersectionListLength = 0;
|
|
|
|
// Add lines to the working list if necessary
|
|
for ( ; index < this->LineListLength; index++ )
|
|
{
|
|
if ( sortedList[4*index + offset1] < axis1 )
|
|
{
|
|
this->WorkingList[this->WorkingListLength] = index;
|
|
this->WorkingListLength++;
|
|
}
|
|
else
|
|
{
|
|
break;
|
|
}
|
|
}
|
|
|
|
// Do the intersections, removing lines from the
|
|
// working list if necessary
|
|
for ( i = 0; i < this->WorkingListLength; i++ )
|
|
{
|
|
line = sortedList + 4*this->WorkingList[i];
|
|
|
|
// Yes, it intersects, add it to the intersection list
|
|
if ( line[offset1] < axis1 && line[offset2] > axis1 )
|
|
{
|
|
// Compute the intersection distance
|
|
// For x lines this is y = y1 + (y2 - y1)*((x - x1)/(x2 - x1))
|
|
// For y lines this is x = x1 + (x2 - x1)*((y - y1)/(y2 - y1))
|
|
this->IntersectionList[this->IntersectionListLength] =
|
|
line[offset3] + (line[offset4] - line[offset3]) *
|
|
((axis1 - line[offset1]) / (line[offset2] - line[offset1] ));
|
|
|
|
// Make sure this distance is sorted
|
|
for ( j = this->IntersectionListLength; j > 0; j-- )
|
|
{
|
|
if ( this->IntersectionList[j] < this->IntersectionList[j-1] )
|
|
{
|
|
tmp = this->IntersectionList[j];
|
|
this->IntersectionList[j] = this->IntersectionList[j-1];
|
|
this->IntersectionList[j-1] = tmp;
|
|
}
|
|
else
|
|
{
|
|
break;
|
|
}
|
|
}
|
|
this->IntersectionListLength++;
|
|
}
|
|
// No, it doesn't intersect, remove it from the working list
|
|
else
|
|
{
|
|
for ( j = i; j < (this->WorkingListLength-1); j++ )
|
|
{
|
|
this->WorkingList[j] = this->WorkingList[j+1];
|
|
}
|
|
this->WorkingListLength--;
|
|
i--;
|
|
}
|
|
}
|
|
|
|
// Now we have all the intersections for the x or y line, in sorted
|
|
// order. Use them to fill in distances (as long as there are
|
|
// any)
|
|
if ( this->IntersectionListLength )
|
|
{
|
|
currSlicePtr = slicePtr + currSlice*increment2;
|
|
currentIntersection = 0;
|
|
// We are starting outside which has a negative distance
|
|
sign = -1.0;
|
|
for ( axis2 = low2; axis2 < high2; axis2 += 1.0 )
|
|
{
|
|
while( currentIntersection < this->IntersectionListLength &&
|
|
this->IntersectionList[currentIntersection] < axis2 )
|
|
|
|
{
|
|
currentIntersection++;
|
|
|
|
// Each time we cross a line we are moving across an
|
|
// inside/outside boundary
|
|
sign *= -1.0;
|
|
}
|
|
// We are now positioned at an x or y value between currentIntersection
|
|
// and currentIntersection - 1 (except at boundaries where we are
|
|
// before intersection 0 or after the last intersection)
|
|
|
|
if ( currentIntersection == 0 )
|
|
{
|
|
d1 = axis2 - this->IntersectionList[currentIntersection];
|
|
*currSlicePtr = (*currSlicePtr > d1 )?(*currSlicePtr):(d1);
|
|
}
|
|
else if ( currentIntersection == this->IntersectionListLength )
|
|
{
|
|
d1 = this->IntersectionList[currentIntersection-1] - axis2;
|
|
*currSlicePtr = (*currSlicePtr > d1 )?(*currSlicePtr):(d1);
|
|
}
|
|
else
|
|
{
|
|
d1 = axis2 - this->IntersectionList[currentIntersection-1];
|
|
d2 = this->IntersectionList[currentIntersection] - axis2;
|
|
d1 = ( d1 < d2 )?(d1):(d2);
|
|
if ( type == 0 )
|
|
{
|
|
*currSlicePtr = sign*d1;
|
|
}
|
|
else
|
|
{
|
|
*currSlicePtr =
|
|
(sign*(*currSlicePtr) < d1 )?(*currSlicePtr):(sign*d1);
|
|
}
|
|
}
|
|
|
|
currSlicePtr += increment1;
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
void vtkVoxelContoursToSurfaceFilter::PushDistances( float *volumePtr,
|
|
int gridSize[3],
|
|
int chunkSize )
|
|
{
|
|
int i, j, k;
|
|
float *vptr;
|
|
|
|
// Push distances along x (both ways) and y (both ways) on each slice
|
|
for ( k = 0; k < chunkSize; k++ )
|
|
{
|
|
// Do the x rows
|
|
for ( j = 0; j < gridSize[1]; j++ )
|
|
{
|
|
vptr = volumePtr + k*gridSize[0]*gridSize[1] + j*gridSize[0];
|
|
vptr++;
|
|
|
|
// first one way
|
|
for ( i = 1; i < gridSize[0]; i++ )
|
|
{
|
|
if ( *vptr > 0 && *(vptr-1) + 1 < *(vptr) )
|
|
{
|
|
*vptr = *(vptr-1) + 1;
|
|
}
|
|
else if ( *vptr < 0 && *(vptr-1) - 1 > *(vptr) )
|
|
{
|
|
*vptr = *(vptr-1) - 1;
|
|
}
|
|
vptr++;
|
|
}
|
|
|
|
vptr -= 2;
|
|
i -= 2;
|
|
|
|
// then the other
|
|
for ( ; i >= 0; i-- )
|
|
{
|
|
if ( *vptr > 0 && *(vptr+1) + 1 < *(vptr) )
|
|
{
|
|
*vptr = *(vptr+1) + 1;
|
|
}
|
|
else if ( *vptr < 0 && *(vptr+1) - 1 > *(vptr) )
|
|
{
|
|
*vptr = *(vptr+1) - 1;
|
|
}
|
|
}
|
|
|
|
}
|
|
|
|
|
|
// Do the y columns
|
|
for ( i = 0; i < gridSize[0]; i++ )
|
|
{
|
|
vptr = volumePtr + k*gridSize[0]*gridSize[1] + i;
|
|
|
|
vptr+=gridSize[0];
|
|
|
|
// first one way
|
|
for ( j = 1; j < gridSize[1]; j++ )
|
|
{
|
|
if ( *vptr > 0 && *(vptr-gridSize[0]) + 1 < *(vptr) )
|
|
{
|
|
*vptr = *(vptr-gridSize[0]) + 1;
|
|
}
|
|
else if ( *vptr < 0 && *(vptr-gridSize[0]) - 1 > *(vptr) )
|
|
{
|
|
*vptr = *(vptr-gridSize[0]) - 1;
|
|
}
|
|
vptr += gridSize[0];
|
|
}
|
|
|
|
vptr -= 2*gridSize[0];
|
|
j -= 2;
|
|
|
|
// then the other
|
|
for ( ; j >= 0; j-- )
|
|
{
|
|
if ( *vptr > 0 && *(vptr+gridSize[0]) + 1 < *(vptr) )
|
|
{
|
|
*vptr = *(vptr+gridSize[0]) + 1;
|
|
}
|
|
else if ( *vptr < 0 && *(vptr+gridSize[0]) - 1 > *(vptr) )
|
|
{
|
|
*vptr = *(vptr+gridSize[0]) - 1;
|
|
}
|
|
}
|
|
|
|
}
|
|
}
|
|
}
|
|
|
|
// Append data sets into single unstructured grid
|
|
int vtkVoxelContoursToSurfaceFilter::RequestData(
|
|
vtkInformation *vtkNotUsed(request),
|
|
vtkInformationVector **inputVector,
|
|
vtkInformationVector *outputVector)
|
|
{
|
|
// get the info objects
|
|
vtkInformation *inInfo = inputVector[0]->GetInformationObject(0);
|
|
vtkInformation *outInfo = outputVector->GetInformationObject(0);
|
|
|
|
// get the input and ouptut
|
|
vtkPolyData *input = vtkPolyData::SafeDownCast(
|
|
inInfo->Get(vtkDataObject::DATA_OBJECT()));
|
|
vtkPolyData *output = vtkPolyData::SafeDownCast(
|
|
outInfo->Get(vtkDataObject::DATA_OBJECT()));
|
|
|
|
vtkCellArray *inputPolys = input->GetPolys();
|
|
int gridSize[3];
|
|
double gridOrigin[3];
|
|
double contourBounds[6];
|
|
int chunkSize;
|
|
int currentSlice, lastSlice, currentIndex;
|
|
int i, j;
|
|
int numberOfInputCells;
|
|
int currentInputCellIndex;
|
|
vtkIdType npts = 0;
|
|
vtkIdType *pts = 0;
|
|
double point1[3], point2[3];
|
|
double currentZ;
|
|
vtkStructuredPoints *volume;
|
|
float *volumePtr, *slicePtr;
|
|
vtkContourFilter *contourFilter;
|
|
vtkPolyData *contourOutput;
|
|
vtkAppendPolyData *appendFilter;
|
|
|
|
vtkDebugMacro(<<"Creating surfaces from contours");
|
|
|
|
// Get the bounds of the input contours
|
|
input->GetBounds( contourBounds );
|
|
|
|
if (contourBounds[0] > contourBounds[1])
|
|
{ // empty input
|
|
return 1;
|
|
}
|
|
|
|
// From the bounds, compute the grid size, and origin
|
|
|
|
// The origin of the grid should be (-0.5, -0.5, 0.0) away from the
|
|
// lower bounds of the contours. This is because we want the grid
|
|
// to lie halfway between integer endpoint locations of the line
|
|
// segments on each plane. Also, we want an extra plane on each end
|
|
// for capping
|
|
gridOrigin[0] = contourBounds[0] - 0.5;
|
|
gridOrigin[1] = contourBounds[2] - 0.5;
|
|
gridOrigin[2] = contourBounds[4] - 1.0;
|
|
|
|
// The difference between the bounds, plus one to account a
|
|
// sample on the first and last location, plus one to account
|
|
// for the larger grid size ( the 0.5 unit border ) On Z, we
|
|
// want to sample exactly on the contours so we don't need to
|
|
// add the extra 1, but we have added two extra planes so we
|
|
// need another 2.
|
|
gridSize[0] = (int) (contourBounds[1] - contourBounds[0] + 2);
|
|
gridSize[1] = (int) (contourBounds[3] - contourBounds[2] + 2);
|
|
gridSize[2] = (int) (contourBounds[5] - contourBounds[4] + 3);
|
|
|
|
// How many slices in a chunk? This will later be decremented
|
|
// by one to account for the fact that the last slice in the
|
|
// previous chuck is copied to the first slice in the next chunk.
|
|
// Stay within memory limit. There are 4 bytes per double.
|
|
chunkSize = this->MemoryLimitInBytes / ( gridSize[0] * gridSize[1] * 4 );
|
|
if ( chunkSize > gridSize[2] )
|
|
{
|
|
chunkSize = gridSize[2];
|
|
}
|
|
|
|
currentSlice = 0;
|
|
currentZ = contourBounds[4] - 1.0;
|
|
currentIndex = 0;
|
|
lastSlice = gridSize[2] - 1;
|
|
numberOfInputCells = inputPolys->GetNumberOfCells();
|
|
currentInputCellIndex = 0;
|
|
|
|
volume = vtkStructuredPoints::New();
|
|
volume->SetDimensions( gridSize[0], gridSize[1], chunkSize );
|
|
volume->SetSpacing( this->Spacing );
|
|
volume->SetScalarType( VTK_FLOAT );
|
|
volume->AllocateScalars();
|
|
volumePtr =
|
|
(float *)(volume->GetPointData()->GetScalars()->GetVoidPointer(0));
|
|
|
|
|
|
contourFilter = vtkContourFilter::New();
|
|
contourFilter->SetInput( volume );
|
|
contourFilter->SetNumberOfContours(1);
|
|
contourFilter->SetValue( 0, 0.0 );
|
|
|
|
appendFilter = vtkAppendPolyData::New();
|
|
|
|
inputPolys->InitTraversal();
|
|
inputPolys->GetNextCell( npts, pts );
|
|
|
|
while ( currentSlice <= lastSlice )
|
|
{
|
|
// Make sure the origin of the volume is in the right
|
|
// place so that the appended polydata all matches up
|
|
// nicely.
|
|
volume->SetOrigin( gridOrigin[0], gridOrigin[1],
|
|
gridOrigin[2] +
|
|
this->Spacing[2] * (currentSlice - (currentSlice!=0)) );
|
|
|
|
for ( i = currentIndex; i < chunkSize; i++ )
|
|
{
|
|
slicePtr = volumePtr + i * gridSize[0] * gridSize[1];
|
|
|
|
// Clear out the slice memory - set it all to a large negative
|
|
// value indicating no surfaces are nearby, and we assume we
|
|
// are outside of any surface
|
|
for ( j = 0; j < gridSize[0] * gridSize[1]; j++ )
|
|
{
|
|
*(slicePtr+j) = -9.99e10;
|
|
}
|
|
|
|
// If we are past the end, don't do anything
|
|
if ( currentSlice > lastSlice )
|
|
{
|
|
continue;
|
|
}
|
|
|
|
this->LineListLength = 0;
|
|
|
|
// Read in the lines for the contours on this slice
|
|
while ( currentInputCellIndex < numberOfInputCells )
|
|
{
|
|
// Check if we are still on the right z slice
|
|
input->GetPoint( pts[0], point1 );
|
|
if ( point1[2] != currentZ )
|
|
{
|
|
break;
|
|
}
|
|
|
|
// This contour is on the right z slice - add the lines
|
|
// to our list
|
|
for ( j = 0; j < npts; j++ )
|
|
{
|
|
input->GetPoint( pts[j], point1 );
|
|
input->GetPoint( pts[(j+1)%npts], point2 );
|
|
this->AddLineToLineList( point1[0], point1[1],
|
|
point2[0], point2[1] );
|
|
}
|
|
|
|
inputPolys->GetNextCell( npts, pts );
|
|
currentInputCellIndex++;
|
|
}
|
|
|
|
// Sort the contours in x and y
|
|
this->SortLineList();
|
|
|
|
// Cast lines in x and y filling in distance
|
|
this->CastLines( slicePtr, gridOrigin, gridSize, 0 );
|
|
this->CastLines( slicePtr, gridOrigin, gridSize, 1 );
|
|
|
|
// Move on to the next slice
|
|
currentSlice++;
|
|
currentIndex++;
|
|
currentZ += 1.0;
|
|
}
|
|
|
|
this->PushDistances( volumePtr, gridSize, chunkSize );
|
|
|
|
// Update the contour filter and grab the output
|
|
// Make a new output for it, then grab the output and
|
|
// add it to the append filter, then delete the output
|
|
// which is ok since it was registered by the appendFilter
|
|
contourOutput = vtkPolyData::New();
|
|
contourFilter->Update();
|
|
contourOutput->ShallowCopy(contourFilter->GetOutput());
|
|
appendFilter->AddInput( contourOutput );
|
|
contourOutput->Delete();
|
|
|
|
|
|
if ( currentSlice <= lastSlice )
|
|
{
|
|
// Copy last slice to first slice
|
|
memcpy( volumePtr, volumePtr + (chunkSize-1)*gridSize[0]*gridSize[1],
|
|
sizeof(float) * gridSize[0] * gridSize[1] );
|
|
|
|
// reset currentIndex to 1
|
|
currentIndex = 1;
|
|
}
|
|
}
|
|
|
|
appendFilter->Update();
|
|
|
|
// Grab the appended data as the output to this filter
|
|
output->SetPoints( appendFilter->GetOutput()->GetPoints() );
|
|
output->SetVerts( appendFilter->GetOutput()->GetVerts() );
|
|
output->SetLines( appendFilter->GetOutput()->GetLines() );
|
|
output->SetPolys( appendFilter->GetOutput()->GetPolys() );
|
|
output->SetStrips( appendFilter->GetOutput()->GetStrips() );
|
|
output->GetPointData()->PassData(appendFilter->GetOutput()->GetPointData());
|
|
|
|
contourFilter->Delete();
|
|
appendFilter->Delete();
|
|
volume->Delete();
|
|
|
|
return 1;
|
|
}
|
|
|
|
|
|
void vtkVoxelContoursToSurfaceFilter::PrintSelf(ostream& os, vtkIndent indent)
|
|
{
|
|
this->Superclass::PrintSelf(os,indent);
|
|
|
|
os << indent << "Memory Limit (in bytes): " <<
|
|
this->MemoryLimitInBytes << endl;
|
|
|
|
os << indent << "Spacing: " << this->Spacing[0] << " " <<
|
|
this->Spacing[1] << " " << this->Spacing[2] << endl;
|
|
}
|
|
|