Extracting Triangles from PolyData while keeping geometry and topology

Support
1

Hello,
I’m trying to extract triangles from a PolyData object depending on the direction in which the triangle faces. These PolyData objects are shaped like a slice of pineapple. Round with a hole in the middle and much wider than high. This shape I want to divide into four areas. The top plane, the bottom plane, the inside wall and the outside wall. The triangles from the top and bottom I get through checking if the triangle’s normal vector faces up or downward. The triangles on the side I get through calculating the distance between on of the triangle points and the center of the slice. If this distance is greater than a threshold the triangle is part of the outer wall. If the distance is smaller than the threshold it is part of the inner wall.
When I try to push the resulting triangles back into a polydata structure the resulting meshes don’t look like anything I would expect. It seems as if the triangles just get stacked on one another. Any ideas what I do wrong, or what I can do?

Best regards, Armin

vtkSmartPointer<vtkPolyData> areaDataTop = vtkSmartPointer<vtkPolyData>::New();
vtkSmartPointer<vtkPolyData> areaDataBottom = vtkSmartPointer<vtkPolyData>::New();
vtkSmartPointer<vtkPolyData> areaDataOuter = vtkSmartPointer<vtkPolyData>::New();
vtkSmartPointer<vtkPolyData> areaDataInner = vtkSmartPointer<vtkPolyData>::New();

void divideIntoAreas(vtkSmartPointer<vtkPolyData> slice, double3 sliceMainAxis, double threshold) {	
	vtkSmartPointer<vtkCellArray> trianglesTop = vtkSmartPointer<vtkCellArray>::New();
	vtkSmartPointer<vtkCellArray> trianglesBottom = vtkSmartPointer<vtkCellArray>::New();
	vtkSmartPointer<vtkCellArray> trianglesOuter = vtkSmartPointer<vtkCellArray>::New();
	vtkSmartPointer<vtkCellArray> trianglesInner = vtkSmartPointer<vtkCellArray>::New();

	vtkSmartPointer<vtkPoints> pointsTop = vtkSmartPointer<vtkPoints>::New();
	vtkSmartPointer<vtkPoints> pointsBottom = vtkSmartPointer<vtkPoints>::New();
	vtkSmartPointer<vtkPoints> pointsOuter = vtkSmartPointer<vtkPoints>::New();
	vtkSmartPointer<vtkPoints> pointsInner = vtkSmartPointer<vtkPoints>::New();

	double3 sliceMidPoint;
	double pointFromSlice[3];

	// calculate the middle of the slice
	for (int i = 0; i < slice->GetNumberOfPoints(); i++) {
		slice->GetPoint(i, pointFromSlice);

		sliceMidPoint.x += pointFromSlice[0];
		sliceMidPoint.y += pointFromSlice[1];
		sliceMidPoint.z += pointFromSlice[2];
	}
	sliceMidPoint.x /= slice->GetNumberOfPoints();
	sliceMidPoint.y /= slice->GetNumberOfPoints();
	sliceMidPoint.z /= slice->GetNumberOfPoints();

	// iterate over all cells in the slice
	for (int i = 0; i < slice->GetNumberOfCells(); i++) {
		double normal[3], angle, p0[3], p1[3], p2[3];
		
		// get normal from points of the triangle
		normal = getNormalFromCell(slice->GetCell(i));

		// calculate the angle between the normal of a cell and the heart main axis
		angle = vtkMath::DegreesFromRadians(vtkMath::AngleBetweenVectors(sliceMainAxis, normal));

		// fetch cell and cast it to a triangle
		vtkSmartPointer<vtkCell> cell = slice->GetCell(i);
		vtkSmartPointer<vtkTriangle> triangle = vtkTriangle::SafeDownCast(cell);

		// get the points of the current cell (a triangle)
		triangle->GetPoints()->GetPoint(0, p0);
		triangle->GetPoints()->GetPoint(1, p1);
		triangle->GetPoints()->GetPoint(2, p2);

		// decide to which area the cell belongs
		if ((angle >= -15.0) && (angle <= 15.0)) {
			// point is on the top of the slice

			// insert into a cellArray
			trianglesTop->InsertNextCell(cell);

			// insert points into a points-collection
			pointsTop->InsertNextPoint(p0);
			pointsTop->InsertNextPoint(p1);
			pointsTop->InsertNextPoint(p2);
		} else if ((angle >= 165.0) || (angle <= -165.0)) {
			// point is on the bottom of the slice

			// insert triangle into a cellArray
			trianglesBottom->InsertNextCell(cell);

			// insert points into a points-collection
			pointsBottom->InsertNextPoint(p0);
			pointsBottom->InsertNextPoint(p1);
			pointsBottom->InsertNextPoint(p2);
		} else if (std::isnan(angle)) {
			// angle is "nan"
			throw std::runtime_error("divideIntoAreas(): The angle is NaN.\n");
		} else {
			// point is on the sides
			
			// decide if inner or outer wall based on distance to the sliceMidPoint

			// convert one of the edge-points to double3
			double3 edgePoint;
			edgePoint.x = p0[0];
			edgePoint.y = p0[1];
			edgePoint.z = p0[2];

			if (calcDistance(edgePoint, sliceMidPoint) > treshold) {
				// cell is part of the outer wall

				// insert triangle into a cellArray
				trianglesOuter->InsertNextCell(cell);

				// insert points into a points-collection
				pointsOuter->InsertNextPoint(p0);
				pointsOuter->InsertNextPoint(p1);
				pointsOuter->InsertNextPoint(p2);
			} else {
				// cell is part of the inner wall

				// insert triangle into a cellArray
				trianglesInner->InsertNextCell(cell);

				// insert points into a points-collection
				pointsInner->InsertNextPoint(p0);
				pointsInner->InsertNextPoint(p1);
				pointsInner->InsertNextPoint(p2);
			}
		} //end if
	} //end for

	// add the geometry and topology to the areas
	areaDataTop->SetPoints(pointsTop);
	areaDataTop->SetPolys(trianglesTop);

	areaDataBottom->SetPoints(pointsBottom);
	areaDataBottom->SetPolys(trianglesBottom);
	
	areaDataOuter->SetPoints(pointsOuter);
	areaDataOuter->SetPolys(trianglesOuter);
	
	areaDataInner->SetPoints(pointsInner);
	areaDataInner->SetPolys(trianglesInner);
 }

