Commit bc793aec authored by Will Schroeder's avatar Will Schroeder Committed by Kitware Robot
Browse files

Merge topic 'FlyingEdges-CornerCases'

43e4b1c6 Handle multiple contour values properly
ac2fb754 New valid regression test image
2f1df1ad

 Resolved merge conflict
Acked-by: Kitware Robot's avatarKitware Robot <kwrobot@kitware.com>
Reviewed-by: default avatarDavid E DeMarle <dave.demarle@kitware.com>
Merge-request: !779
parents f0011556 43e4b1c6
......@@ -48,6 +48,7 @@ int TestFlyingEdges(int argc, char *argv[])
vtkNew<vtkPolyDataMapper> mapper;
mapper->SetInputConnection(flyingEdges->GetOutputPort());
mapper->SetScalarRange(128,225);
vtkNew<vtkActor> actor;
actor->SetMapper(mapper.GetPointer());
vtkNew<vtkRenderer> ren;
......
88377f459808402531ecc5a555ee8d20
b59c1892d5b4743818906e989b11324f
......@@ -18,6 +18,7 @@ vtk_add_test_python(
TestElevationFilter.py
TestFlyingEdges2D.py
TestFlyingEdges3D.py
TestFlyingEdgesExtents.py
TestFlyingEdgesPlaneCutter.py
TestGridSynchronizedTemplates3D.py
TestMarchingSquares.py
......
#!/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()
ren2 = vtk.vtkRenderer()
renWin = vtk.vtkRenderWindow()
renWin.SetMultiSamples(0)
renWin.AddRenderer(ren1)
renWin.AddRenderer(ren2)
ren1.SetViewport(0,0,0.5,1)
ren2.SetViewport(0.5,0,1,1)
iren = vtk.vtkRenderWindowInteractor()
iren.SetRenderWindow(renWin)
# Test the negative extents
source = vtk.vtkRTAnalyticSource()
iso = vtk.vtkFlyingEdges3D()
iso.SetInputConnection(source.GetOutputPort())
iso.SetValue(0,150)
isoMapper = vtk.vtkPolyDataMapper()
isoMapper.SetInputConnection(iso.GetOutputPort())
isoMapper.ScalarVisibilityOff()
isoActor = vtk.vtkActor()
isoActor.SetMapper(isoMapper)
isoActor.GetProperty().SetColor(1,1,1)
isoActor.GetProperty().SetOpacity(1)
outline = vtk.vtkOutlineFilter()
outline.SetInputConnection(source.GetOutputPort())
outlineMapper = vtk.vtkPolyDataMapper()
outlineMapper.SetInputConnection(outline.GetOutputPort())
outlineActor = vtk.vtkActor()
outlineActor.SetMapper(outlineMapper)
outlineProp = outlineActor.GetProperty()
# Test the non-x-edge-intersecting contours and cuts
res = 100
plane = vtk.vtkPlane()
plane.SetOrigin(0,0,0)
plane.SetNormal(0,1,0)
sample = vtk.vtkSampleFunction()
sample.SetImplicitFunction(plane)
sample.SetModelBounds(-10,10, -10,10, -10,10)
sample.SetSampleDimensions(res,res,res)
sample.SetOutputScalarTypeToFloat()
sample.Update()
iso2 = vtk.vtkFlyingEdges3D()
iso2.SetInputConnection(sample.GetOutputPort())
iso2.SetValue(0,0.0)
isoMapper2 = vtk.vtkPolyDataMapper()
isoMapper2.SetInputConnection(iso2.GetOutputPort())
isoMapper2.ScalarVisibilityOff()
isoActor2 = vtk.vtkActor()
isoActor2.SetMapper(isoMapper2)
isoActor2.GetProperty().SetColor(1,1,1)
isoActor2.GetProperty().SetOpacity(1)
outline2 = vtk.vtkOutlineFilter()
outline2.SetInputConnection(sample.GetOutputPort())
outlineMapper2 = vtk.vtkPolyDataMapper()
outlineMapper2.SetInputConnection(outline2.GetOutputPort())
outlineActor2 = vtk.vtkActor()
outlineActor2.SetMapper(outlineMapper)
# Extract a slice in the desired orientation
center = [0.0, 0.0, 0.0]
axial = vtk.vtkMatrix4x4()
axial.DeepCopy((1, 0, 0, center[0],
0, 1, 0, center[1],
0, 0, 1, center[2],
0, 0, 0, 1))
reslice = vtk.vtkImageReslice()
reslice.SetInputConnection(sample.GetOutputPort())
reslice.SetOutputDimensionality(2)
reslice.SetResliceAxes(axial)
reslice.SetInterpolationModeToLinear()
iso3 = vtk.vtkFlyingEdges2D()
iso3.SetInputConnection(reslice.GetOutputPort())
iso3.SetValue(0,0.0)
tube = vtk.vtkTubeFilter()
tube.SetInputConnection(iso3.GetOutputPort())
tube.SetRadius(0.25)
mapper3 = vtk.vtkPolyDataMapper()
mapper3.SetInputConnection(tube.GetOutputPort())
actor3 = vtk.vtkActor()
actor3.SetMapper(mapper3)
# Add the actors to the renderer, set the background and size
#
renWin.SetSize(600,300)
ren1.AddActor(outlineActor)
ren1.AddActor(isoActor)
ren1.SetBackground(0,0,0)
ren1.ResetCamera()
cam = vtk.vtkCamera()
cam.SetPosition(0,1,0)
cam.SetFocalPoint(0,0,0)
cam.SetViewUp(0,0,1)
ren2.SetActiveCamera(cam)
ren2.AddActor(outlineActor2)
ren2.AddActor(isoActor2)
ren2.AddActor(actor3)
ren2.SetBackground(0,0,0)
ren2.ResetCamera()
iren.Initialize()
renWin.Render()
# --- end of script --
......@@ -37,7 +37,7 @@ class vtkFlyingEdges2DAlgorithm
{
public:
// Edge case table values.
enum {
enum EdgeClass {
Below = 0, //below isovalue
Above = 1, //above isovalue
LeftAbove = 1, //left vertex is above isovalue
......@@ -46,7 +46,7 @@ public:
};
// Dealing with boundary situations when processing images.
enum {
enum CellClass {
Interior = 0,
MinBoundary = 1,
MaxBoundary = 2
......@@ -81,9 +81,9 @@ public:
// Internal variables used by the various algorithm methods. Interfaces VTK
// image data in a form more convenient to the algorithm.
vtkIdType Dims[2];
double *Origin;
double *Spacing;
double X;
double Origin[3];
double Spacing[3];
double Z;
int Axis0;
int Min0;
int Max0;
......@@ -103,6 +103,14 @@ public:
// Instantiate and initialize key data members.
vtkFlyingEdges2DAlgorithm();
// Adjust the origin to the lower-left corner of the volume (if necessary)
void AdjustOrigin(int updateExt[6])
{
this->Origin[0] = this->Origin[0] + this->Spacing[0]*updateExt[0];
this->Origin[1] = this->Origin[1] + this->Spacing[1]*updateExt[2];
this->Origin[2] = this->Origin[2] + this->Spacing[2]*updateExt[4];;
}
// The three passes of the algorithm.
void ProcessXEdge(double value, T* inPtr, vtkIdType row); //PASS 1
void ProcessYEdges(vtkIdType row); //PASS 2
......@@ -164,7 +172,7 @@ public:
float *x = this->NewPoints + 3*vId;
x[0] = x0[0] + t*(x1[0]-x0[0]);
x[1] = x0[1] + t*(x1[1]-x0[1]);
x[2] = this->X;
x[2] = this->Z;
}
// Interpolate along an arbitrary edge, typically one that may be on the
......@@ -343,7 +351,7 @@ InterpolateEdge(double value, T *s, float x[3],unsigned char edgeNum,
float *xPtr = this->NewPoints + 3*vId;
xPtr[0] = x0[0] + t*(x1[0]-x0[0]);
xPtr[1] = x0[1] + t*(x1[1]-x0[1]);
xPtr[2] = this->X;
xPtr[2] = this->Z;
}
//----------------------------------------------------------------------------
......@@ -446,7 +454,7 @@ template <class T> void vtkFlyingEdges2DAlgorithm<T>::
ProcessYEdges(vtkIdType row)
{
// Grab the two edge cases bounding this pixel x-row.
unsigned char *ePtr0, *ePtr1, ec0, ec1;
unsigned char *ePtr0, *ePtr1, ec0, ec1, xInts=1;
ePtr0 = this->XCases + row*(this->Dims[0]-1);
ePtr1 = ePtr0 + this->Dims[0]-1;
......@@ -461,7 +469,11 @@ ProcessYEdges(vtkIdType row)
{
if ( *ePtr0 == *ePtr1 )
{
return; //there are no y-ints, thus no contour, skip pixel row
return; //there are no x- or y-ints, thus no contour, skip pixel row
}
else
{
xInts = 0; //there are y- edge ints however
}
}
......@@ -474,24 +486,32 @@ ProcessYEdges(vtkIdType row)
// the top and bottom rows of x-edges (without intersecting x-edges).
vtkIdType xL = ( (eMD0[3] < eMD1[3]) ? eMD0[3] : eMD1[3]);
vtkIdType xR = ( (eMD0[4] > eMD1[4]) ? eMD0[4] : eMD1[4]);
if ( xL > 0 )
if ( xInts )
{
ec0 = *(ePtr0 + xL);
ec1 = *(ePtr1 + xL);
if ( (ec0 & 0x1) != (ec1 & 0x1) )
if ( xL > 0 )
{
xL = eMD0[3] = 0; //reset left trim
ec0 = *(ePtr0 + xL);
ec1 = *(ePtr1 + xL);
if ( (ec0 & 0x1) != (ec1 & 0x1) )
{
xL = eMD0[3] = 0; //reset left trim
}
}
}
if ( xR < (this->Dims[0]-1) )
{
ec0 = *(ePtr0 + xR);
ec1 = *(ePtr1 + xR);
if ( (ec0 & 0x2) != (ec1 & 0x2) )
if ( xR < (this->Dims[0]-1) )
{
xR = eMD0[4] = this->Dims[0] - 1; //reset right trim
ec0 = *(ePtr0 + xR);
ec1 = *(ePtr1 + xR);
if ( (ec0 & 0x2) != (ec1 & 0x2) )
{
xR = eMD0[4] = this->Dims[0] - 1; //reset right trim
}
}
}
else //contour cuts through without intersecting x-edges, reset trim edges
{
xL = eMD0[3] = 0;
xR = eMD0[4] = this->Dims[0]-1;
}
// Okay run along the x-pixels and count the number of
// y-intersections. Here we are just checking y edges that make up the
......@@ -568,7 +588,7 @@ GenerateOutput(double value, T* rowPtr, vtkIdType row)
T *sPtr;
float x[3];
x[1] = this->Origin[this->Axis1] + row*this->Spacing[this->Axis1];
x[2] = this->X;
x[2] = this->Z;
for (i=xL; i < xR; ++i)
{
if ( (numLines=this->GetNumberOfPrimitives(eCase)) > 0 )
......@@ -622,8 +642,9 @@ ContourImage(vtkFlyingEdges2D *self, T *scalars, vtkPoints *newPts,
// Figure out which 2D plane the image lies in. Capture information for
// subsequent processing.
vtkFlyingEdges2DAlgorithm<T> algo;
algo.Origin = input->GetOrigin();
algo.Spacing = input->GetSpacing();
input->GetOrigin(algo.Origin);
input->GetSpacing(algo.Spacing);
algo.AdjustOrigin(updateExt);
if (updateExt[4] == updateExt[5])
{ // z collapsed
algo.Axis0 = 0;
......@@ -634,7 +655,7 @@ ContourImage(vtkFlyingEdges2D *self, T *scalars, vtkPoints *newPts,
algo.Min1 = updateExt[2];
algo.Max1 = updateExt[3];
algo.Inc1 = incs[1];
algo.X = algo.Origin[2] + (updateExt[4]*algo.Spacing[2]);
algo.Z = algo.Origin[2] + (updateExt[4]*algo.Spacing[2]);
algo.Axis2 = 2;
}
else if (updateExt[2] == updateExt[3])
......@@ -647,7 +668,7 @@ ContourImage(vtkFlyingEdges2D *self, T *scalars, vtkPoints *newPts,
algo.Min1 = updateExt[4];
algo.Max1 = updateExt[5];
algo.Inc1 = incs[2];
algo.X = algo.Origin[1] + (updateExt[2]*algo.Spacing[1]);
algo.Z = algo.Origin[1] + (updateExt[2]*algo.Spacing[1]);
algo.Axis2 = 1;
}
else if (updateExt[0] == updateExt[1])
......@@ -660,7 +681,7 @@ ContourImage(vtkFlyingEdges2D *self, T *scalars, vtkPoints *newPts,
algo.Min1 = updateExt[4];
algo.Max1 = updateExt[5];
algo.Inc1 = incs[2];
algo.X = algo.Origin[0] + (updateExt[0]*algo.Spacing[0]);
algo.Z = algo.Origin[0] + (updateExt[0]*algo.Spacing[0]);
algo.Axis2 = 0;
}
else
......@@ -735,21 +756,25 @@ ContourImage(vtkFlyingEdges2D *self, T *scalars, vtkPoints *newPts,
}
// Output can now be allocated.
newPts->GetData()->WriteVoidPointer(0,3*(numOutXPts+numOutYPts));
algo.NewPoints = static_cast<float*>(newPts->GetVoidPointer(0));
newLines->WritePointer(numOutLines,3*numOutLines);
algo.NewLines = static_cast<vtkIdType*>(newLines->GetPointer());
if (newScalars)
vtkIdType totalPts = numOutXPts + numOutYPts;
if ( totalPts > 0 )
{
newScalars->WriteVoidPointer(0,numOutXPts+numOutYPts);
algo.NewScalars = static_cast<T*>(newScalars->GetVoidPointer(0));
T TValue = static_cast<T>(value);
std::fill_n(algo.NewScalars, numOutXPts+numOutYPts, TValue);
}
newPts->GetData()->WriteVoidPointer(0,3*totalPts);
algo.NewPoints = static_cast<float*>(newPts->GetVoidPointer(0));
newLines->WritePointer(numOutLines,3*numOutLines);
algo.NewLines = static_cast<vtkIdType*>(newLines->GetPointer());
if (newScalars)
{
newScalars->WriteVoidPointer(0,numOutXPts+numOutYPts);
algo.NewScalars = static_cast<T*>(newScalars->GetVoidPointer(0));
T TValue = static_cast<T>(value);
std::fill_n(algo.NewScalars, totalPts, TValue);
}
// PASS 4: Now process each x-row and produce the output primitives.
Pass4<T> pass4(&algo,value);
vtkSMPTools::For(0,algo.Dims[1]-1, pass4);
// PASS 4: Now process each x-row and produce the output primitives.
Pass4<T> pass4(&algo,value);
vtkSMPTools::For(0,algo.Dims[1]-1, pass4);
}//if output generated
// Handle multiple contours
startXPts = numOutXPts;
......
......@@ -42,7 +42,7 @@ class vtkFlyingEdges3DAlgorithm
{
public:
// Edge case table values.
enum {
enum EdgeClass {
Below = 0, //below isovalue
Above = 1, //above isovalue
LeftAbove = 1, //left vertex is above isovalue
......@@ -51,11 +51,11 @@ public:
};
// Dealing with boundary situations when processing volumes.
enum {
enum CellClass {
Interior = 0,
MinBoundary = 1,
MaxBoundary = 2
} vtkBoundarySituations;
};
// Edge-based case table to generate output triangle primitives. It is
// equivalent to the vertex-based Marching Cubes case table but provides
......@@ -100,8 +100,8 @@ public:
// image data in a form more convenient to the algorithm.
T *Scalars;
vtkIdType Dims[3];
double *Origin;
double *Spacing;
double Origin[3];
double Spacing[3];
vtkIdType NumberOfEdges;
vtkIdType SliceOffset;
int Min0;
......@@ -125,6 +125,14 @@ public:
// Setup algorithm
vtkFlyingEdges3DAlgorithm();
// Adjust the origin to the lower-left corner of the volume (if necessary)
void AdjustOrigin()
{
this->Origin[0] = this->Origin[0] + this->Spacing[0]*this->Min0;
this->Origin[1] = this->Origin[1] + this->Spacing[1]*this->Min1;
this->Origin[2] = this->Origin[2] + this->Spacing[2]*this->Min2;;
}
// The three main passes of the algorithm.
void ProcessXEdge(double value, T const * const inPtr, vtkIdType row, vtkIdType slice); //PASS 1
void ProcessYZEdges(vtkIdType row, vtkIdType slice); //PASS 2
......@@ -844,7 +852,7 @@ template <class T> void vtkFlyingEdges3DAlgorithm<T>::
ProcessYZEdges(vtkIdType row, vtkIdType slice)
{
// Grab the four edge cases bounding this voxel x-row.
unsigned char *ePtr[4], ec0, ec1, ec2, ec3;
unsigned char *ePtr[4], ec0, ec1, ec2, ec3, xInts=1;
ePtr[0] = this->XCases + slice*this->SliceOffset + row*(this->Dims[0]-1);
ePtr[1] = ePtr[0] + this->Dims[0]-1;
ePtr[2] = ePtr[0] + this->SliceOffset;
......@@ -867,6 +875,10 @@ ProcessYZEdges(vtkIdType row, vtkIdType slice)
{
return; //there are no y- or z-ints, thus no contour, skip voxel row
}
else
{
xInts = 0; //there are y- or z- edge ints however
}
}
// Determine proximity to the boundary of volume. This information is used
......@@ -883,33 +895,41 @@ ProcessYZEdges(vtkIdType row, vtkIdType slice)
// row trim edges, need to check all four x-edges.
vtkIdType xL=eMD[0][4], xR=eMD[0][5];
vtkIdType i;
for (i=1; i < 4; ++i)
if ( xInts )
{
xL = ( eMD[i][4] < xL ? eMD[i][4] : xL);
xR = ( eMD[i][5] > xR ? eMD[i][5] : xR);
}
for (i=1; i < 4; ++i)
{
xL = ( eMD[i][4] < xL ? eMD[i][4] : xL);
xR = ( eMD[i][5] > xR ? eMD[i][5] : xR);
}
if ( xL > 0 ) //if trimmed in the -x direction
{
ec0 = *(ePtr[0]+xL); ec1 = *(ePtr[1]+xL);
ec2 = *(ePtr[2]+xL); ec3 = *(ePtr[3]+xL);
if ( (ec0 & 0x1) != (ec1 & 0x1) || (ec1 & 0x1) != (ec2 & 0x1) ||
(ec2 & 0x1) != (ec3 & 0x1) )
if ( xL > 0 ) //if trimmed in the -x direction
{
xL = eMD[0][4] = 0; //reset left trim
ec0 = *(ePtr[0]+xL); ec1 = *(ePtr[1]+xL);
ec2 = *(ePtr[2]+xL); ec3 = *(ePtr[3]+xL);
if ( (ec0 & 0x1) != (ec1 & 0x1) || (ec1 & 0x1) != (ec2 & 0x1) ||
(ec2 & 0x1) != (ec3 & 0x1) )
{
xL = eMD[0][4] = 0; //reset left trim
}
}
}
if ( xR < (this->Dims[0]-1) ) //if trimmed in the +x direction
{
ec0 = *(ePtr[0]+xR); ec1 = *(ePtr[1]+xR);
ec2 = *(ePtr[2]+xR); ec3 = *(ePtr[3]+xR);
if ( (ec0 & 0x2) != (ec1 & 0x2) || (ec1 & 0x2) != (ec2 & 0x2) ||
(ec2 & 0x2) != (ec3 & 0x2) )
if ( xR < (this->Dims[0]-1) ) //if trimmed in the +x direction
{
xR = eMD[0][5] = this->Dims[0]-1; //reset right trim
ec0 = *(ePtr[0]+xR); ec1 = *(ePtr[1]+xR);
ec2 = *(ePtr[2]+xR); ec3 = *(ePtr[3]+xR);
if ( (ec0 & 0x2) != (ec1 & 0x2) || (ec1 & 0x2) != (ec2 & 0x2) ||
(ec2 & 0x2) != (ec3 & 0x2) )
{
xR = eMD[0][5] = this->Dims[0]-1; //reset right trim
}
}
}
else //contour cuts through without intersecting x-edges, reset trim edges
{
xL = eMD[0][4] = 0;
xR = eMD[0][5] = this->Dims[0]-1;
}
// Okay run along the x-voxels and count the number of y- and
// z-intersections. Here we are just checking y,z edges that make up the
......@@ -1065,8 +1085,8 @@ Contour(vtkFlyingEdges3D *self, vtkImageData *input, int extent[6],
// subsequent processing.
vtkFlyingEdges3DAlgorithm<T> algo;
algo.Scalars = scalars;
algo.Origin = input->GetOrigin();
algo.Spacing = input->GetSpacing();
input->GetOrigin(algo.Origin);
input->GetSpacing(algo.Spacing);
algo.Min0 = extent[0];
algo.Max0 = extent[1];
algo.Inc0 = incs[0];
......@@ -1076,6 +1096,7 @@ Contour(vtkFlyingEdges3D *self, vtkImageData *input, int extent[6],
algo.Min2 = extent[4];
algo.Max2 = extent[5];
algo.Inc2 = incs[2];
algo.AdjustOrigin();
// Now allocate working arrays. The XCases array tracks x-edge cases.
algo.Dims[0] = algo.Max0 - algo.Min0 + 1;
......@@ -1148,35 +1169,41 @@ Contour(vtkFlyingEdges3D *self, vtkImageData *input, int extent[6],
}
// Output can now be allocated.
newPts->GetData()->WriteVoidPointer(0,3*(numOutXPts+numOutYPts+numOutZPts));
algo.NewPoints = static_cast<float*>(newPts->GetVoidPointer(0));
newTris->WritePointer(numOutTris,4*numOutTris);
algo.NewTris = static_cast<vtkIdType*>(newTris->GetPointer());
if (newScalars)
{
newScalars->WriteVoidPointer(0,(numOutXPts+numOutYPts+numOutZPts));
algo.NewScalars = static_cast<T*>(newScalars->GetVoidPointer(0));
T TValue = static_cast<T>(value);
std::fill_n(algo.NewScalars, numOutXPts+numOutYPts+numOutZPts, TValue);
}
if (newGradients)
{
newGradients->WriteVoidPointer(0,3*(numOutXPts+numOutYPts+numOutZPts));
algo.NewGradients = static_cast<float*>(newGradients->GetVoidPointer(0));
}
if (newNormals)
vtkIdType totalPts = numOutXPts + numOutYPts + numOutZPts;
if ( totalPts > 0 )
{
newNormals->WriteVoidPointer(0,3*(numOutXPts+numOutYPts+numOutZPts));
algo.NewNormals = static_cast<float*>(newNormals->GetVoidPointer(0));
}
algo.NeedGradients = (algo.NewGradients || algo.NewNormals ? 1 : 0);
// PASS 4: Fourth and final pass: Process voxel rows and generate output.
// Note that we are simultaneously generating triangles and interpolating
// points. These could be split into separate, parallel operations for
// maximum performance.
Pass4<T> pass4(&algo,value);
vtkSMPTools::For(0,algo.Dims[2]-1, pass4);
newPts->GetData()->WriteVoidPointer(0,3*totalPts);
algo.NewPoints = static_cast<float*>(newPts->GetVoidPointer(0));
newTris->WritePointer(numOutTris,4*numOutTris);
algo.NewTris = static_cast<vtkIdType*>(newTris->GetPointer());
if (newScalars)
{
vtkIdType numPrevPts = newScalars->GetNumberOfTuples();
vtkIdType numNewPts = totalPts - numPrevPts;
newScalars->WriteVoidPointer(0,totalPts);
algo.NewScalars = static_cast<T*>(newScalars->GetVoidPointer(0));
T TValue = static_cast<T>(value);
std::fill_n(algo.NewScalars+numPrevPts, numNewPts, TValue);
}
if (newGradients)
{
newGradients->WriteVoidPointer(0,3*totalPts);
algo.NewGradients = static_cast<float*>(newGradients->GetVoidPointer(0));
}
if (newNormals)
{
newNormals->WriteVoidPointer(0,3*totalPts);
algo.NewNormals = static_cast<float*>(newNormals->GetVoidPointer(0));
}
algo.NeedGradients = (algo.NewGradients || algo.NewNormals ? 1 : 0);
// PASS 4: Fourth and final pass: Process voxel rows and generate output.
// Note that we are simultaneously generating triangles and interpolating
// points. These could be split into separate, parallel operations for
// maximum performance.
Pass4<T> pass4(&algo,value);
vtkSMPTools::For(0,algo.Dims[2]-1, pass4);
}//if anything generated
// Handle multiple contours
startXPts = numOutXPts;
......