# Extracting Triangles from PolyData while keeping geometry and topology

**URL:** https://discourse.vtk.org/t/extracting-triangles-from-polydata-while-keeping-geometry-and-topology/1557
**Category:** Support
**Created:** [August 12, 2019, 12:59pm UTC](https://discourse.vtk.org/t/extracting-triangles-from-polydata-while-keeping-geometry-and-topology/1557 "2019-08-12T12:59:13Z")
**Posts on this page:** 8
**Page:** 1

<div class="post-metadata">

### Author: ![burgadon](https://discourse.vtk.org/letter_avatar_proxy/v4/letter/b/9e8a1a/32.png) [@burgadon](https://discourse.vtk.org/u/burgadon)
#### Post date: [August 12, 2019, 12:59pm UTC](https://discourse.vtk.org/t/extracting-triangles-from-polydata-while-keeping-geometry-and-topology/1557/1 "2019-08-12T12:59:13Z")

</div>

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](https://filebin.net/u4fjq7by5cqeucps)

---

<div class="post-metadata">

### Author: ![lorensen](https://discourse.vtk.org/user_avatar/discourse.vtk.org/lorensen/32/20_2.png) [@lorensen](https://discourse.vtk.org/u/lorensen)
#### Post date: [August 12, 2019, 6:15pm UTC](https://discourse.vtk.org/t/extracting-triangles-from-polydata-while-keeping-geometry-and-topology/1557/2 "2019-08-12T18:15:15Z")

</div>

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

Thanks,

Bill

---

<div class="post-metadata">

### Author: ![lassoan](https://discourse.vtk.org/user_avatar/discourse.vtk.org/lassoan/32/50_2.png) [@lassoan](https://discourse.vtk.org/u/lassoan)
#### Post date: [August 13, 2019, 1:54am UTC](https://discourse.vtk.org/t/extracting-triangles-from-polydata-while-keeping-geometry-and-topology/1557/3 "2019-08-13T01:54:48Z")

</div>

For situations like this, I usually compute point scalars using [vtkArrayCalculator](https://vtk.org/doc/nightly/html/classvtkArrayCalculator.html) and extract relevant parts of the mesh using [vtkClipPolyData](https://vtk.org/doc/nightly/html/classvtkClipPolyData.html).

---

<div class="post-metadata">

### Author: ![burgadon](https://discourse.vtk.org/letter_avatar_proxy/v4/letter/b/9e8a1a/32.png) [@burgadon](https://discourse.vtk.org/u/burgadon)
#### Post date: [August 13, 2019, 10:43am UTC](https://discourse.vtk.org/t/extracting-triangles-from-polydata-while-keeping-geometry-and-topology/1557/4 "2019-08-13T10:43:19Z")

</div>

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

---

<div class="post-metadata">

### Author: ![burgadon](https://discourse.vtk.org/letter_avatar_proxy/v4/letter/b/9e8a1a/32.png) [@burgadon](https://discourse.vtk.org/u/burgadon)
#### Post date: [August 13, 2019, 10:51am UTC](https://discourse.vtk.org/t/extracting-triangles-from-polydata-while-keeping-geometry-and-topology/1557/5 "2019-08-13T10:51:12Z")

</div>

@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

---

<div class="post-metadata">

### Author: ![lassoan](https://discourse.vtk.org/user_avatar/discourse.vtk.org/lassoan/32/50_2.png) [@lassoan](https://discourse.vtk.org/u/lassoan)
#### Post date: [August 14, 2019, 6:09am UTC](https://discourse.vtk.org/t/extracting-triangles-from-polydata-while-keeping-geometry-and-topology/1557/6 "2019-08-14T06:09:13Z")

</div>

> [@burgadon](#):
>
> Do you have an example on how to do it at hand?

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

```python
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()]

```

---

<div class="post-metadata">

### Author: ![burgadon](https://discourse.vtk.org/letter_avatar_proxy/v4/letter/b/9e8a1a/32.png) [@burgadon](https://discourse.vtk.org/u/burgadon)
#### Post date: [August 14, 2019, 11:36am UTC](https://discourse.vtk.org/t/extracting-triangles-from-polydata-while-keeping-geometry-and-topology/1557/7 "2019-08-14T11:36:19Z")

</div>

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

---

<div class="post-metadata">

### Author: ![burgadon](https://discourse.vtk.org/letter_avatar_proxy/v4/letter/b/9e8a1a/32.png) [@burgadon](https://discourse.vtk.org/u/burgadon)
#### Post date: [August 22, 2019, 10:29am UTC](https://discourse.vtk.org/t/extracting-triangles-from-polydata-while-keeping-geometry-and-topology/1557/8 "2019-08-22T10:29:13Z")

</div>

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