[vtkusers] Creating a "upper sum" geometry

[Tobias Meisen]

Hi everybody,

I would like to create some kind of "upper sum" geometry using vtk. The
“upper sum” geometry is created for a given input geometry. It has to be
an image data dataset with a spacing defined by the user. It has two
dataarrays: one cellarray called “cellinside” and one pointarray called
“pointinside”. The pointdata is “1” if the point is inside the input
geometry, otherwise “0”. The celldata is “1” if one point of the cell is
inside the input geometry.
I thought I can use the findCell Method of vtkDataset to check whether a
point is inside the input geometry or not, but it failed in some test
cases declaring some points as outside the geometry that are inside.
Here is the extract of the algorithm I am using (saving and loading of
the input and output geometry is not part of the code extract):

// get bounds
vtkDataSet* inputSet = reader->GetOutput();
double* bounds = inputSet->GetBounds();

// create output
vtkImageData* outputSet = vtkImageData::New();
outputSet->SetOrigin(bounds[0], bounds[2], bounds[4]);
double dim_x = ceil((bounds[1] - bounds[0]) / spacing_x) + 1;
double dim_y = ceil((bounds[3] - bounds[2]) / spacing_y) + 1;
double dim_z = ceil((bounds[5] - bounds[4]) / spacing_z) + 1;
outputSet->SetDimensions(dim_x, dim_y, dim_z);
outputSet->SetSpacing(spacing_x, spacing_y, spacing_z);
outputSet->Update();

// create data arrays
vtkIntArray* cellvalue = vtkIntArray::New();
vtkIntArray* pointvalue = vtkIntArray::New();
cellvalue->Allocate(outputSet->GetNumberOfCells());
cellvalue->SetName("cellinside");
cellvalue->SetNumberOfValues(outputSet->GetNumberOfCells());
cellvalue->SetNumberOfComponents(1);
pointvalue->Allocate(outputSet->GetNumberOfPoints());
pointvalue->SetName("pointinside");
pointvalue->SetNumberOfValues(outputSet->GetNumberOfPoints());
pointvalue->SetNumberOfComponents(1);
for (int cellId = 0; cellId < outputSet->GetNumberOfCells(); cellId++)
	cellvalue->SetValue(cellId, 0);

// declare variables
vtkIdList* cellIds = vtkIdList::New();
int subId;
double pcoords[3];
double* coords;
double* weights = new double[inputSet->GetMaxCellSize()];

// check points of upper sum geometry
for (int pointId = 0; pointId < outputSet->GetNumberOfPoints(); pointId++)
{
	// get point data
	coords = outputSet->GetPoint(pointId);

	// check if point is inside a cell of the input geometry
	int foundCellId = inputSet->FindCell(coords, NULL, 0, 					0.0, subId,
pcoords, weights);
	if (foundCellId > -1)
	{
		pointvalue->SetValue(pointId, 1);
		outputSet->GetPointCells(pointId, cellIds);
		for (int index = 0; index < cellIds->GetNumberOfIds();
				index++)
		cellvalue->SetValue(cellIds->GetId(index), 1);
	}
	else
		pointvalue->SetValue(pointId, 0);
}

// attach data to output set
cellvalue->Squeeze();
pointvalue->Squeeze();
outputSet->GetCellData()->AddArray(cellvalue);
outputSet->GetPointData()->AddArray(pointvalue);

Anyone can give me a hint, what i am doing wrong.

Thx.

-------------- next part --------------
A non-text attachment was scrubbed...
Name: signature.asc
Type: application/pgp-signature
Size: 251 bytes
Desc: OpenPGP digital signature
URL: <http://www.vtk.org/pipermail/vtkusers/attachments/20091112/2b98b216/attachment.pgp>