Commit 8a2dc609 authored by Steven Hahn's avatar Steven Hahn

Minimize integer division in surface and triangle filters.

Recalculating the cell indices from the cell id is slow
and in these particular cases unnecessary.
parent ef07338e
......@@ -110,6 +110,11 @@ public:
* THIS METHOD IS NOT THREAD SAFE.
*/
virtual vtkCell *GetCell(vtkIdType cellId) = 0;
virtual vtkCell *GetCell(int vtkNotUsed(i), int vtkNotUsed(j), int vtkNotUsed(k))
{
vtkErrorMacro("ijk indices are only valid with structured data!");
return NULL;
}
/**
* Get cell with cellId such that: 0 <= cellId < NumberOfCells.
......
......@@ -315,6 +315,7 @@ public:
* Get cell with cellId such that: 0 <= cellId < NumberOfCells.
* THIS METHOD IS NOT THREAD SAFE.
*/
using vtkDataSet::GetCell;
vtkCell *GetCell(vtkIdType cellId) VTK_OVERRIDE;
/**
......
......@@ -244,6 +244,7 @@ public:
* Get cell with cellId such that: 0 <= cellId < NumberOfCells.
* THIS METHOD IS NOT THREAD SAFE.
*/
using vtkDataSet::GetCell;
vtkCell* GetCell( vtkIdType ) VTK_OVERRIDE;
/**
......
......@@ -276,6 +276,119 @@ vtkCell *vtkImageData::GetCell(vtkIdType cellId)
return cell;
}
vtkCell *vtkImageData::GetCell(int iMin, int jMin, int kMin) {
vtkCell *cell = NULL;
int loc[3];
vtkIdType idx, npts;
int iMax = 0, jMax = 0, kMax = 0;
double x[3];
const double *origin = this->Origin;
const double *spacing = this->Spacing;
const int *extent = this->Extent;
// Use vtkIdType to avoid overflow on large images
vtkIdType cellDims[3];
cellDims[0] = extent[1] - extent[0];
cellDims[1] = extent[3] - extent[2];
cellDims[2] = extent[5] - extent[4];
vtkIdType dims[3];
dims[0] = cellDims[0] + 1;
dims[1] = cellDims[1] + 1;
dims[2] = cellDims[2] + 1;
vtkIdType d01 = dims[0] * dims[1];
if (dims[0] == 0 || dims[1] == 0 || dims[2] == 0) {
vtkErrorMacro("Requesting a cell from an empty image.");
return NULL;
}
switch (this->DataDescription) {
case VTK_EMPTY:
// cell = this->EmptyCell;
return NULL;
case VTK_SINGLE_POINT: // cellId can only be = 0
cell = this->Vertex;
break;
case VTK_X_LINE:
iMax = iMin + 1;
jMax = jMin = 0;
kMax = kMin = 0;
cell = this->Line;
break;
case VTK_Y_LINE:
iMax = iMin = 0;
jMax = jMin + 1;
kMax = kMin = 0;
cell = this->Line;
break;
case VTK_Z_LINE:
iMax = iMin = 0;
jMax = jMin = 0;
kMax = kMin + 1;
cell = this->Line;
break;
case VTK_XY_PLANE:
iMax = iMin + 1;
jMax = jMin + 1;
kMax = kMin = 0;
cell = this->Pixel;
break;
case VTK_YZ_PLANE:
iMax = iMin = 0;
jMax = jMin + 1;
kMax = kMin + 1;
cell = this->Pixel;
break;
case VTK_XZ_PLANE:
iMax = iMin + 1;
jMax = jMin = 0;
kMax = kMin + 1;
cell = this->Pixel;
break;
case VTK_XYZ_GRID:
iMax = iMin + 1;
jMax = jMin + 1;
kMax = kMin + 1;
cell = this->Voxel;
break;
default:
vtkErrorMacro("Invalid DataDescription.");
return NULL;
}
// Extract point coordinates and point ids
// Ids are relative to extent min.
npts = 0;
for (loc[2] = kMin; loc[2] <= kMax; loc[2]++)
{
x[2] = origin[2] + (loc[2] + extent[4]) * spacing[2];
for (loc[1] = jMin; loc[1] <= jMax; loc[1]++)
{
x[1] = origin[1] + (loc[1] + extent[2]) * spacing[1];
for (loc[0] = iMin; loc[0] <= iMax; loc[0]++)
{
x[0] = origin[0] + (loc[0] + extent[0]) * spacing[0];
idx = loc[0] + loc[1] * dims[0] + loc[2] * d01;
cell->PointIds->SetId(npts, idx);
cell->Points->SetPoint(npts++, x);
}
}
}
return cell;
}
//----------------------------------------------------------------------------
void vtkImageData::GetCell(vtkIdType cellId, vtkGenericCell *cell)
{
......
......@@ -64,6 +64,7 @@ public:
double *GetPoint(vtkIdType ptId) VTK_OVERRIDE;
void GetPoint(vtkIdType id, double x[3]) VTK_OVERRIDE;
vtkCell *GetCell(vtkIdType cellId) VTK_OVERRIDE;
vtkCell *GetCell(int i, int j, int k) VTK_OVERRIDE;
void GetCell(vtkIdType cellId, vtkGenericCell *cell) VTK_OVERRIDE;
void GetCellBounds(vtkIdType cellId, double bounds[6]) VTK_OVERRIDE;
virtual vtkIdType FindPoint(double x, double y, double z)
......
......@@ -161,6 +161,7 @@ public:
void CopyStructure(vtkDataSet *pd) VTK_OVERRIDE;
void ShallowCopy(vtkDataObject *src) VTK_OVERRIDE;
vtkIdType GetNumberOfCells() VTK_OVERRIDE;
using vtkDataSet::GetCell;
vtkCell* GetCell(vtkIdType cellId) VTK_OVERRIDE;
void GetCell(vtkIdType cellId, vtkGenericCell *cell) VTK_OVERRIDE;
int GetCellType(vtkIdType cellId) VTK_OVERRIDE;
......
......@@ -85,6 +85,7 @@ public:
* vtkPath doesn't use cells. These methods return trivial values.
*/
vtkIdType GetNumberOfCells() VTK_OVERRIDE { return 0; }
using vtkDataSet::GetCell;
vtkCell *GetCell(vtkIdType) VTK_OVERRIDE { return NULL; }
void GetCell(vtkIdType, vtkGenericCell *) VTK_OVERRIDE;
int GetCellType(vtkIdType) VTK_OVERRIDE { return 0; }
......
......@@ -99,6 +99,7 @@ public:
* Standard vtkDataSet interface.
*/
vtkIdType GetNumberOfCells() VTK_OVERRIDE;
using vtkDataSet::GetCell;
vtkCell *GetCell(vtkIdType cellId) VTK_OVERRIDE;
void GetCell(vtkIdType cellId, vtkGenericCell *cell) VTK_OVERRIDE;
int GetCellType(vtkIdType cellId) VTK_OVERRIDE;
......
......@@ -241,6 +241,101 @@ vtkCell *vtkRectilinearGrid::GetCell(vtkIdType cellId)
return cell;
}
//----------------------------------------------------------------------------
vtkCell *vtkRectilinearGrid::GetCell(int iMin, int jMin, int kMin) {
vtkCell *cell = NULL;
vtkIdType idx, npts;
int loc[3];
int iMax, jMax, kMax;
int d01 = this->Dimensions[0] * this->Dimensions[1];
double x[3];
iMin = iMax = jMin = jMax = kMin = kMax = 0;
switch (this->DataDescription) {
case VTK_EMPTY:
// return this->EmptyCell;
return NULL;
case VTK_SINGLE_POINT: // cellId can only be = 0
cell = this->Vertex;
break;
case VTK_X_LINE:
iMax = iMin + 1;
jMin = jMax = 0;
kMin = kMax = 0;
cell = this->Line;
break;
case VTK_Y_LINE:
iMin = iMax = 0;
jMax = jMin + 1;
kMin = kMax = 0;
cell = this->Line;
break;
case VTK_Z_LINE:
iMin = iMax = 0;
jMin = jMax = 0;
kMax = kMin + 1;
cell = this->Line;
break;
case VTK_XY_PLANE:
iMax = iMin + 1;
jMax = jMin + 1;
kMin = kMax = 0;
cell = this->Pixel;
break;
case VTK_YZ_PLANE:
iMin = iMax = 0;
jMax = jMin + 1;
kMax = kMin + 1;
cell = this->Pixel;
break;
case VTK_XZ_PLANE:
iMax = iMin + 1;
jMin = kMax = 0;
kMax = kMin + 1;
cell = this->Pixel;
break;
case VTK_XYZ_GRID:
iMax = iMin + 1;
jMax = jMin + 1;
kMax = kMin + 1;
cell = this->Voxel;
break;
default:
vtkErrorMacro(<< "Invalid DataDescription.");
return NULL;
}
// Extract point coordinates and point ids
for (npts = 0, loc[2] = kMin; loc[2] <= kMax; loc[2]++)
{
x[2] = this->ZCoordinates->GetComponent(loc[2], 0);
for (loc[1] = jMin; loc[1] <= jMax; loc[1]++)
{
x[1] = this->YCoordinates->GetComponent(loc[1], 0);
for (loc[0] = iMin; loc[0] <= iMax; loc[0]++)
{
x[0] = this->XCoordinates->GetComponent(loc[0], 0);
idx = loc[0] + loc[1] * this->Dimensions[0] + loc[2] * d01;
cell->PointIds->SetId(npts, idx);
cell->Points->SetPoint(npts++, x);
}
}
}
return cell;
}
//----------------------------------------------------------------------------
void vtkRectilinearGrid::GetCell(vtkIdType cellId, vtkGenericCell *cell)
{
......
......@@ -80,6 +80,7 @@ public:
double *GetPoint(vtkIdType ptId) VTK_OVERRIDE;
void GetPoint(vtkIdType id, double x[3]) VTK_OVERRIDE;
vtkCell *GetCell(vtkIdType cellId) VTK_OVERRIDE;
vtkCell *GetCell(int i, int j, int k) VTK_OVERRIDE;
void GetCell(vtkIdType cellId, vtkGenericCell *cell) VTK_OVERRIDE;
void GetCellBounds(vtkIdType cellId, double bounds[6]) VTK_OVERRIDE;
vtkIdType FindPoint(double x, double y, double z) { return this->vtkDataSet::FindPoint(x, y, z);};
......
......@@ -269,7 +269,128 @@ vtkCell *vtkStructuredGrid::GetCell(vtkIdType cellId)
idx = cell->PointIds->GetId(i);
cell->Points->SetPoint(i,this->Points->GetPoint(idx));
}
return cell;
}
//----------------------------------------------------------------------------
vtkCell *vtkStructuredGrid::GetCell(int i, int j, int k) {
vtkIdType cellId =
i + (j + (k * (this->Dimensions[1] - 1))) * (this->Dimensions[0] - 1);
vtkCell *cell = NULL;
vtkIdType idx;
int d01, offset1, offset2;
// Make sure data is defined
if (!this->Points)
{
vtkErrorMacro(<< "No data");
return NULL;
}
// see whether the cell is blanked
if (!this->IsCellVisible(cellId))
{
return this->EmptyCell;
}
// Update dimensions
this->GetDimensions();
switch (this->DataDescription)
{
case VTK_EMPTY:
return this->EmptyCell;
case VTK_SINGLE_POINT: // cellId can only be = 0
cell = this->Vertex;
cell->PointIds->SetId(0, 0);
break;
case VTK_X_LINE:
cell = this->Line;
cell->PointIds->SetId(0, cellId);
cell->PointIds->SetId(1, cellId + 1);
break;
case VTK_Y_LINE:
cell = this->Line;
cell->PointIds->SetId(0, cellId);
cell->PointIds->SetId(1, cellId + 1);
break;
case VTK_Z_LINE:
cell = this->Line;
cell->PointIds->SetId(0, cellId);
cell->PointIds->SetId(1, cellId + 1);
break;
case VTK_XY_PLANE:
cell = this->Quad;
idx = i + j * this->Dimensions[0];
offset1 = 1;
offset2 = this->Dimensions[0];
cell->PointIds->SetId(0, idx);
cell->PointIds->SetId(1, idx + offset1);
cell->PointIds->SetId(2, idx + offset1 + offset2);
cell->PointIds->SetId(3, idx + offset2);
break;
case VTK_YZ_PLANE:
cell = this->Quad;
idx = j + k * this->Dimensions[1];
offset1 = 1;
offset2 = this->Dimensions[1];
cell->PointIds->SetId(0, idx);
cell->PointIds->SetId(1, idx + offset1);
cell->PointIds->SetId(2, idx + offset1 + offset2);
cell->PointIds->SetId(3, idx + offset2);
break;
case VTK_XZ_PLANE:
cell = this->Quad;
idx = i + k * this->Dimensions[0];
offset1 = 1;
offset2 = this->Dimensions[0];
cell->PointIds->SetId(0, idx);
cell->PointIds->SetId(1, idx + offset1);
cell->PointIds->SetId(2, idx + offset1 + offset2);
cell->PointIds->SetId(3, idx + offset2);
break;
case VTK_XYZ_GRID:
cell = this->Hexahedron;
d01 = this->Dimensions[0] * this->Dimensions[1];
idx = i + j * this->Dimensions[0] + k * d01;
offset1 = 1;
offset2 = this->Dimensions[0];
cell->PointIds->SetId(0, idx);
cell->PointIds->SetId(1, idx + offset1);
cell->PointIds->SetId(2, idx + offset1 + offset2);
cell->PointIds->SetId(3, idx + offset2);
idx += d01;
cell->PointIds->SetId(4, idx);
cell->PointIds->SetId(5, idx + offset1);
cell->PointIds->SetId(6, idx + offset1 + offset2);
cell->PointIds->SetId(7, idx + offset2);
break;
default:
vtkErrorMacro(<< "Invalid DataDescription.");
return NULL;
}
// Extract point coordinates and point ids. NOTE: the ordering of the vtkQuad
// and vtkHexahedron cells are tricky.
int NumberOfIds = cell->PointIds->GetNumberOfIds();
for (i = 0; i < NumberOfIds; i++)
{
idx = cell->PointIds->GetId(i);
cell->Points->SetPoint(i, this->Points->GetPoint(idx));
}
return cell;
}
......
......@@ -83,6 +83,7 @@ public:
void GetPoint(vtkIdType ptId, double p[3]) VTK_OVERRIDE
{this->vtkPointSet::GetPoint(ptId,p);}
vtkCell *GetCell(vtkIdType cellId) VTK_OVERRIDE;
vtkCell *GetCell(int i, int j, int k) VTK_OVERRIDE;
void GetCell(vtkIdType cellId, vtkGenericCell *cell) VTK_OVERRIDE;
void GetCellBounds(vtkIdType cellId, double bounds[6]) VTK_OVERRIDE;
int GetCellType(vtkIdType cellId) VTK_OVERRIDE;
......
......@@ -313,6 +313,125 @@ vtkCell *vtkUniformGrid::GetCell(vtkIdType cellId)
return cell;
}
//----------------------------------------------------------------------------
vtkCell *vtkUniformGrid::GetCell(int iMin, int jMin, int kMin)
{
vtkIdType cellId = iMin + (jMin + (kMin * (this->Dimensions[1] - 1))) * (this->Dimensions[0] - 1);
vtkCell *cell = NULL;
int loc[3];
vtkIdType idx, npts;
int iMax = 0, jMax = 0, kMax = 0;
double x[3];
double *origin = this->GetOrigin();
double *spacing = this->GetSpacing();
int extent[6];
this->GetExtent(extent);
int dims[3];
dims[0] = extent[1] - extent[0] + 1;
dims[1] = extent[3] - extent[2] + 1;
dims[2] = extent[5] - extent[4] + 1;
int d01 = dims[0]*dims[1];
if (dims[0] == 0 || dims[1] == 0 || dims[2] == 0)
{
vtkErrorMacro("Requesting a cell from an empty image.");
return this->GetEmptyCell();
}
// see whether the cell is blanked
if (!this->IsCellVisible(cellId) )
{
return this->GetEmptyCell();
}
switch (this->GetDataDescription())
{
case VTK_EMPTY:
return this->GetEmptyCell();
case VTK_SINGLE_POINT: // cellId can only be = 0
cell = this->Vertex;
break;
case VTK_X_LINE:
iMax = iMin + 1;
jMax = jMin = 0;
kMax = kMin = 0;
cell = this->Line;
break;
case VTK_Y_LINE:
iMax = iMin = 0;
jMax = jMin + 1;
kMax = kMin = 0;
cell = this->Line;
break;
case VTK_Z_LINE:
iMin = iMax = 0;
jMin = jMax = 0;
kMax = kMin + 1;
cell = this->Line;
break;
case VTK_XY_PLANE:
iMax = iMin + 1;
jMax = jMin + 1;
kMin = kMax = 0;
cell = this->Pixel;
break;
case VTK_YZ_PLANE:
iMin = iMax = 0;
jMax = jMin + 1;
kMax = kMin + 1;
cell = this->Pixel;
break;
case VTK_XZ_PLANE:
iMax = iMin + 1;
jMin = kMax = 0;
kMax = kMin + 1;
cell = this->Pixel;
break;
case VTK_XYZ_GRID:
iMax = iMin + 1;
jMax = jMin + 1;
kMax = kMin + 1;
cell = this->Voxel;
break;
default:
vtkErrorMacro(<<"Invalid DataDescription.");
return NULL;
}
// Extract point coordinates and point ids
// Ids are relative to extent min.
npts = 0;
for (loc[2]=kMin; loc[2]<=kMax; loc[2]++)
{
x[2] = origin[2] + (loc[2]+extent[4]) * spacing[2];
for (loc[1]=jMin; loc[1]<=jMax; loc[1]++)
{
x[1] = origin[1] + (loc[1]+extent[2]) * spacing[1];
for (loc[0]=iMin; loc[0]<=iMax; loc[0]++)
{
x[0] = origin[0] + (loc[0]+extent[0]) * spacing[0];
idx = loc[0] + loc[1]*dims[0] + loc[2]*d01;
cell->PointIds->SetId(npts,idx);
cell->Points->SetPoint(npts++,x);
}
}
}
return cell;
}
//----------------------------------------------------------------------------
void vtkUniformGrid::GetCell(vtkIdType cellId, vtkGenericCell *cell)
{
......
......@@ -58,6 +58,7 @@ public:
/**
* Standard vtkDataSet API methods. See vtkDataSet for more information.
*/
vtkCell *GetCell(int i, int j, int k) VTK_OVERRIDE;
vtkCell *GetCell(vtkIdType cellId) VTK_OVERRIDE;
void GetCell(vtkIdType cellId, vtkGenericCell *cell) VTK_OVERRIDE;
vtkIdType FindCell(
......
......@@ -134,6 +134,7 @@ public:
void Reset();
void CopyStructure(vtkDataSet *ds) VTK_OVERRIDE;
vtkIdType GetNumberOfCells() VTK_OVERRIDE;
using vtkDataSet::GetCell;
vtkCell *GetCell(vtkIdType cellId) VTK_OVERRIDE;
void GetCell(vtkIdType cellId, vtkGenericCell *cell) VTK_OVERRIDE;
void GetCellBounds(vtkIdType cellId, double bounds[6]) VTK_OVERRIDE;
......
......@@ -146,7 +146,7 @@ void vtkDataSetTriangleFilter::StructuredExecute(vtkDataSet *input,
for (i = 0; i < dimensions[0]; i++)
{
inId = i+(j+(k*dimensions[1]))*dimensions[0];
vtkCell *cell = input->GetCell(inId);
vtkCell *cell = input->GetCell(i, j, k);
if ((i+j+k)%2 == 0)
{
cell->Triangulate(0, cellPtIds, cellPts);
......
......@@ -203,7 +203,7 @@ int vtkDataSetSurfaceFilter::RequestData(
vtkStructuredGrid *grid = vtkStructuredGrid::SafeDownCast(input);
if (grid->HasAnyBlankCells())
{
return this->DataSetExecute(grid, output);
return this->StructuredWithBlankingExecute(grid, output);
}
else
{
......@@ -1038,6 +1038,172 @@ void vtkDataSetSurfaceFilter::ExecuteFaceQuads(vtkDataSet *input,
}
}
//----------------------------------------------------------------------------
int vtkDataSetSurfaceFilter::StructuredWithBlankingExecute(vtkStructuredGrid *input,
vtkPolyData *output)
{
vtkIdType newCellId;
vtkIdType numPts=input->GetNumberOfPoints();
vtkIdType numCells=input->GetNumberOfCells();
vtkCell *face;
double x[3];
vtkIdList *cellIds;
vtkIdList *pts;
vtkPoints *newPts;
vtkIdType ptId, pt;
int npts;
vtkPointData *pd = input->GetPointData();
vtkCellData *cd = input->GetCellData();
vtkPointData *outputPD = output->GetPointData();
vtkCellData *outputCD = output->GetCellData();
if (numCells == 0)
{
vtkWarningMacro(<<"Number of cells is zero, no data to process.");
return 1;
}
if (this->PassThroughCellIds)
{
this->OriginalCellIds = vtkIdTypeArray::New();
this->OriginalCellIds->SetName(this->GetOriginalCellIdsName());
this->OriginalCellIds->SetNumberOfComponents(1);
this->OriginalCellIds->Allocate(numCells);
outputCD->AddArray(this->OriginalCellIds);
}
if (this->PassThroughPointIds)
{
this->OriginalPointIds = vtkIdTypeArray::New();
this->OriginalPointIds->SetName(this->GetOriginalPointIdsName());
this->OriginalPointIds->SetNumberOfComponents(1);
this->OriginalPointIds->Allocate(numPts);
outputPD->AddArray(this->OriginalPointIds);
}
cellIds = vtkIdList::New();
pts = vtkIdList::New();
vtkDebugMacro(<<"Executing geometry filter");
// Allocate
//
newPts = vtkPoints::New();
// we don't know what type of data the input points are so
// we keep the output points to have the default type (float)
newPts->Allocate(numPts,numPts/2);
output->Allocate(4*numCells,numCells/2);
outputPD->CopyGlobalIdsOn();
outputPD->CopyAllocate(pd,numPts,numPts/2);
outputCD->CopyGlobalIdsOn();
outputCD->CopyAllocate(cd,numCells,numCells/2);
// Traverse cells to extract geometry
//
int abort=0;
int dims[3];
input->GetCellDims(dims);
vtkIdType d01 = static_cast<vtkIdType>(dims[0])*dims[1];
for (int k = 0; k < dims[2] && !abort; ++k)
{
vtkDebugMacro(<< "Process cell #" << d01*k);
this->UpdateProgress(k / dims[2]);
abort = this->GetAbortExecute();
for (int j = 0; j < dims[1]; ++j)
{
for (int i = 0; i < dims[0]; ++i)
{
vtkIdType cellId = d01*k + dims[0]*j + i;
if (!input->IsCellVisible(cellId))
{
continue;
}
vtkCell *cell = input->GetCell(i,j,k);
switch (cell->GetCellDimension())