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.
536 lines
15 KiB
536 lines
15 KiB
/*=========================================================================
|
|
|
|
Program: Visualization Toolkit
|
|
Module: $RCSfile: vtkSynchronizedTemplates2D.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 "vtkSynchronizedTemplates2D.h"
|
|
|
|
#include "vtkCellArray.h"
|
|
#include "vtkCellData.h"
|
|
#include "vtkCharArray.h"
|
|
#include "vtkDoubleArray.h"
|
|
#include "vtkFloatArray.h"
|
|
#include "vtkImageData.h"
|
|
#include "vtkInformation.h"
|
|
#include "vtkInformationVector.h"
|
|
#include "vtkIntArray.h"
|
|
#include "vtkLongArray.h"
|
|
#include "vtkObjectFactory.h"
|
|
#include "vtkPointData.h"
|
|
#include "vtkPolyData.h"
|
|
#include "vtkShortArray.h"
|
|
#include "vtkStreamingDemandDrivenPipeline.h"
|
|
#include "vtkUnsignedCharArray.h"
|
|
#include "vtkUnsignedIntArray.h"
|
|
#include "vtkUnsignedLongArray.h"
|
|
#include "vtkUnsignedShortArray.h"
|
|
|
|
#include <math.h>
|
|
|
|
vtkCxxRevisionMacro(vtkSynchronizedTemplates2D, "$Revision: 1.7 $");
|
|
vtkStandardNewMacro(vtkSynchronizedTemplates2D);
|
|
|
|
//----------------------------------------------------------------------------
|
|
// Description:
|
|
// Construct object with initial scalar range (0,1) and single contour value
|
|
// of 0.0. The ImageRange are set to extract the first k-plane.
|
|
vtkSynchronizedTemplates2D::vtkSynchronizedTemplates2D()
|
|
{
|
|
this->ContourValues = vtkContourValues::New();
|
|
this->ComputeScalars = 1;
|
|
this->ArrayComponent = 0;
|
|
|
|
// by default process active point scalars
|
|
this->SetInputArrayToProcess(0,0,0,vtkDataObject::FIELD_ASSOCIATION_POINTS,
|
|
vtkDataSetAttributes::SCALARS);
|
|
}
|
|
|
|
vtkSynchronizedTemplates2D::~vtkSynchronizedTemplates2D()
|
|
{
|
|
this->ContourValues->Delete();
|
|
}
|
|
|
|
//----------------------------------------------------------------------------
|
|
// Description:
|
|
// Overload standard modified time function. If contour values are modified,
|
|
// then this object is modified as well.
|
|
unsigned long vtkSynchronizedTemplates2D::GetMTime()
|
|
{
|
|
unsigned long mTime=this->Superclass::GetMTime();
|
|
unsigned long mTime2=this->ContourValues->GetMTime();
|
|
|
|
mTime = ( mTime2 > mTime ? mTime2 : mTime );
|
|
return mTime;
|
|
}
|
|
|
|
//----------------------------------------------------------------------------
|
|
//
|
|
// Contouring filter specialized for images
|
|
//
|
|
template <class T>
|
|
void vtkContourImage(vtkSynchronizedTemplates2D *self,
|
|
T *scalars, vtkPoints *newPts,
|
|
vtkDataArray *newScalars, vtkCellArray *lines,
|
|
vtkImageData *input, int *updateExt)
|
|
{
|
|
double *values = self->GetValues();
|
|
int numContours = self->GetNumberOfContours();
|
|
T *inPtr, *rowPtr;
|
|
double x[3];
|
|
double *origin = input->GetOrigin();
|
|
double *spacing = input->GetSpacing();
|
|
double y, t;
|
|
int *isect1Ptr, *isect2Ptr;
|
|
vtkIdType ptIds[2];
|
|
int *tablePtr;
|
|
int v0, v1 = 0, v2;
|
|
int idx, vidx;
|
|
double s0, s1, s2, value;
|
|
int i, j;
|
|
int lineCases[64];
|
|
|
|
// The update extent may be different than the extent of the image.
|
|
// The only problem with using the update extent is that one or two
|
|
// sources enlarge the update extent. This behavior is slated to be
|
|
// eliminated.
|
|
vtkIdType *incs = input->GetIncrements();
|
|
int *ext = input->GetExtent();
|
|
int axis0, axis1;
|
|
int min0, max0, dim0;
|
|
int min1, max1;
|
|
int inc0, inc1;
|
|
|
|
// Figure out which plane the image lies in.
|
|
if (updateExt[4] == updateExt[5])
|
|
{ // z collapsed
|
|
axis0 = 0;
|
|
min0 = updateExt[0];
|
|
max0 = updateExt[1];
|
|
inc0 = incs[0];
|
|
axis1 = 1;
|
|
min1 = updateExt[2];
|
|
max1 = updateExt[3];
|
|
inc1 = incs[1];
|
|
x[2] = origin[2] + (updateExt[4]*spacing[2]);
|
|
}
|
|
else if (updateExt[2] == updateExt[3])
|
|
{ // y collapsed
|
|
axis0 = 0;
|
|
min0 = updateExt[0];
|
|
max0 = updateExt[1];
|
|
inc0 = incs[0];
|
|
axis1 = 2;
|
|
min1 = updateExt[4];
|
|
max1 = updateExt[5];
|
|
inc1 = incs[2];
|
|
x[1] = origin[1] + (updateExt[2]*spacing[1]);
|
|
}
|
|
else if (updateExt[0] == updateExt[1])
|
|
{ // x collapsed
|
|
axis0 = 1;
|
|
min0 = updateExt[2];
|
|
max0 = updateExt[3];
|
|
inc0 = incs[1];
|
|
axis1 = 2;
|
|
min1 = updateExt[4];
|
|
max1 = updateExt[5];
|
|
inc1 = incs[2];
|
|
x[0] = origin[0] + (updateExt[0]*spacing[0]);
|
|
}
|
|
else
|
|
{
|
|
vtkGenericWarningMacro("Expecting 2D data.");
|
|
return;
|
|
}
|
|
dim0 = max0-min0+1;
|
|
|
|
|
|
// setup the table entries
|
|
for (i = 0; i < 64; i++)
|
|
{
|
|
lineCases[i] = -1;
|
|
}
|
|
|
|
lineCases[12] = 3;
|
|
lineCases[13] = dim0*2;
|
|
|
|
lineCases[20] = 1;
|
|
lineCases[21] = dim0*2;
|
|
|
|
lineCases[24] = 1;
|
|
lineCases[25] = 3;
|
|
|
|
lineCases[36] = 0;
|
|
lineCases[37] = dim0*2;
|
|
|
|
lineCases[40] = 0;
|
|
lineCases[41] = 3;
|
|
|
|
lineCases[48] = 0;
|
|
lineCases[49] = 1;
|
|
|
|
lineCases[60] = 0;
|
|
lineCases[61] = 1;
|
|
lineCases[62] = 3;
|
|
lineCases[63] = dim0*2;
|
|
|
|
// allocate storage arrays
|
|
int *isect1 = new int [dim0*4];
|
|
isect1[dim0*2-2] = -1;
|
|
isect1[dim0*2-1] = -1;
|
|
isect1[dim0*4-2] = -1;
|
|
isect1[dim0*4-1] = -1;
|
|
|
|
|
|
// Compute the staring location. We may be operating
|
|
// on a part of the image.
|
|
scalars += incs[0]*(updateExt[0]-ext[0])
|
|
+ incs[1]*(updateExt[2]-ext[2])
|
|
+ incs[2]*(updateExt[4]-ext[4]) + self->GetArrayComponent();
|
|
|
|
// for each contour
|
|
for (vidx = 0; vidx < numContours; vidx++)
|
|
{
|
|
rowPtr = scalars;
|
|
inPtr = rowPtr;
|
|
|
|
lineCases[13] = dim0*2;
|
|
lineCases[21] = dim0*2;
|
|
lineCases[37] = dim0*2;
|
|
lineCases[63] = dim0*2;
|
|
|
|
value = values[vidx];
|
|
|
|
// Traverse all pixel cells, generating line segements using templates
|
|
for (j = min1; j <= max1; j++)
|
|
{
|
|
inPtr = rowPtr;
|
|
rowPtr += inc1;
|
|
|
|
// set the y coordinate
|
|
y = origin[axis1] + j*spacing[axis1];
|
|
// first compute the intersections
|
|
s1 = *inPtr;
|
|
|
|
// swap the buffers
|
|
if (j%2)
|
|
{
|
|
lineCases[13] = dim0*2;
|
|
lineCases[21] = dim0*2;
|
|
lineCases[37] = dim0*2;
|
|
lineCases[63] = dim0*2;
|
|
isect1Ptr = isect1;
|
|
isect2Ptr = isect1 + dim0*2;
|
|
}
|
|
else
|
|
{
|
|
lineCases[13] = -dim0*2;
|
|
lineCases[21] = -dim0*2;
|
|
lineCases[37] = -dim0*2;
|
|
lineCases[63] = -dim0*2;
|
|
isect1Ptr = isect1 + dim0*2;
|
|
isect2Ptr = isect1;
|
|
}
|
|
|
|
for (i = min0; i < max0; i++)
|
|
{
|
|
s0 = s1;
|
|
s1 = *(inPtr + inc0);
|
|
// compute in/out for verts
|
|
v0 = (s0 < value ? 0 : 1);
|
|
v1 = (s1 < value ? 0 : 1);
|
|
// if we have an intersection
|
|
*isect2Ptr = -1;
|
|
*(isect2Ptr + 1) = -1;
|
|
if (v0 ^ v1)
|
|
{
|
|
// watch for degenerate points
|
|
if (s0 == value)
|
|
{
|
|
if (i > min0 && *(isect2Ptr-2) > -1)
|
|
{
|
|
*isect2Ptr = *(isect2Ptr-2);
|
|
}
|
|
else if (j > min1 && *(isect1Ptr+1) > -1)
|
|
{
|
|
*isect2Ptr = *(isect1Ptr+1);
|
|
}
|
|
}
|
|
else if (s1 == value && j > min1 && *(isect1Ptr+3) > -1)
|
|
{
|
|
*isect2Ptr = *(isect1Ptr+3);
|
|
}
|
|
// if the edge has not been set yet then it is a new point
|
|
if (*isect2Ptr == -1)
|
|
{
|
|
t = (value - s0) / (s1 - s0);
|
|
x[axis0] = origin[axis0] + spacing[axis0]*(i+t);
|
|
x[axis1] = y;
|
|
*isect2Ptr = newPts->InsertNextPoint(x);
|
|
if (newScalars)
|
|
{
|
|
newScalars->InsertNextTuple(&value);
|
|
}
|
|
}
|
|
}
|
|
if (j < max1)
|
|
{
|
|
s2 = *(inPtr + inc1);
|
|
v2 = (s2 < value ? 0 : 1);
|
|
if (v0 ^ v2)
|
|
{
|
|
if (s0 == value)
|
|
{
|
|
if (*isect2Ptr > -1)
|
|
{
|
|
*(isect2Ptr + 1) = *isect2Ptr;
|
|
}
|
|
else if (j > min1 && *(isect1Ptr+1) > -1)
|
|
{
|
|
*(isect2Ptr + 1) = *(isect1Ptr+1);
|
|
}
|
|
else if (i > min0 && *(isect2Ptr-2) > -1)
|
|
{
|
|
*(isect2Ptr + 1) = *(isect2Ptr-2);
|
|
}
|
|
}
|
|
// if the edge has not been set yet then it is a new point
|
|
if (*(isect2Ptr + 1) == -1)
|
|
{
|
|
t = (value - s0) / (s2 - s0);
|
|
x[axis0] = origin[axis0] + spacing[axis0]*i;
|
|
x[axis1] = y + spacing[axis1]*t;
|
|
*(isect2Ptr + 1) = newPts->InsertNextPoint(x);
|
|
if (newScalars)
|
|
{
|
|
newScalars->InsertNextTuple(&value);
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
if (j > min1)
|
|
{
|
|
// now add any lines that need to be added
|
|
// basically look at the isect values,
|
|
// form an index and lookup the lines
|
|
idx = (*isect1Ptr > -1 ? 8 : 0);
|
|
idx = idx + (*(isect1Ptr +1) > -1 ? 4 : 0);
|
|
idx = idx + (*(isect1Ptr +3) > -1 ? 2 : 0);
|
|
idx = idx + (*isect2Ptr > -1 ? 1 : 0);
|
|
tablePtr = lineCases + idx*4;
|
|
|
|
if (*tablePtr != -1)
|
|
{
|
|
ptIds[0] = *(isect1Ptr + *tablePtr);
|
|
tablePtr++;
|
|
ptIds[1] = *(isect1Ptr + *tablePtr);
|
|
if (ptIds[0] != ptIds[1])
|
|
{
|
|
// insert non-degenerate lines
|
|
lines->InsertNextCell(2,ptIds);
|
|
}
|
|
tablePtr++;
|
|
if (*tablePtr != -1)
|
|
{
|
|
ptIds[0] = *(isect1Ptr + *tablePtr);
|
|
tablePtr++;
|
|
ptIds[1] = *(isect1Ptr + *tablePtr);
|
|
if (ptIds[0] != ptIds[1])
|
|
{
|
|
lines->InsertNextCell(2,ptIds);
|
|
}
|
|
}
|
|
}
|
|
}
|
|
inPtr += inc0;
|
|
isect2Ptr += 2;
|
|
isect1Ptr += 2;
|
|
}
|
|
// now compute the last column, use s2 since it is around
|
|
if (j < max1)
|
|
{
|
|
s2 = *(inPtr + dim0);
|
|
v2 = (s2 < value ? 0 : 1);
|
|
*(isect2Ptr + 1) = -1;
|
|
if (v1 ^ v2)
|
|
{
|
|
if (s1 == value && *(isect2Ptr-2) > -1)
|
|
{
|
|
*(isect2Ptr + 1) = *(isect2Ptr-2);
|
|
}
|
|
else if (s1 == value && *(isect1Ptr+1) > -1)
|
|
{
|
|
*(isect2Ptr + 1) = *(isect1Ptr+1);
|
|
}
|
|
else
|
|
{
|
|
t = (value - s1) / (s2 - s1);
|
|
x[axis0] = origin[axis0] + spacing[axis0]*max0;
|
|
x[axis1] = y + spacing[axis1]*t;
|
|
*(isect2Ptr + 1) = newPts->InsertNextPoint(x);
|
|
if (newScalars)
|
|
{
|
|
newScalars->InsertNextTuple(&value);
|
|
}
|
|
}
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
delete [] isect1;
|
|
}
|
|
|
|
//----------------------------------------------------------------------------
|
|
//
|
|
// Contouring filter specialized for images (or slices from images)
|
|
//
|
|
int vtkSynchronizedTemplates2D::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
|
|
vtkImageData *input = vtkImageData::SafeDownCast(
|
|
inInfo->Get(vtkDataObject::DATA_OBJECT()));
|
|
vtkPolyData *output = vtkPolyData::SafeDownCast(
|
|
outInfo->Get(vtkDataObject::DATA_OBJECT()));
|
|
|
|
vtkPoints *newPts;
|
|
vtkCellArray *newLines;
|
|
vtkDataArray *inScalars;
|
|
vtkDataArray *newScalars = NULL;
|
|
int *ext;
|
|
int dims[3];
|
|
int dataSize, estimatedSize;
|
|
|
|
|
|
vtkDebugMacro(<< "Executing 2D structured contour");
|
|
|
|
ext =
|
|
inInfo->Get(vtkStreamingDemandDrivenPipeline::UPDATE_EXTENT());
|
|
inScalars = this->GetInputArrayToProcess(0,inputVector);
|
|
if ( inScalars == NULL )
|
|
{
|
|
vtkErrorMacro(<<"Scalars must be defined for contouring");
|
|
return 1;
|
|
}
|
|
|
|
int numComps = inScalars->GetNumberOfComponents();
|
|
if (this->ArrayComponent >= numComps)
|
|
{
|
|
vtkErrorMacro("Scalars have " << numComps << " components. "
|
|
"ArrayComponent must be smaller than " << numComps);
|
|
return 1;
|
|
}
|
|
|
|
// We have to compute the dimenisons from the update extent because
|
|
// the extent may be larger.
|
|
dims[0] = ext[1]-ext[0]+1;
|
|
dims[1] = ext[3]-ext[2]+1;
|
|
dims[2] = ext[5]-ext[4]+1;
|
|
|
|
//
|
|
// Check dimensionality of data and get appropriate form
|
|
//
|
|
dataSize = dims[0] * dims[1] * dims[2];
|
|
|
|
//
|
|
// Allocate necessary objects
|
|
//
|
|
estimatedSize = (int) (sqrt((double)(dataSize)));
|
|
if (estimatedSize < 1024)
|
|
{
|
|
estimatedSize = 1024;
|
|
}
|
|
newPts = vtkPoints::New();
|
|
newPts->Allocate(estimatedSize,estimatedSize);
|
|
newLines = vtkCellArray::New();
|
|
newLines->Allocate(newLines->EstimateSize(estimatedSize,2));
|
|
|
|
//
|
|
// Check data type and execute appropriate function
|
|
//
|
|
|
|
void *scalars = inScalars->GetVoidPointer(0);
|
|
if (this->ComputeScalars)
|
|
{
|
|
newScalars = inScalars->NewInstance();
|
|
newScalars->SetNumberOfComponents(inScalars->GetNumberOfComponents());
|
|
newScalars->SetName(inScalars->GetName());
|
|
newScalars->Allocate(5000,25000);
|
|
}
|
|
switch (inScalars->GetDataType())
|
|
{
|
|
vtkTemplateMacro(
|
|
vtkContourImage(this,(VTK_TT *)scalars, newPts,
|
|
newScalars, newLines, input, ext));
|
|
}//switch
|
|
|
|
// Lets set the name of the scalars here.
|
|
if (newScalars)
|
|
{
|
|
newScalars->SetName(inScalars->GetName());
|
|
}
|
|
|
|
vtkDebugMacro(<<"Created: "
|
|
<< newPts->GetNumberOfPoints() << " points, "
|
|
<< newLines->GetNumberOfCells() << " lines");
|
|
//
|
|
// Update ourselves. Because we don't know up front how many lines
|
|
// we've created, take care to reclaim memory.
|
|
//
|
|
output->SetPoints(newPts);
|
|
newPts->Delete();
|
|
|
|
output->SetLines(newLines);
|
|
newLines->Delete();
|
|
|
|
if (newScalars)
|
|
{
|
|
int idx = output->GetPointData()->AddArray(newScalars);
|
|
output->GetPointData()->SetActiveAttribute(idx, vtkDataSetAttributes::SCALARS);
|
|
newScalars->Delete();
|
|
}
|
|
|
|
output->Squeeze();
|
|
return 1;
|
|
}
|
|
|
|
//----------------------------------------------------------------------------
|
|
int vtkSynchronizedTemplates2D::FillInputPortInformation(int, vtkInformation *info)
|
|
{
|
|
info->Set(vtkAlgorithm::INPUT_REQUIRED_DATA_TYPE(), "vtkImageData");
|
|
return 1;
|
|
}
|
|
|
|
//----------------------------------------------------------------------------
|
|
void vtkSynchronizedTemplates2D::PrintSelf(ostream& os, vtkIndent indent)
|
|
{
|
|
this->Superclass::PrintSelf(os,indent);
|
|
|
|
this->ContourValues->PrintSelf(os,indent.GetNextIndent());
|
|
if (this->ComputeScalars)
|
|
{
|
|
os << indent << "ComputeScalarsOn\n";
|
|
}
|
|
else
|
|
{
|
|
os << indent << "ComputeScalarsOff\n";
|
|
}
|
|
os << indent << "ArrayComponent: " << this->ArrayComponent << endl;
|
|
}
|
|
|