A input slice file and the resulting areas (written to file with a vtkXMLPolyDataWriter) are found here: https://filebin.net/u4fjq7by5cqeucps

2

Can you post a complete compilable example and a CMakeLists file? Also test data.

Thanks,

Bill

3

For situations like this, I usually compute point scalars using vtkArrayCalculator and extract relevant parts of the mesh using vtkClipPolyData.

4

Unfortunately, I can’t post a complete compilable example because this code is part of greater project. But I can add some test-input data in the above post.

Is the code understandable as it stands there? If not, let me know which part isn’t clear and I try to explain it.
(double3 is a struct to store double[3] in because we use an older c++ version, which doesn’t support the passing of double[3] as a return value: struct double3 {double x, y, z;}; )

Best, Armin

5

@lassoan Do you have an example on how to do it at hand?

I have now assigned each cell a double-value as an attribute which shows me, if the cell is part of the top/bottom/inner/outer.
How do I create a clipper which uses this information? I can’t find an example which covers this exact usage of an vtkClipPolyData

vtkSmartPointer<vtkDoubleArray> decisionArray = vtkSmartPointer<vtkDoubleArray>::New();
decisionArray->SetName("decArr");

	// decide to which area the cell belongs
	if ((angle >= -30.0) && (angle <= 30.0)) {
		// point is on the top of the slice
		decisionArray->InsertNextValue(0.0);
	} else if ((angle >= 150.0) || (angle <= -150.0)) {
		// point is on the bottom of the slice
		decisionArray->InsertNextValue(1.0);
	} else if (std::isnan(angle)) {
		// angle is "nan"
		throw std::runtime_error(divideIntoAreas(): The angle is NaN.\n");
	} else {
		// point is on the sides
		// decide if inner or outer based on distance to the sliceMidPoint

		// convert one of the edge-points to double3
		double3 edgePoint;
		edgePoint.x = p0[0];
		edgePoint.y = p0[1];
		edgePoint.z = p0[2];

		if (calcDistance(edgePoint, sliceMidPoint) > threshold) {
			// cell is part of the outer wall
			decisionArray->InsertNextValue(2.0);
		} else {
			// cell is part of the inner wall
			decisionArray->InsertNextValue(3.0);
		}
	} //end if

	slice->GetCellData()->AddArray(decisionArray);

	// TODO: create and use vtkClipPolyData here

Best, Armin

6

This is a code snippet to split surface by a closed loop defined by selectionPoints.

surface = vtk.vtkPolyData()
selectionPoints = vtk.vtkPoints()
...
loop = vtk.vtkSelectPolyData()
loop.SetLoop(selectionPoints)
loop.GenerateSelectionScalarsOn()
loop.SetInputData(surface)
loop.SetSelectionModeToLargestRegion()

clip = vtk.vtkClipPolyData()
clip.InsideOutOn()
clip.SetInputConnection(loop.GetOutputPort())
clip.GenerateClippedOutputOn()
clip.Update()
clippedOutputs = [clip.GetOutput(), clip.GetClippedOutput()]
7

Thanks for the snippet. I tried with the loop, but that doesn’t really work in my case as it seems. I get an error that the algorithm can’t follow the edges.
Also, I still don’t get how to use the array I created in post before. I just need to extract those triangles from the inputPolyData that satisfy a condition. With the additional array I thought there might be a function where I can say something like: extract the triangle if triangle->GetValueOfDecisionArrayForCurrentTriangle == 0.0 (or 1.0/2.0/3.0). Is there no way to do that?

My code snippet in the first post extracts some triangles and i’m able to put them into a PolyData-object, but I’m unable to keep the geometry and topology intact, so all points are stacked on top of each other.

Best, Armin

8

I found a solution and wanted to share it because I feel that it might be useful for someone.

First I created vtkIdTypeArrays:

// arrays to store the selected ids in
vtkSmartPointer<vtkIdTypeArray> idsTop = vtkSmartPointer<vtkIdTypeArray>::New();
vtkSmartPointer<vtkIdTypeArray> idsBottom = vtkSmartPointer<vtkIdTypeArray>::New();
vtkSmartPointer<vtkIdTypeArray> idsOuter = vtkSmartPointer<vtkIdTypeArray>::New();
vtkSmartPointer<vtkIdTypeArray> idsInner = vtkSmartPointer<vtkIdTypeArray>::New();
    	
// id-arrays are all one-dimensional
idsTop->SetNumberOfComponents(1);
idsBottom->SetNumberOfComponents(1);
idsOuter->SetNumberOfComponents(1);
idsInner->SetNumberOfComponents(1);

Then I just “collected” all the IDs of those cells that satisfied my condition. This should work with any condition.

// iterate over all cells in the slice
for (int i = 0; i < slice->GetNumberOfCells(); i++) {
	double cellNormal[3], angle, p0[3], p1[3], p2[3];
	double3 normal;
	
	// get normal from points of the triangle
	normal = getNormalFromCell(slice->GetCell(i));
	
	// calculate the angle between the normal of a cell and the heart main axis
	angle = vtkMath::DegreesFromRadians(vtkMath::AngleBetweenVectors(sliceMainAxis, normal));

	// fetch cell and cast it to a triangle
	vtkSmartPointer<vtkCell> cell = slice->GetCell(i);
	vtkSmartPointer<vtkTriangle> triangle = vtkTriangle::SafeDownCast(cell);

	// get the points of the current cell (a triangle)
	triangle->GetPoints()->GetPoint(0, p0);
	triangle->GetPoints()->GetPoint(1, p1);
	triangle->GetPoints()->GetPoint(2, p2);
	
	// decide to which area the cell belongs
	if ((angle >= -15.0) && (angle <= 15.0)) {
		// point is on the top of the slice
		idsTop->InsertNextValue(i);
	} else if ((angle >= 165.0) || (angle <= -165.0)) {
		// point is on the bottom of the slice
		idsBottom->InsertNextValue(i);
	} else if (std::isnan(angle)) {
		// angle is "nan"
		throw std::runtime_error("CBCompareSlices::DivideIntoAreas(): The angle between the heart's main axis and the cell-normal is nan.\n");
	} else {
		// point is on the sides

		// convert one of the edge-points to double3
		double3 edgePoint;
		edgePoint.x = p0[0];
		edgePoint.y = p0[1];
		edgePoint.z = p0[2];

		// decide if endo or epi based on distance to the sliceMidPoint
		if (calcDistance(edgePoint, sliceMidPoint) > threshold) {
			// cell is part of the outer wall
			idsOuter->InsertNextValue(i);
		} else {
			// cell is part of the inner wall
			idsInner->InsertNextValue(i);
		}
	} //end if
} //end for

//call the extraction method
vtkSmartPointer<vtkUnstructuredGrid> resultTop = performCellExtraction(slice, idsTop);
vtkSmartPointer<vtkUnstructuredGrid> resultBottom = performCellExtraction(slice, idsBottom);
vtkSmartPointer<vtkUnstructuredGrid> resultOuter = performCellExtraction(slice, idsOuter);
vtkSmartPointer<vtkUnstructuredGrid> resultInner = performCellExtraction(slice, idsInner);

At this point I have four lists containing the IDs of those cells that satisfy the conditions.
All I have to do now, is extracting these cells based on their IDs. I therefore wrote an additional method which just does this extraction:

vtkSmartPointer<vtkUnstructuredGrid> performCellExtraction(vtkSmartPointer<vtkPolyData> slice, vtkSmartPointer<vtkIdTypeArray> idArray) {
	vtkSmartPointer<vtkSelectionNode> selectionNode = vtkSmartPointer<vtkSelectionNode>::New();
	vtkSmartPointer<vtkSelection> selection = vtkSmartPointer<vtkSelection>::New();
	vtkSmartPointer<vtkExtractSelection> extractSelection = vtkSmartPointer<vtkExtractSelection>::New();
	vtkSmartPointer<vtkUnstructuredGrid> selected = vtkSmartPointer<vtkUnstructuredGrid>::New();

    // configure selection nodes
    selectionNode->SetFieldType(vtkSelectionNode::CELL);
    selectionNode->SetContentType(vtkSelectionNode::INDICES);
    selectionNode->SetSelectionList(idArray);

    // add the selection node to the selction
    selection->AddNode(selectionNode);

    // set input data for the extractor and perform the extraction
    extractSelection->SetInputData(0, slice);
    extractSelection->SetInputData(1, selection);
    extractSelection->Update();

    // retrieve extraction result
    selected->ShallowCopy(extractSelection->GetOutput());

	return selected;
}

I then have my unstructured grid divided in the areas I wanted it to be divided to and can work further on these parts.

I hope that this helps.
Best, Armin