Commit 74575a55 authored by Will Schroeder's avatar Will Schroeder

Now interpolating point attribute data on request

Flying edges is now interpolating point attribute data when the
flag InterpolateAttributes is set.

Also made some minor performance tweaks to make up for the extra work.
parent 1521bd0a
......@@ -9,6 +9,7 @@ set(Module_SRCS
vtkAnnotation.cxx
vtkAnnotationLayers.cxx
vtkArrayData.cxx
vtkArrayListTemplate.txx
vtkAttributesErrorMetric.cxx
vtkBiQuadraticQuad.cxx
vtkBiQuadraticQuadraticHexahedron.cxx
......@@ -213,6 +214,7 @@ set(Module_SRCS
)
set(${vtk-module}_HDRS
vtkArrayListTemplate.h
vtkCellType.h
vtkMappedUnstructuredGrid.h
vtkMappedUnstructuredGridCellIterator.h
......@@ -260,6 +262,7 @@ set_source_files_properties(
set_source_files_properties(
vtkAMRBox
vtkArrayListTemplate.txx
vtkAtom
vtkBond
vtkBoundingBox
......
......@@ -26,10 +26,10 @@
// vtkDataSetAttributes::CopyInterpolate() or InterpolateAllocate() which
// performs the initial magic of constructing input and output arrays. Then
// the input attributes, and output attributes, are passed to initialize the
// internal structures. Internally, pairs of typed arrays are created; the
// operations (e.g., interpolate) occur on these typed arrays using a
// typeless, virtual dispatch base class.
// internal structures. Essentially these internal structures are pairs of
// arrays of the same type, which can be efficently accessed and
// assigned. The operations on these array pairs (e.g., interpolation) occur
// using a typeless, virtual dispatch base class.
// .SECTION See Also
// vtkFieldData vtkDataSetAttributes vtkPointData vtkCellData
......@@ -46,8 +46,10 @@ struct BaseArrayPair
{
vtkIdType Num;
int NumComp;
vtkDataArray *OutputArray;
BaseArrayPair(vtkIdType num, int numComp) : Num(num), NumComp(numComp)
BaseArrayPair(vtkIdType num, int numComp, vtkDataArray *outArray) :
Num(num), NumComp(numComp), OutputArray(outArray)
{
}
virtual ~BaseArrayPair()
......@@ -56,8 +58,11 @@ struct BaseArrayPair
virtual void Copy(vtkIdType inId, vtkIdType outId) = 0;
virtual void Interpolate(int numWeights, const vtkIdType *ids,
const double *weights, vtkIdType outPtId) = 0;
const double *weights, vtkIdType outId) = 0;
virtual void InterpolateEdge(vtkIdType v0, vtkIdType v1,
double t, vtkIdType outId) = 0;
virtual void AssignNullValue(vtkIdType outId) = 0;
virtual void Realloc(vtkIdType sze) = 0;
};
// Type specific interpolation on a matched pair of data arrays
......@@ -68,8 +73,8 @@ struct ArrayPair : public BaseArrayPair
T *Output;
T NullValue;
ArrayPair(T *in, T *out, vtkIdType num, int numComp, T null) :
BaseArrayPair(num,numComp), Input(in), Output(out), NullValue(null)
ArrayPair(T *in, T *out, vtkIdType num, int numComp, vtkDataArray *outArray, T null) :
BaseArrayPair(num,numComp,outArray), Input(in), Output(out), NullValue(null)
{
}
virtual ~ArrayPair() //calm down some finicky compilers
......@@ -85,7 +90,7 @@ struct ArrayPair : public BaseArrayPair
}
virtual void Interpolate(int numWeights, const vtkIdType *ids,
const double *weights, vtkIdType outPtId)
const double *weights, vtkIdType outId)
{
for (int j=0; j < this->NumComp; ++j)
{
......@@ -94,7 +99,19 @@ struct ArrayPair : public BaseArrayPair
{
v += weights[i] * static_cast<double>(this->Input[ids[i]*this->NumComp+j]);
}
this->Output[outPtId*this->NumComp+j] = static_cast<T>(v);
this->Output[outId*this->NumComp+j] = static_cast<T>(v);
}
}
virtual void InterpolateEdge(vtkIdType v0, vtkIdType v1, double t, vtkIdType outId)
{
double v;
vtkIdType numComp=this->NumComp;
for (int j=0; j < numComp; ++j)
{
v = this->Input[v0*numComp+j] +
t * (this->Input[v1*numComp+j] - this->Input[v0*numComp+j]);
this->Output[outId*numComp+j] = static_cast<T>(v);
}
}
......@@ -106,6 +123,12 @@ struct ArrayPair : public BaseArrayPair
}
}
virtual void Realloc(vtkIdType sze)
{
this->OutputArray->WriteVoidPointer(0,sze*this->NumComp);
this->Output = static_cast<T*>(this->OutputArray->GetVoidPointer(0));
}
};
// Forward declarations. This makes working with vtkTemplateMacro easier.
......@@ -137,16 +160,26 @@ struct ArrayList
}
// Loop over the arrays and have them interpolate themselves
void Interpolate(int numWeights, const vtkIdType *ids, const double *weights, vtkIdType outPtId)
void Interpolate(int numWeights, const vtkIdType *ids, const double *weights, vtkIdType outId)
{
for (std::vector<BaseArrayPair*>::iterator it = Arrays.begin();
it != Arrays.end(); ++it)
{
(*it)->Interpolate(numWeights, ids, weights, outPtId);
(*it)->Interpolate(numWeights, ids, weights, outId);
}
}
// Loop over the arrays and have them interpolate themselves
// Loop over the arrays perform edge interpolation
void InterpolateEdge(vtkIdType v0, vtkIdType v1, double t, vtkIdType outId)
{
for (std::vector<BaseArrayPair*>::iterator it = Arrays.begin();
it != Arrays.end(); ++it)
{
(*it)->InterpolateEdge(v0, v1, t, outId);
}
}
// Loop over the arrays and assign the null value
void AssignNullValue(vtkIdType outId)
{
for (std::vector<BaseArrayPair*>::iterator it = Arrays.begin();
......@@ -156,6 +189,16 @@ struct ArrayList
}
}
// Extend (realloc) the arrays
void Realloc(vtkIdType sze)
{
for (std::vector<BaseArrayPair*>::iterator it = Arrays.begin();
it != Arrays.end(); ++it)
{
(*it)->Realloc(sze);
}
}
// Only you can prevent memory leaks!
~ArrayList()
{
......
......@@ -22,9 +22,9 @@
// Sort of a little object factory (in conjunction w/ vtkTemplateMacro())
template <typename T>
void CreateArrayPair(ArrayList *list, T *inData, T *outData,
vtkIdType numPts, int numComp, T nullValue)
vtkIdType numPts, int numComp, vtkDataArray *outArray, T nullValue)
{
ArrayPair<T> *pair = new ArrayPair<T>(inData,outData,numPts,numComp, nullValue);
ArrayPair<T> *pair = new ArrayPair<T>(inData,outData,numPts,numComp,outArray,nullValue);
list->Arrays.push_back(pair);
}
......@@ -64,7 +64,7 @@ AddArrays(vtkIdType numOutPts, vtkDataSetAttributes *inPD, vtkDataSetAttributes
{
vtkTemplateMacro(CreateArrayPair(this, static_cast<VTK_TT *>(iD),
static_cast<VTK_TT *>(oD),numOutPts,oNumComp,
static_cast<VTK_TT>(nullValue)));
oArray,static_cast<VTK_TT>(nullValue)));
}//over all VTK types
}//if matching types
}//if matching input array
......
......@@ -18,6 +18,7 @@ vtk_add_test_python(
TestElevationFilter.py
TestFlyingEdges2D.py
TestFlyingEdges3D.py
TestFlyingEdges3DWithInterpolation.py
TestFlyingEdgesExtents.py
TestFlyingEdgesPlaneCutter.py
TestGridSynchronizedTemplates3D.py
......
......@@ -31,6 +31,8 @@ iso.SetInputConnection(sample.GetOutputPort())
iso.SetValue(0,0.25)
iso.ComputeNormalsOn()
iso.ComputeGradientsOn()
iso.ComputeScalarsOn()
iso.InterpolateAttributesOff()
isoMapper = vtk.vtkPolyDataMapper()
isoMapper.SetInputConnection(iso.GetOutputPort())
......
#!/usr/bin/env python
import vtk
from vtk.test import Testing
from vtk.util.misc import vtkGetDataRoot
VTK_DATA_ROOT = vtkGetDataRoot()
# Create the RenderWindow, Renderer and both Actors
#
ren1 = vtk.vtkRenderer()
renWin = vtk.vtkRenderWindow()
renWin.SetMultiSamples(0)
renWin.AddRenderer(ren1)
iren = vtk.vtkRenderWindowInteractor()
iren.SetRenderWindow(renWin)
# Create a synthetic source
sphere = vtk.vtkSphere()
sphere.SetCenter( 0.0,0.0,0.0)
sphere.SetRadius(0.25)
# Iso-surface to create geometry. Demonstrate the ability to
# interpolate data attributes.
sample = vtk.vtkSampleFunction()
sample.SetImplicitFunction(sphere)
sample.SetModelBounds(-0.5,0.5, -0.5,0.5, -0.5,0.5)
sample.SetSampleDimensions(100,100,100)
# Now create some new attributes
cyl = vtk.vtkCylinder()
cyl.SetRadius(0.1)
cyl.SetAxis(1,1,1)
attr = vtk.vtkSampleImplicitFunctionFilter()
attr.SetInputConnection(sample.GetOutputPort())
attr.SetImplicitFunction(cyl)
attr.ComputeGradientsOn()
attr.Update()
iso = vtk.vtkFlyingEdges3D()
iso.SetInputConnection(attr.GetOutputPort())
iso.SetInputArrayToProcess(0, 0, 0, vtk.vtkDataObject.FIELD_ASSOCIATION_POINTS, "scalars")
iso.SetValue(0,0.25)
iso.ComputeNormalsOn()
iso.ComputeGradientsOn()
iso.ComputeScalarsOn()
iso.InterpolateAttributesOn()
# Time execution
timer = vtk.vtkTimerLog()
timer.StartTimer()
iso.Update()
timer.StopTimer()
time = timer.GetElapsedTime()
print("Flying edges with attributes: {0}".format(time))
isoMapper = vtk.vtkPolyDataMapper()
isoMapper.SetInputConnection(iso.GetOutputPort())
isoMapper.ScalarVisibilityOn()
isoMapper.SetScalarModeToUsePointFieldData()
isoMapper.SelectColorArray("Implicit scalars")
isoMapper.SetScalarRange(0,.3)
isoActor = vtk.vtkActor()
isoActor.SetMapper(isoMapper)
isoActor.GetProperty().SetColor(1,1,1)
isoActor.GetProperty().SetOpacity(1)
outline = vtk.vtkOutlineFilter()
outline.SetInputConnection(sample.GetOutputPort())
outlineMapper = vtk.vtkPolyDataMapper()
outlineMapper.SetInputConnection(outline.GetOutputPort())
outlineActor = vtk.vtkActor()
outlineActor.SetMapper(outlineMapper)
outlineProp = outlineActor.GetProperty()
# Add the actors to the renderer, set the background and size
#
ren1.AddActor(outlineActor)
ren1.AddActor(isoActor)
ren1.SetBackground(0,0,0)
renWin.SetSize(300,300)
ren1.ResetCamera()
iren.Initialize()
renWin.Render()
# --- end of script --
......@@ -26,6 +26,7 @@
#include "vtkFloatArray.h"
#include "vtkStreamingDemandDrivenPipeline.h"
#include "vtkMarchingCubesTriangleCases.h"
#include "vtkArrayListTemplate.h"
#include "vtkSMPTools.h"
#include <cmath>
......@@ -33,6 +34,9 @@
vtkStandardNewMacro(vtkFlyingEdges3D);
//----------------------------------------------------------------------------
namespace {
#include "vtkArrayListTemplate.h" // For processing attribute data
}
// This templated class implements the heart of the algorithm.
// vtkFlyingEdges3D populates the information in this class and
......@@ -121,6 +125,8 @@ public:
float *NewGradients;
float *NewNormals;
bool NeedGradients;
bool InterpolateAttributes;
ArrayList Arrays;
// Setup algorithm
vtkFlyingEdges3DAlgorithm();
......@@ -212,7 +218,8 @@ public:
const int incs[3],
float x1[3],
vtkIdType vId,
vtkIdType ijk[3],
vtkIdType ijk0[3],
vtkIdType ijk1[3],
float g0[3])
{
......@@ -224,7 +231,7 @@ public:
if ( this->NeedGradients )
{
float gTmp[3], g1[3];
this->ComputeGradient(loc,ijk,
this->ComputeGradient(loc,ijk1,
s + incs[0], s - incs[0],
s + incs[1], s - incs[1],
s + incs[2], s - incs[2],
......@@ -244,6 +251,13 @@ public:
vtkMath::Normalize(n);
}
}//if normals or gradients required
if ( this->InterpolateAttributes )
{
vtkIdType v0=ijk0[0] + ijk0[1]*incs[1] + ijk0[2]*incs[2];
vtkIdType v1=ijk1[0] + ijk1[1]*incs[1] + ijk1[2]*incs[2];;
this->Arrays.InterpolateEdge(v0,v1,t,vId);
}
}
// Compute the gradient on a point which may be on the boundary of the volume.
......@@ -380,7 +394,7 @@ public:
// Interface between VTK and templated functions
static void Contour(vtkFlyingEdges3D *self, vtkImageData *input,
int extent[6], vtkIdType *incs, T *scalars,
vtkPoints *newPts, vtkCellArray *newTris,
vtkPolyData *output, vtkPoints *newPts, vtkCellArray *newTris,
vtkDataArray *newScalars,vtkFloatArray *newNormals,
vtkFloatArray *newGradients);
};
......@@ -683,6 +697,13 @@ InterpolateEdge(double value, vtkIdType ijk[3],
vtkMath::Normalize(n);
}
}//if normals or gradients required
if ( this->InterpolateAttributes )
{
vtkIdType v0=ijk0[0] + ijk0[1]*incs[1] + ijk0[2]*incs[2];
vtkIdType v1=ijk1[0] + ijk1[1]*incs[1] + ijk1[2]*incs[2];;
this->Arrays.InterpolateEdge(v0,v1,t,vId);
}
}
//----------------------------------------------------------------------------
......@@ -719,7 +740,7 @@ GeneratePoints(double value, unsigned char loc, vtkIdType ijk[3],
T const * const sPtr2 = (sPtr+incs[i]);
double t = (value - *sPtr) / (*sPtr2 - *sPtr);
this->InterpolateAxesEdge(t, loc, x, sPtr2, incs, x1, eIds[i*4], ijk1, g0);
this->InterpolateAxesEdge(t, loc, x, sPtr2, incs, x1, eIds[i*4], ijk, ijk1, g0);
}
}
......@@ -794,15 +815,14 @@ ProcessXEdge(double value, T const* const inPtr, vtkIdType row, vtkIdType slice)
vtkIdType nxcells=this->Dims[0]-1;
vtkIdType minInt=nxcells, maxInt = 0;
vtkIdType *edgeMetaData;
unsigned char *ePtr = this->XCases + slice*this->SliceOffset + row*nxcells;
unsigned char edgeCase, *ePtr=this->XCases+slice*this->SliceOffset+row*nxcells;
double s0, s1 = static_cast<double>(*inPtr);
vtkIdType sum = 0;
//run along the entire x-edge computing edge cases
edgeMetaData = this->EdgeMetaData + (slice*this->Dims[1] + row)*6;
std::fill_n(edgeMetaData, 6, 0);
vtkIdType sum = 0;
//pull this out help reduce false sharing
vtkIdType inc0 = this->Inc0;
......@@ -811,11 +831,14 @@ ProcessXEdge(double value, T const* const inPtr, vtkIdType row, vtkIdType slice)
s0 = s1;
s1 = static_cast<double>(*(inPtr + (i+1)*inc0));
unsigned char edgeCase = vtkFlyingEdges3DAlgorithm::Below;
if (s0 >= value)
{
edgeCase = vtkFlyingEdges3DAlgorithm::LeftAbove;
}
else
{
edgeCase = vtkFlyingEdges3DAlgorithm::Below;
}
if( s1 >= value)
{
edgeCase |= vtkFlyingEdges3DAlgorithm::RightAbove;
......@@ -828,7 +851,10 @@ ProcessXEdge(double value, T const* const inPtr, vtkIdType row, vtkIdType slice)
edgeCase == vtkFlyingEdges3DAlgorithm::RightAbove )
{
++sum; //increment number of intersections along x-edge
minInt = ( i < minInt ? i : minInt);
if ( i < minInt )
{
minInt = i;
}
maxInt = i + 1;
}//if contour interacts with this x-edge
}//for all x-cell edges along this x-edge
......@@ -936,6 +962,7 @@ ProcessYZEdges(vtkIdType row, vtkIdType slice)
// voxel axes. Also check the number of primitives generated.
unsigned char *edgeUses, eCase, numTris;
ePtr[0] += xL; ePtr[1] += xL; ePtr[2] += xL; ePtr[3] += xL;
const vtkIdType dim0Wall = this->Dims[0]-2;
for (i=xL; i < xR; ++i) //run along the trimmed x-voxels
{
eCase = this->GetEdgeCase(ePtr);
......@@ -950,7 +977,7 @@ ProcessYZEdges(vtkIdType row, vtkIdType slice)
edgeUses = this->GetEdgeUses(eCase);
eMD[0][1] += edgeUses[4]; //y-voxel axes edge always counted
eMD[0][2] += edgeUses[8]; //z-voxel axes edge always counted
loc = yzLoc | (i >= (this->Dims[0]-2) ? MaxBoundary : Interior);
loc = yzLoc | (i >= dim0Wall ? MaxBoundary : Interior);
if ( loc != 0 )
{
this->CountBoundaryYZInts(loc,edgeUses,eMD);
......@@ -1030,6 +1057,8 @@ GenerateOutput(double value, T* rowPtr, vtkIdType row, vtkIdType slice)
//load the inc0/inc1/inc2 into local memory
const int incs[3] = { this->Inc0, this->Inc1, this->Inc2 };
const T* sPtr = rowPtr + xL*incs[0];
const double xSpace = this->Spacing[0];
const vtkIdType dim0Wall = this->Dims[0]-2;
for (i=xL; i < xR; ++i)
{
......@@ -1042,7 +1071,7 @@ GenerateOutput(double value, T* rowPtr, vtkIdType row, vtkIdType slice)
// Now generate point(s) along voxel axes if needed. Remember to take
// boundary into account.
loc = yzLoc | (i < 1 ? MinBoundary :
(i >= (this->Dims[0]-2) ? MaxBoundary : Interior));
(i >= dim0Wall ? MaxBoundary : Interior));
if ( this->CaseIncludesAxes(eCase) || loc != Interior )
{
unsigned char const * const edgeUses = this->GetEdgeUses(eCase);
......@@ -1059,7 +1088,7 @@ GenerateOutput(double value, T* rowPtr, vtkIdType row, vtkIdType slice)
++ijk[0];
sPtr += incs[0];
x[0]+= this->Spacing[0];
x[0] += xSpace;
} //for all non-trimmed cells along this x-edge
}
......@@ -1069,9 +1098,9 @@ GenerateOutput(double value, T* rowPtr, vtkIdType row, vtkIdType slice)
// class. It also invokes the three passes of the Flying Edges algorithm.
template <class T> void vtkFlyingEdges3DAlgorithm<T>::
Contour(vtkFlyingEdges3D *self, vtkImageData *input, int extent[6],
vtkIdType *incs, T *scalars, vtkPoints *newPts, vtkCellArray *newTris,
vtkDataArray *newScalars, vtkFloatArray *newNormals,
vtkFloatArray *newGradients)
vtkIdType *incs, T *scalars, vtkPolyData *output, vtkPoints *newPts,
vtkCellArray *newTris, vtkDataArray *newScalars,
vtkFloatArray *newNormals, vtkFloatArray *newGradients)
{
double value, *values = self->GetValues();
int numContours = self->GetNumberOfContours();
......@@ -1114,6 +1143,12 @@ Contour(vtkFlyingEdges3D *self, vtkImageData *input, int extent[6],
// for computational trimming).
algo.EdgeMetaData = new vtkIdType [algo.NumberOfEdges*6];
// Interpolating attributes and other stuff. Interpolate extra attributes only if they
// exist and the user requests it.
algo.NeedGradients = (newGradients || newNormals);
algo.InterpolateAttributes = (self->GetInterpolateAttributes() &&
input->GetPointData()->GetNumberOfArrays() > 1) ? true : false;
// Loop across each contour value. This encompasses all three passes.
for (vidx = 0; vidx < numContours; vidx++)
{
......@@ -1195,7 +1230,19 @@ Contour(vtkFlyingEdges3D *self, vtkImageData *input, int extent[6],
newNormals->WriteVoidPointer(0,3*totalPts);
algo.NewNormals = static_cast<float*>(newNormals->GetVoidPointer(0));
}
algo.NeedGradients = (algo.NewGradients || algo.NewNormals);
if ( algo.InterpolateAttributes )
{
if ( vidx == 0 ) //first contour
{
// output->GetPointData()->CopyScalarsOff();
output->GetPointData()->InterpolateAllocate(input->GetPointData(),totalPts);
algo.Arrays.AddArrays(totalPts,input->GetPointData(),output->GetPointData());
}
else
{
algo.Arrays.Realloc(totalPts);
}
}
// PASS 4: Fourth and final pass: Process voxel rows and generate output.
// Note that we are simultaneously generating triangles and interpolating
......@@ -1226,6 +1273,7 @@ vtkFlyingEdges3D::vtkFlyingEdges3D()
this->ComputeNormals = 1;
this->ComputeGradients = 0;
this->ComputeScalars = 1;
this->InterpolateAttributes = 0;
this->ArrayComponent = 0;
// by default process active point scalars
......@@ -1365,7 +1413,7 @@ int vtkFlyingEdges3D::RequestData(
{
vtkTemplateMacro(vtkFlyingEdges3DAlgorithm<VTK_TT>::
Contour(this, input, exExt, incs, (VTK_TT *)ptr,
newPts, newTris, newScalars, newNormals,
output, newPts, newTris, newScalars, newNormals,
newGradients));
}
......@@ -1422,5 +1470,6 @@ void vtkFlyingEdges3D::PrintSelf(ostream& os, vtkIndent indent)
os << indent << "Compute Normals: " << (this->ComputeNormals ? "On\n" : "Off\n");
os << indent << "Compute Gradients: " << (this->ComputeGradients ? "On\n" : "Off\n");
os << indent << "Compute Scalars: " << (this->ComputeScalars ? "On\n" : "Off\n");
os << indent << "Interpolate Attributes: " << (this->InterpolateAttributes ? "On\n" : "Off\n");
os << indent << "ArrayComponent: " << this->ArrayComponent << endl;
}
......@@ -100,6 +100,15 @@ public:
vtkGetMacro(ComputeScalars,int);
vtkBooleanMacro(ComputeScalars,int);
// Description:
// Indicate whether to interpolate other attribute data. That is, as the
// isosurface is generated, interpolate all point attribute data across
// the edge. This is independent of scalar interpolation, which is
// controlled by the ComputeScalars flag.
vtkSetMacro(InterpolateAttributes,int);
vtkGetMacro(InterpolateAttributes,int);
vtkBooleanMacro(InterpolateAttributes,int);
// Description:
// Set a particular contour value at contour number i. The index i ranges
// between 0<=i<NumberOfContours.
......@@ -157,6 +166,7 @@ protected:
int ComputeNormals;
int ComputeGradients;
int ComputeScalars;
int InterpolateAttributes;
int ArrayComponent;
vtkContourValues *ContourValues;
......
set(Module_SRCS
vtkArrayListTemplate.txx
vtkEllipsoidalGaussianKernel.cxx
vtkGaussianKernel.cxx
vtkInterpolationKernel.cxx
......@@ -10,17 +9,9 @@ set(Module_SRCS
vtkVoronoiKernel.cxx
)
set(${vtk-module}_HDRS
vtkArrayListTemplate.h
)
set_source_files_properties(
vtkInterpolationKernel
ABSTRACT
)
set_source_files_properties(
vtkArrayListTemplate.txx
WRAP_EXCLUDE
)
vtk_module_library(vtkFiltersPoints ${Module_SRCS})
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment