# Remove displayed points delimited by 6 planes

**URL:** https://discourse.vtk.org/t/remove-displayed-points-delimited-by-6-planes/5183
**Category:** Support
**Created:** [February 13, 2021, 6:15pm UTC](https://discourse.vtk.org/t/remove-displayed-points-delimited-by-6-planes/5183 "2021-02-13T18:15:33Z")
**Posts on this page:** 8
**Page:** 1

<div class="post-metadata">

### Author: ![Luis\_Carlos\_Carneiro](https://discourse.vtk.org/user_avatar/discourse.vtk.org/luis_carlos_carneiro/32/1360_2.png) [@Luis\_Carlos\_Carneiro](https://discourse.vtk.org/u/Luis_Carlos_Carneiro)
#### Post date: [February 13, 2021, 6:15pm UTC](https://discourse.vtk.org/t/remove-displayed-points-delimited-by-6-planes/5183/1 "2021-02-13T18:15:33Z")

</div>

Dear All,

Which the best way to remove displayed points of a volume, delimited by 6 planes?

I want to give 6x3 points (x,y,z each point) defining the planes and remove the points inside. Each plane defined by 3 points.

Thanks,

Luís Gonçalves

---

<div class="post-metadata">

### Author: ![Paulo\_Carvalho](https://discourse.vtk.org/user_avatar/discourse.vtk.org/paulo_carvalho/32/370_2.png) [@Paulo\_Carvalho](https://discourse.vtk.org/u/Paulo_Carvalho)
#### Post date: [February 13, 2021, 10:24pm UTC](https://discourse.vtk.org/t/remove-displayed-points-delimited-by-6-planes/5183/2 "2021-02-13T22:24:41Z")

</div>

Hello, Luís,

Here is the code I use to remove points from an arbitrary hexahedron. Actually this code can be used for any polyhedron provided all faces are counter-clockwise when viewed from inside:

```cpp
-------------------------vector3d.h-----------------------------
#ifndef VECTOR3D_H
#define VECTOR3D_H

class Vector3D {
public:
	double x, y, z;

	Vector3D operator-(Vector3D p) const;

	Vector3D operator+(Vector3D p) const;

	Vector3D cross(Vector3D p) const;

	double dot(Vector3D p) const;

	double norm() const;
};

Vector3D operator *( double scalar, const Vector3D& vector );

//make a point type with the same structure of the vector for clarity
typedef Vector3D Vertex3D;

#endif // VECTOR3D_H

```

```cpp
---------------------------vector3d.cpp---------------------------------------------

#include "vector3d.h"
#include <cmath>

Vector3D Vector3D::operator-(Vector3D p) const {
	return Vector3D{x - p.x, y - p.y, z - p.z};
}

Vector3D Vector3D::operator+(Vector3D p) const
{
	return Vector3D{x + p.x, y + p.y, z + p.z};
}

Vector3D Vector3D::cross(Vector3D p) const {
	return Vector3D{
		y * p.z - p.y * z,
				z * p.x - p.z * x,
				x * p.y - p.x * y
	};
}

double Vector3D::dot(Vector3D p) const {
	return x * p.x + y * p.y + z * p.z;
}

double Vector3D::norm() const {
	return std::sqrt(x*x + y*y + z*z);
}

Vector3D operator *(double scalar, const Vector3D & vector)
{
	Vector3D result;
	result.x = scalar * vector.x;
	result.y = scalar * vector.y;
	result.z = scalar * vector.z;
	return result;
}

```

```cpp
-------------------------face3d.h----------------------------

#ifndef FACE3D_H
#define FACE3D_H

#include "geometry/vector3d.h"
#include <vector>

class Face3D {
public:
	Face3D(); //making the face with 4 vertexes by default

	std::vector<Vertex3D> v;

	Vector3D normal() const;

	/** Computes the distance from point to the plane that contains this face. */
	double distance( const Vertex3D& point ) const;
};

#endif // FACE3D_H

```

```cpp
-----------------------face3d.cpp--------------------------------------

#include "face3d.h"
#include <cassert>

Face3D::Face3D() : v(4){}

// TODO: treat degenerate cases (e.g. a face may collapse into a triangle, edge or point).
// current GeoGrid constructors do not generate degenerate cells, but in the future, they might.

Vector3D Face3D::normal() const {
	assert(v.size() > 2);
	Vector3D dir1 = v[1] - v[0];
	Vector3D dir2 = v[2] - v[0];
	// dir1 cross dir2 or dir2 cross dir1 depends on whether the 3D coordinate system is right-handed or left-handed
	Vector3D n = dir2.cross(dir1);
	double d = n.norm();
	return Vector3D{n.x / d, n.y / d, n.z / d};
}

double Face3D::distance(const Vertex3D & point) const
{
	//Get the vector normal (n) to the plane of this face.
	Vector3D n = normal();
	//Plane of triangle:
	// n . u = n . v[0], where u is any (x,y,z) contained in the plane (. == dot product).
	// v[0] is a vertex of the face.
	//Line that contains the point with direction given by the normal (n):
	// u = point + tn (t is a scalar)
	//Intersection point satisfies:
	// n . ( point + tn ) = n . v[0]
	// n . point + t = n . v[0]
	double t = n.dot( v[0] ) - n.dot( point );
	//Hence, the point that intersects the line and the plane, p0, is:
	Vertex3D p0 = point + t * n;
	//and the distance is point-to-p0 distance
	return ( point - p0 ).norm();
}

```

```cpp

--------------------------------THE METHOD----------------------------------------

/**
 * This method assumes all faces are COUNTER-CLOCKWISE when viewed from INSIDE.
 *
 * Vertex id elements in the vId[8] array forming the geometry of a cell:
 *
 * J
 * |
 * |
 *
 * 3-----------2 face 0 = 0 1 2 3
 * /| /| face 1 = 4 7 6 5
 * / | / | face 2 = 0 3 7 4
 * / | / | face 3 = 1 5 6 2
 * 7--------- 6 | face 4 = 0 4 5 1
 * | | | | face 5 = 3 2 6 7
 * | 0----|------1 --- I
 * | / | /
 * | / | /
 * |/ | /
 * 4--------- 5
 *
 * /
 * /
 * K
 *
 *
 *
 */

bool Util::isInside(const Vertex3D & p, const std::vector<Face3D> & fs)
{
	double bound = -1e-15; // use -1e-15 to exclude boundaries
	for (const Face3D& f : fs) {
		Vector3D p2f = f.v[0] - p; // any point of f cloud be used
		double d = p2f.dot( f.normal() );
		d /= p2f.norm(); // for numeric stability
		if ( d < bound )
			return false;
	}
    return true;
}

```

I hope this helps.

abraço,

Paulo

---

<div class="post-metadata">

### Author: ![Paulo\_Carvalho](https://discourse.vtk.org/user_avatar/discourse.vtk.org/paulo_carvalho/32/370_2.png) [@Paulo\_Carvalho](https://discourse.vtk.org/u/Paulo_Carvalho)
#### Post date: [February 13, 2021, 10:33pm UTC](https://discourse.vtk.org/t/remove-displayed-points-delimited-by-6-planes/5183/3 "2021-02-13T22:33:43Z")

</div>

Oops… Vertex3D = Vector3D: `typedef Vector3D Vertex3D;`

---

<div class="post-metadata">

### Author: ![Luis\_Carlos\_Carneiro](https://discourse.vtk.org/user_avatar/discourse.vtk.org/luis_carlos_carneiro/32/1360_2.png) [@Luis\_Carlos\_Carneiro](https://discourse.vtk.org/u/Luis_Carlos_Carneiro)
#### Post date: [February 13, 2021, 10:35pm UTC](https://discourse.vtk.org/t/remove-displayed-points-delimited-by-6-planes/5183/4 "2021-02-13T22:35:09Z")

</div>

Thanks, for the code.

But I have similar code in C+CUDA for doing for unstructured images (Meshes) but for 6 random planes/faces (each with the information of the side of the plane to remove). For structured images (voxels) it is similar and easier.

What I want is if there is some filter in VTK to do something similar to a displayed image but in Python (VTK).

---

<div class="post-metadata">

### Author: ![will.schroeder](https://discourse.vtk.org/user_avatar/discourse.vtk.org/will.schroeder/32/233_2.png) [@will.schroeder](https://discourse.vtk.org/u/will.schroeder)
#### Post date: [February 14, 2021, 12:05pm UTC](https://discourse.vtk.org/t/remove-displayed-points-delimited-by-6-planes/5183/5 "2021-02-14T12:05:19Z")

</div>

Off the top of my head, vtkExtractGeometry (which I believe takes an unstructured grid as input) should work. You have to provide it with a vtkImplicitFunction, in this case vtkPlanes. vtkPlanes must define a convex space. Also, each plane is defined by a point and a normal, so you’ll have to convert the three points into this form.

---

<div class="post-metadata">

### Author: ![Paulo\_Carvalho](https://discourse.vtk.org/user_avatar/discourse.vtk.org/paulo_carvalho/32/370_2.png) [@Paulo\_Carvalho](https://discourse.vtk.org/u/Paulo_Carvalho)
#### Post date: [February 14, 2021, 12:30pm UTC](https://discourse.vtk.org/t/remove-displayed-points-delimited-by-6-planes/5183/6 "2021-02-14T12:30:39Z")

</div>

Complementing Will’s answer, you can get the normal vector from the three points defining the plane by using the cross product between two non-parallel vectors. Chose one point as the common origin of the two vectors. Each vector is then defined by this common origin and another point.

In the code I provided above, the method `Vector3D::cross()` can be used to compute the normal (translating it to Python should be easy). Again, ordering is important here too:

![image](https://discourse.vtk.org/uploads/default/original/2X/6/692cac4255edd1d9276942ea202942ef7c4d514a.png)

---

<div class="post-metadata">

### Author: ![Luis\_Carlos\_Carneiro](https://discourse.vtk.org/user_avatar/discourse.vtk.org/luis_carlos_carneiro/32/1360_2.png) [@Luis\_Carlos\_Carneiro](https://discourse.vtk.org/u/Luis_Carlos_Carneiro)
#### Post date: [February 14, 2021, 6:23pm UTC](https://discourse.vtk.org/t/remove-displayed-points-delimited-by-6-planes/5183/7 "2021-02-14T18:23:41Z")

</div>

Tanks, Will.

---

<div class="post-metadata">

### Author: ![Luis\_Carlos\_Carneiro](https://discourse.vtk.org/user_avatar/discourse.vtk.org/luis_carlos_carneiro/32/1360_2.png) [@Luis\_Carlos\_Carneiro](https://discourse.vtk.org/u/Luis_Carlos_Carneiro)
#### Post date: [February 23, 2021, 1:23pm UTC](https://discourse.vtk.org/t/remove-displayed-points-delimited-by-6-planes/5183/8 "2021-02-23T13:23:11Z")

</div>

Dear All,

Can somebody advise me?  
I have 6 planes defined by 3 points each. And a 4th point defining the side of the plane I want to select.  
I do the cross product to find the normal to the plane. I do the dot product between the (4th point minus onepointofplane) and the normal of the plane to find the right direction of the normal.  
I want to remove the polygon volume defined by those 6 planes from a preexistent volume (self.volume). For that I delete self.volume and create a new one supposed without the polygon volume. I used self.reader.GetOutputPort() that contains the original image. I have stability errors e.g system instability. Perhaps problems with OpenGL.  
Need I to delete the original volume? There is a better way to delete the polygon volume in the original image directly? What the errors below?

```
            plane1 = vtk.vtkPlane()
            plane1.SetOrigin(self.x1[0],self.y1[0],self.z1[0])
            p3p1=np.array([self.x1[2]-self.x1[0],self.y1[2]-self.y1[0],self.z1[2]-self.z1[0]])
            p2p1=np.array([self.x1[1]-self.x1[0],self.y1[1]-self.y1[0],self.z1[1]-self.z1[0]])
            cross1=np.cross(p3p1,p2p1)
            cross1=cross1/np.linalg.norm(cross1)
            dot1=cross1[0]*(self.sidex[0]-self.x1[0])+cross1[1]*(self.sidey[0]-self.y1[0])+cross1[2]*(self.sidez[0]-self.z1[0])
            if (dot1 < 0.0):
                cross1=-cross1
            plane1.SetNormal(cross1[0],cross1[1],cross1[2])

```

# code for plane2 to plane6

```
            union = vtk.vtkImplicitBoolean()
            union.AddFunction(plane1)
            union.AddFunction(plane2)
            union.AddFunction(plane3)
            union.AddFunction(plane4)
            union.AddFunction(plane5)
            union.AddFunction(plane6)
            union.SetOperationType(0) #union
            extract = vtk.vtkExtractGeometry()
            extract.SetInputConnection(self.reader.GetOutputPort())
            extract.SetImplicitFunction(union)              
            volumeMapper = vtk.vtkOpenGLGPUVolumeRayCastMapper()

            volumeMapper.SetInputConnection(extract.GetOutputPort())

            volumeColor = vtk.vtkColorTransferFunction()
            volumeColor.AddRGBPoint(0,0.0,0.0,0.0)
            volumeColor.AddRGBPoint(255,0.0,0.0,1.0)

            volumeScalarOpacity = vtk.vtkPiecewiseFunction()
            volumeScalarOpacity.AddPoint(0,0.0)
            volumeScalarOpacity.AddPoint(255,1.0)

            volumeGradientOpacity = vtk.vtkPiecewiseFunction()
            volumeGradientOpacity.AddPoint(0,1.0)
            volumeGradientOpacity.AddPoint(255,1.0)
            volumeProperty = vtk.vtkVolumeProperty()
            volumeProperty.SetColor(volumeColor)
            volumeProperty.SetScalarOpacity(volumeScalarOpacity)
            volumeProperty.SetGradientOpacity(volumeGradientOpacity)
            volumeProperty.SetInterpolationTypeToLinear()
            volumeProperty.ShadeOff()
            volumeProperty.SetAmbient(0.6)
            volumeProperty.SetDiffuse(0.6)
            volumeProperty.SetSpecular(0.1)
            renderer.RemoveVolume(self.volume)
            renderWindow.Render()
            self.volume = vtk.vtkVolume()
            self.volume.SetMapper(volumeMapper)
            self.volume.SetProperty(volumeProperty)
            renderer.AddVolume(self.volume)
            renderWindow.Render()

```

Thanks,

Luís Gonçalves