[vtkusers] help with Contouring oddities using vtkContourFilter
Subha
smahaade at essence.utsi.edu
Wed Apr 9 12:42:13 EDT 2003
Hello all, ( am re-sending this mail bcos the earlier mail that I sent out
did not attach my file properly).
I tried a small experiment using vtkContourFilter ( intended for use
somewhere else).
I used the code from vtkContourFilter::Execute() function in my main program
(literally copied and pasted the code from the vtkContourFilter class in my
program ) to generate an isosurface for a Structured Grid dataset and timed
it. My program flows like this (program files also attached)
method-1 (my program)
=========
myMainProgram.cxx
{
reader->update()
TIME START
code from vtkContourFilter::Execute() with proper inputs, initializations
TIME END
above generated vtkPolyData
polyDataMapper->SetInput(above generated polyData)
renderWindow->render()
}
I also rendered the isosurface in the way it is usually done (I have timed
this method between setting the input to contour filter and the render
command because of the way the vtk pipeline works ).
method-2
========
TIME START
contour->SetInput(structuredGrid)
polydataMapper->SetInput(contour->GetOutput())
actor->setMapper(polyDataMapper)
renderWindow->render()
TIME END
1) Under the same load conditions on a dedicated PC, method-1 was
significantly slower than method-2.
Method-1 takes as much as 10 seconds more than method-2 to render a 7000 cell
structured grid (which translates to that atleast 100 times the number of
machine cycles on a 733 mhz and 512 mb ram linux machine).
2) Most of the time seems to have been taken by the
vtkHexahedron::Contour(...) function ( this function is called from the
vtkContourFilter::Execute() function. It interpolates & creates the triangles
and new points for the isosurface )
I ran the programs in gdb to look through all the classes where the
vtkHexahedron::Contour(...) function leads, but am unable to say "Why" there
is a significant delay is caused when using method-1. I am missing something
very obvious here.
Will any of you experienced folks who know what might be happening, please
make time to tell me ?
I am grateful for any pointers/ thoughts that you can share with me.
Thanks & Regards,
Subha.
//**************************************************************************//
#include "vtkRenderer.h"
#include "vtkRenderWindow.h"
#include "vtkRenderWindowInteractor.h"
#include "vtkActor.h"
#include "vtkStructuredGridReader.h"
#include "vtkPolyDataMapper.h"
#include "../../patented/vtkMarchingContourFilter.h"
#include "vtkLight.h"
#include "vtkRenderWindowInteractor.h"
#include "vtkMergePoints.h"
#include "vtkMarchingCubesCases.h"
#include <sys/times.h>
#include <unistd.h>
#include <iostream.h>
#include <stdio.h>
#include <iostream.h>
#include "SaveImage.h"
void main( int argc , char *argv[] )
{
char a;
int totalPoints, totalCells;
//float value;
ofstream kout("marchIndex.txt", ios::out);
static int count=0;
clock_t start, end, proStart, proEnd;
struct tms tmsstart, tmsend, tmsProStart, tmsProEnd;
int status;
static long clktck = 0;
vtkRenderer *renderer = vtkRenderer::New();
renderer->SetBackground(1,1,1);
vtkRenderWindow *renWin = vtkRenderWindow::New();
renWin->AddRenderer(renderer);
renWin->SetSize(300,300);
vtkRenderWindowInteractor *iren = vtkRenderWindowInteractor::New();
iren->SetRenderWindow(renWin);
vtkLight *light1 = vtkLight::New();
renderer->AddLight(light1);
vtkActor *actor = vtkActor::New();
cout<< "reading from file" <<endl;
vtkStructuredGridReader *reader = vtkStructuredGridReader::New();
reader->SetFileName(argv[1]);
reader->SetScalarsName("sc2");
reader->Update();
vtkStructuredGrid *structGrid = reader->GetOutput();
vtkPointData *pointData = structGrid->GetPointData();
vtkScalars *scalars = pointData->GetScalars();
vtkPolyDataMapper *volMapper = vtkPolyDataMapper::New();
totalCells = structGrid->GetNumberOfCells();
float value[10];
value[0] = 0.0;
value[1] = 0.1;
value[2] = 0.2;
value[3] = 0.3;
value[4] = 0.5;
value[5] = 0.6;
value[6] = 0.7;
value[7] = 0.8;
value[8] = 0.9;
value[9] = 1.0;
/**** CODE FROM vtkContourFilter.cxx::Execute() **********/
static int CASE_MASK[8] = {1,2,4,8,16,32,64,128};
TRIANGLE_CASES *triCase;
EDGE_LIST *edge;
int i, j, index, *vert;
static int vertMap[8] = { 0, 1, 3, 2, 4, 5, 7, 6 };
if((start = times(&tmsstart)) == -1)
cout << "times error " <<endl;
int estimatedSize = (int) pow ((double) totalCells, .75);
estimatedSize = estimatedSize / 1024 * 1024;
if(estimatedSize < 1024 )
estimatedSize = 1024;
vtkCell *cell;
vtkPoints * newPts = vtkPoints::New();
vtkPointLocator *locator = vtkMergePoints::New();
locator->InitPointInsertion(newPts, structGrid->GetBounds(),
estimatedSize);
vtkCellArray * verts = vtkCellArray::New();
verts->Allocate(estimatedSize, estimatedSize);
vtkCellArray * lines = vtkCellArray::New();
lines->Allocate(estimatedSize, estimatedSize);
vtkCellArray * polys = vtkCellArray::New();
polys->Allocate(estimatedSize, estimatedSize);
vtkPolyData *output = vtkPolyData::New();
vtkScalars *cellScalars = vtkScalars::New();
cellScalars->Allocate(VTK_CELL_SIZE);
//vtkPointData inPd is just pointData here
vtkPointData *outPd = output->GetPointData();
vtkCellData *inCd = structGrid->GetCellData();
vtkCellData *outCd = output->GetCellData();
outPd->CopyScalarsOn();
outPd->InterpolateAllocate(pointData, estimatedSize, estimatedSize);
outCd->CopyAllocate(inCd, estimatedSize, estimatedSize);
int numContours = 10;
vtkIdList *cellPts;
//vtkIdList *NcellList = vtkIdList::New();
static int count1=0;
//vtkPolyData *output1 = vtkPolyData::New();
int cellNumber;
for (int k=0; k < 10; k++)
{
for (int cellId=0; cellId < totalCells; cellId++)
{
prevCells = polys->GetNumberOfCells();
cell = structGrid->GetCell(cellId);
cellPts = cell->GetPointIds();
scalars->GetScalars(cellPts,cellScalars);
for ( i=0, index = 0; i < 8; i++)
{
if (cellScalars->GetScalar(vertMap[i]) >= value[k])
{
index |= CASE_MASK[i];
}
}
cell->Contour(value[k], cellScalars, locator,
verts, lines, polys, pointData, outPd,
inCd, cellId, outCd);
currentCells = polys->GetNumberOfCells();
// if(prevCells < currentCells) {
// kout << cellId << " " << index << " " << k << endl;
// count1 ++;
}
}// for all cells
} // for all contour values
// Update ourselves. Because we don't know up front how many verts,
lines,
// polys we've created, take care to reclaim memory.
//
output->SetPoints(newPts);
newPts->Delete();
if (verts->GetNumberOfCells())
{
output->SetVerts(verts);
}
verts->Delete();
if (lines->GetNumberOfCells())
{
output->SetLines(lines);
}
lines->Delete();
if (polys->GetNumberOfCells())
{
output->SetPolys(polys);
}
polys->Delete();
locator->Initialize();//releases leftover memory
output->Squeeze();
/*** Usual code below ********/
volMapper->SetInput(output);
actor->SetMapper(volMapper);
renderer->AddActor(actor);
renWin->Render();
cout << "polys created " << output->GetNumberOfCells() << endl;
if( ( end = times(&tmsend)) == -1 )
cout << "times error tmsend " << endl;
if( 0 == clktck)
if( (clktck = sysconf(_SC_CLK_TCK)) < 0)
cout << "sysconf error " <<endl;
fprintf(stderr, "processing real: %7.2f\n", (end - start) / (double)
clktck);
fprintf(stderr, "processing user: %7.2f\n", (tmsend.tms_utime -
tmsstart.tms_utime) / (double) clktck);
fprintf(stderr, "processing sys: %7.2f\n", (tmsend.tms_stime -
tmsstart.tms_stime) / (double) clktck);
fprintf(stderr, "processing child user: %7.2f\n", (tmsend.tms_cutime -
tmsstart.tms_cutime) / (double) clktck);
fprintf(stderr, "processing child sys: %7.2f\n", (tmsend.tms_cstime -
tmsstart.tms_cstime) / (double) clktck);
//kout.close();
SAVEIMAGE( renWin );
iren->Start();
/** enable user
interaction using mouse **/
//clean up
cout<<" reclaiming memory "<<endl;
renderer->Delete();
renWin->Delete();
iren->Delete();
reader->Delete();
cellScalars->Delete();
locator->Delete();
output->Delete();
// structGrid->Delete();
// contour->Delete();
// geom->Delete();
volMapper->Delete();
// NcellList->Delete();
light1->Delete();
actor->Delete();
}
*******************************
Subha Mahaadevan
------v--
/ __ |_*\\
|_| |_| \\
*******************************
More information about the vtkusers
mailing list