From 4584c83f2289b12cd69aacdc1b7c433cd9e86bb9 Mon Sep 17 00:00:00 2001
From: joannacirillo <cirillojoanna@gmail.com>
Date: Mon, 18 Jul 2022 16:10:33 +0200
Subject: [PATCH 1/2] feat(polygon, testPolygon): add functions to polygon

Translation from cpp implementation and add test to those functions
Create functions as static functions and in publicAPI
Fix pointInPolygon and triangulate to return the triangles
Updated ts definition
---
 Sources/Common/Core/PriorityQueue/index.js    |   14 +-
 Sources/Common/DataModel/Cell/index.js        |    8 +-
 Sources/Common/DataModel/Polygon/Constants.js |   13 +-
 Sources/Common/DataModel/Polygon/index.d.ts   |  257 +++-
 Sources/Common/DataModel/Polygon/index.js     | 1071 ++++++++++++++---
 .../DataModel/Polygon/test/testPolygon.js     |  545 +++++++++
 .../General/ContourTriangulator/helper.js     |    9 +-
 .../General/ContourTriangulator/index.js      |    4 +-
 Sources/Filters/General/PaintFilter/index.js  |   10 +-
 .../Filters/General/TriangleFilter/index.js   |    9 +-
 10 files changed, 1723 insertions(+), 217 deletions(-)
 create mode 100644 Sources/Common/DataModel/Polygon/test/testPolygon.js

diff --git a/Sources/Common/Core/PriorityQueue/index.js b/Sources/Common/Core/PriorityQueue/index.js
index ed22e222e00..acc41524981 100644
--- a/Sources/Common/Core/PriorityQueue/index.js
+++ b/Sources/Common/Core/PriorityQueue/index.js
@@ -15,18 +15,22 @@ function vtkPriorityQueue(publicAPI, model) {
     model.elements.splice(i, 0, { priority, element });
   };
 
+  publicAPI.deleteById = (id) => {
+    model.elements = model.elements.filter(({ element }) => element !== id);
+  };
+
   publicAPI.pop = () => {
     if (model.elements.length > 0) {
-      return model.elements.shift().element;
+      const element = model.elements.reduce((prev, current) =>
+        prev.priority > current.priority ? prev : current
+      );
+      publicAPI.deleteById(element.element);
+      return element.element;
     }
 
     return null;
   };
 
-  publicAPI.deleteById = (id) => {
-    model.elements = model.elements.filter(({ element }) => element.id !== id);
-  };
-
   publicAPI.length = () => model.elements.length;
 }
 
diff --git a/Sources/Common/DataModel/Cell/index.js b/Sources/Common/DataModel/Cell/index.js
index c1f537ff6e9..885cfd0521c 100644
--- a/Sources/Common/DataModel/Cell/index.js
+++ b/Sources/Common/DataModel/Cell/index.js
@@ -104,8 +104,12 @@ function vtkCell(publicAPI, model) {
     cell.initialize(model.points, model.pointsIds);
   };
 
-  publicAPI.getCellDimension = () => {}; // virtual
-  publicAPI.intersectWithLine = (p1, p2, tol, t, x, pcoords, subId) => {}; // virtual
+  publicAPI.getCellDimension = () => {
+    macro.vtkErrorMacro('vtkCell.getCellDimension is not implemented.');
+  }; // virtual
+  publicAPI.intersectWithLine = (p1, p2, tol, t, x, pcoords, subId) => {
+    macro.vtkErrorMacro('vtkCell.intersectWithLine is not implemented.');
+  }; // virtual
   publicAPI.evaluatePosition = (
     x,
     closestPoint,
diff --git a/Sources/Common/DataModel/Polygon/Constants.js b/Sources/Common/DataModel/Polygon/Constants.js
index 71cc68d9f11..b14c87f72d5 100644
--- a/Sources/Common/DataModel/Polygon/Constants.js
+++ b/Sources/Common/DataModel/Polygon/Constants.js
@@ -6,6 +6,15 @@ export const PolygonWithPointIntersectionState = {
   FAILURE: -1,
   OUTSIDE: 0,
   INSIDE: 1,
-  INTERSECTION: 2,
-  ON_LINE: 3,
 };
+export const VTK_DBL_EPSILON = 2.2204460492503131e-16;
+
+export const PolygonWithCellIntersectionState = {
+  NO_INTERSECTION: 0,
+  POINT_INTERSECTION: 1,
+  LINE_INTERSECTION: 2,
+  OVERLAP: 3,
+  INCLUDED: 4,
+};
+
+export default { PolygonWithPointIntersectionState };
diff --git a/Sources/Common/DataModel/Polygon/index.d.ts b/Sources/Common/DataModel/Polygon/index.d.ts
index 9842609ed8e..7a0f1322227 100644
--- a/Sources/Common/DataModel/Polygon/index.d.ts
+++ b/Sources/Common/DataModel/Polygon/index.d.ts
@@ -1,66 +1,262 @@
 import { vtkObject } from "../../../interfaces";
 import { Bounds, TypedArray, Vector3 } from "../../../types";
+import vtkPoints from "../../Core/Points";
+import vtkCell from "../Cell";
 
 export interface IPolygonInitialValues {
-	firstPoint?: Vector3,
-	pointCount?: number,
-	tris?: Vector3[],
+	pointCount?: number;
+	tris?: Vector3[];
 }
 
 /**
  * Different states which pointInPolygon could return.
  */
-export enum PolygonIntersectionState {
+export enum PolygonWithPointIntersectionState {
 	FAILURE,
 	OUTSIDE,
 	INSIDE,
-	INTERSECTION,
-	ON_LINE,
+}
+
+/**
+ * Different states that intersectWith2DConvexCell could return.
+ */
+export enum PolygonWithCellIntersectionState {
+	NO_INTERSECTION,
+	LINE_INTERSECTION,
+	POINT_INTERSECTION,
+	OVERLAP,
+	INCLUDED
+}
+
+interface IIntersectWithLine {
+	intersection: boolean;
+	betweenPoints: boolean;
+	t: number;
+	x: Vector3;
+}
+
+interface IDistanceToPolygon {
+	t: number,
+	distance: number
 }
 
 export interface vtkPolygon extends vtkObject {
+	/**
+	 * Set the polygon's points
+	 * Points must be ordered in counterclockwise order
+	 * @param {Vector3[]|Array<number>} points The polygon's points.
+	 * @param {Array<number>} pointIds pointIds 
+	 */
+	setPoints(points: Vector3[]|Array<number>, pointIds?: Array<number>): void;
+
+	/**
+ 	 * Get the bounds for this polygon as [xmin, xmax, ymin, ymax, zmin, zmax].
+	 * 	@return {Bounds} bounds
+	 */
+	getBounds(): Bounds
+
+	/**
+	 * Computes the polygon normal
+	 * @return {number} norm of normal (before normalization)
+	 */
+	computeNormal(): number;
 
 	/**
-	 * Get the array of triangles that triangulate the polygon.
+	 * Determine whether a point is inside a polygon. The function uses a winding
+	 * number calculation generalized to the 3D plane one which the polygon
+	 * resides. Returns OUTSIDE if point is not in the polygon; INSIDE if it is inside. Can
+	 * also return FAILURE to indicate a degenerate polygon (points non coplanar or on a line).
+	 * This implementation is inspired by Dan Sunday's algorithm found in the book Practical
+	 * Geometry Algorithms.
+	 * @param {Vector3} point Point to check
+	 * @return {PolygonWithPointIntersectionState} type of intersection
 	 */
-	getPointArray(): Vector3[];
+	pointInPolygon(point: Vector3): PolygonWithPointIntersectionState;
 
 	/**
-	 * Set the polygon's points.
-	 * @param {Vector3[]} points The polygon's points.
+	 * Compute ear triangulation of current polygon
+	 * The polygon must be convex and have at least 3 points
+	 * @return {boolean} whether triangulation failed or not
 	 */
-	setPoints(points: Vector3[]): void;
+	triangulate(): boolean;
 
 	/**
-	 * Triangulate this polygon. 
-	 * The output data must be accessed through `getPointArray`.
-	 * The output data contains points by group of three: each three-group
-	 * defines one triangle.
+	 * Returns the centroid of this polygon
+	 * @return {Vector3} centroid
 	 */
-	triangulate(): void;
+	computeCentroid(): Vector3;
 
+	/**
+	 * Returns the area of the polygon
+	 * @return {number} area
+	 */
+	computeArea(): number;
+
+	/**
+	 * determine the distance from a point to a polygon.
+ 	* @param {Vector3} x
+ 	* @param {Vector3} closest filled with the closest point in the polygon
+ 	* @return {IDistanceToPolygon} object containing the distance (distance) and the tolerance with wich the distance is given (t)
+ 	*/
+	distanceToPolygon(x: Vector3, closest: Vector3): IDistanceToPolygon;
+
+	/**
+	 * Returns whether the polygon is convex or not
+	 * Returns false for degenerate polygon
+	 * @return {boolean} is convex or not
+	 */
+	isConvex(): boolean;
+
+	/**
+	 * Interpolates functions with polygon points
+	 * @param {Vector3} point point to compute the interpolation on
+	 * @param {boolean} useMVCInterpolation
+	 * @return weights corresponding to each point of polygon parametrizing the given point
+	 */
+	interpolateFunctions(
+		point: Vector3,
+		useMVCInterpolation: boolean
+	): number[];
+
+	/**
+	 * Computes intersection of polygon with a line defined by two points
+	 * @param {Vector3} x1 first point of line
+	 * @param {Vector3} x2 second point of line
+	 * @return intersection point coordinates
+	 */
+	intersectWithLine(x1: Vector3, x2: Vector3): IIntersectWithLine;
+
+	/**
+	 * Computes intersection of polygon with another cell.
+	 * It can be a line, a point, no intersection or coincident
+	 * Note: Expects both polygons/cell to be convex
+	 * @param {vtkCell} cell polygon or any object extending from vtkCell with which to compute intersection
+	 * Note : the function intersectWithLine need to be implemented on the class of the cell given
+	 * @return {PolygonWithCellIntersectionState} type of intersection
+	 */
+	intersectConvex2DCells(
+		cell: vtkCell
+	): PolygonWithCellIntersectionState;
 }
 
+// ---------------------------------------------------
+/**
+ * Compute the normal of a polygon and return its squared norm.
+ * @param {vtkPoints} points
+ * @param {Vector3} normal
+ * @return {number}
+ */
+export function getNormal(
+	points: vtkPoints,
+	normal: Vector3
+): number;
+
+/**
+ * Get the bounds for these points as [xmin, xmax, ymin, ymax,zmin, zmax].
+ * @param {vtkPoints} points
+ * @return {Bounds}
+ */
+export function getBounds(points: vtkPoints): Bounds;
+
+/**
+ * Determines whether a polygon is convex
+ * @param {vtkPoints} points vtkPoints defining the polygon
+ * @return {boolean} whether the polygon is convex or not
+ */
+export function isConvex(points: vtkPoints): boolean;
+
+/**
+ * Given a set of points, computes the centroid of the corresponding polygon
+ * @param {vtkPoints} points vtkPoints defining the polygon
+ * @param {Vector3} normal normal to the polygon of which the centroid is computed
+ * @return {Vector3} centroid. Returns null for degenerate polygon
+ */
+export function computeCentroid(points: vtkPoints, normal: Vector3): Vector3;
+
+/**
+ * Given a set of points, computes the area of the corresponding polygon
+ * @param {vtkPoints} points vtkPoints defining the polygon
+ * @param {Vector3} normal normal to the polygon of which the centroid is computed
+ * @return {number} area of polygon
+ */
+export function computeArea(points: vtkPoints, normal: Vector3): number;
+
+/**
+ * Given a set of points, determine the distance from a point to a polygon.
+ * @param {Vector3} x
+ * @param {vtkPoints} points vtkPoints defining the polygon
+ * @param {Vector3} closest filled with the closest point in the polygon
+ * @return {IDistanceToPolygon} object containing the distance (distance) and the tolerance with wich the distance is given (t)
+ */
+export function distanceToPolygon(x: Vector3, points: vtkPoints, closest: Vector3): IDistanceToPolygon;
+
 /**
  * Determine whether a point is inside a polygon. The function uses a winding
  * number calculation generalized to the 3D plane one which the polygon
- * resides. Returns 0 if point is not in the polygon; 1 if it is inside. Can
- * also return -1 to indicate a degenerate polygon. This implementation is
+ * resides. Returns OUTSIDE if point is not in the polygon; INSIDE if it is inside. Can
+ * also return FAILURE to indicate a degenerate polygon. This implementation is
  * inspired by Dan Sunday's algorithm found in the book Practical Geometry
  * Algorithms.
  *
  * @param {Vector3} point Point to check
- * @param {Array<Number>|TypedArray} vertices Vertices of the polygon
+ * @param {Array<number>|TypedArray} vertices Vertices of the polygon
  * @param {Bounds} bounds Bounds of the vertices
  * @param {Vector3} normal Normal vector of the polygon
- * @returns {PolygonIntersectionState} Integer indicating the type of intersection
+ * @return {PolygonWithPointIntersectionState} Integer indicating the type of intersection
  */
 export function pointInPolygon(
-  point: Vector3,
-  vertices: Array<number>|TypedArray,
-  bounds: Bounds,
-  normal: Vector3
-): PolygonIntersectionState;
+	point: Vector3,
+	vertices: Array<number> | TypedArray,
+	bounds: Bounds,
+	normal: Vector3
+): PolygonWithPointIntersectionState;
+
+/**
+ * Given a set of points that define a polygon, determines whether a line defined
+ * by two points intersect with the polygon. There can be no intersection, a point
+ * intersection or a line intersection.
+ * @param {Vector3} p1 first point of the line
+ * @param {Vector3} p2 second point of the line
+ * @param {vtkPoints} points points defining the polygon
+ * @param {Vector3} normal normal to the polygon
+ * @return {IIntersectWithLine} type of intersection
+ */
+export function intersectWithLine(
+	p1: Vector3,
+	p2: Vector3,
+	points: vtkPoints,
+	normal: Vector3
+): IIntersectWithLine;
+
+/**
+ * Given a set of points that define a polygon and another polygon, computes their
+ * intersection. It can be a line, a point, no intersection or coincident
+ * Note: Expects both polygons need to be convex
+ * @param {vtkCell} cell polygon or any object extending from vtkCell with which to compute intersection
+ * Note : the function intersectWithLine need to be implemented on the class of the cell given
+ * @param {vtkPoints} points points defining the polygon
+ * @param {Vector3} normal normal to the polygon
+ * @return {PolygonWithCellIntersectionState} type of intersection
+ */
+export function intersectConvex2DCells(
+	cell: vtkCell,
+	points: vtkPoints,
+	normal: Vector3
+): PolygonWithCellIntersectionState;
+
+/**
+ * Given a set of points, computes the weights corresponding to the interpolation of the
+ * given point with regard to the points of the polygon. The returned array corresponds to
+ * the weights and therefore its size is the number of points in the polygon
+ * @param {Vector3} point point we want the interpolation of
+ * @param {vtkPoints} points points defining the polygon
+ * @param {boolean} useMVCInterpolation whether to use MVC interpolation
+ */
+export function interpolateFunctions(
+	point: Vector3,
+	points: vtkPoints,
+	useMVCInterpolation: boolean
+): Array<number>;
 
 /**
  * Method used to decorate a given object (publicAPI+model) with vtkPolygon characteristics.
@@ -69,7 +265,11 @@ export function pointInPolygon(
  * @param model object on which data structure will be bounds (protected)
  * @param {IPolygonInitialValues} [initialValues] (default: {})
  */
-export function extend(publicAPI: object, model: object, initialValues?: IPolygonInitialValues): void;
+export function extend(
+	publicAPI: object,
+	model: object,
+	initialValues?: IPolygonInitialValues
+): void;
 
 /**
  * Method used to create a new instance of vtkPolygon.
@@ -79,15 +279,14 @@ export function newInstance(initialValues?: IPolygonInitialValues): vtkPolygon;
 
 /**
  * vtkPolygon represents a 2D n-sided polygon.
- * 
+ *
  * The polygons cannot have any internal holes, and cannot self-intersect.
  * Define the polygon with n-points ordered in the counter-clockwise direction.
  * Do not repeat the last point.
  */
 export declare const vtkPolygon: {
-	newInstance: typeof newInstance,
+	newInstance: typeof newInstance;
 	extend: typeof extend;
 	// static
-
 };
 export default vtkPolygon;
diff --git a/Sources/Common/DataModel/Polygon/index.js b/Sources/Common/DataModel/Polygon/index.js
index 7bc937e507c..f01de7aaf83 100644
--- a/Sources/Common/DataModel/Polygon/index.js
+++ b/Sources/Common/DataModel/Polygon/index.js
@@ -3,14 +3,23 @@ import * as vtkMath from 'vtk.js/Sources/Common/Core/Math';
 import vtkLine from 'vtk.js/Sources/Common/DataModel/Line';
 import vtkPlane from 'vtk.js/Sources/Common/DataModel/Plane';
 import vtkPriorityQueue from 'vtk.js/Sources/Common/Core/PriorityQueue';
-import vtkBoundingBox from 'vtk.js/Sources/Common/DataModel/BoundingBox';
 import { IntersectionState } from 'vtk.js/Sources/Common/DataModel/Line/Constants';
+import vtkBoundingBox from 'vtk.js/Sources/Common/DataModel/BoundingBox';
+import vtkCell from 'vtk.js/Sources/Common/DataModel/Cell';
+import vtkPoints from 'vtk.js/Sources/Common/Core/Points';
+
 import {
-  PolygonWithPointIntersectionState,
   EPSILON,
   FLOAT_EPSILON,
   TOLERANCE,
-} from './Constants';
+  PolygonWithPointIntersectionState,
+  VTK_DBL_EPSILON,
+  PolygonWithCellIntersectionState,
+} from 'vtk.js/Sources/Common/DataModel/Polygon/Constants';
+
+// ----------------------------------------------------------------------------
+// vtkPolygon methods
+// ----------------------------------------------------------------------------
 
 // ----------------------------------------------------------------------------
 // Global methods
@@ -29,7 +38,6 @@ function pointLocation(axis0, axis1, p0, p1, point) {
 }
 
 //------------------------------------------------------------------------------
-
 function pointInPolygon(point, vertices, bounds, normal) {
   // Do a quick bounds check to throw out trivial cases.
   // winding plane.
@@ -64,13 +72,24 @@ function pointInPolygon(point, vertices, bounds, normal) {
 
   for (let i = 0; i < vertices.length; ) {
     // Check coincidence to polygon vertices
+    // double* p0 = pts + 3 * i;
     p0[0] = vertices[i++];
     p0[1] = vertices[i++];
     p0[2] = vertices[i++];
+
     if (vtkMath.distance2BetweenPoints(point, p0) <= tol2) {
       return PolygonWithPointIntersectionState.INSIDE;
     }
 
+    if (i < vertices.length) {
+      p1[0] = vertices[i];
+      p1[1] = vertices[i + 1];
+      p1[2] = vertices[i + 2];
+    } else {
+      p1[0] = vertices[0];
+      p1[1] = vertices[1];
+      p1[2] = vertices[2];
+    }
     // Check coincidence to polygon edges
     const { distance, t } = vtkLine.distanceToLine(point, p0, p1);
     if (distance <= tol2 && t > 0.0 && t < 1.0) {
@@ -81,7 +100,7 @@ function pointInPolygon(point, vertices, bounds, normal) {
   // If here, begin computation of the winding number. This method works for
   // points/polygons arbitrarily oriented in 3D space.  Hence a projection
   // onto one of the x-y-z coordinate planes using the maximum normal
-  // component. The computation will be performed in the (axis0, axis1) plane.
+  // component. The computation will be performed in the (axis0,axis1) plane.
   let axis0;
   let axis1;
   if (Math.abs(normal[0]) > Math.abs(normal[1])) {
@@ -142,77 +161,735 @@ function pointInPolygon(point, vertices, bounds, normal) {
     : PolygonWithPointIntersectionState.INSIDE;
 }
 
-// ---------------------------------------------------
-/**
- * Simple utility method for computing polygon bounds.
- * Returns the sum of the squares of the dimensions.
- * Requires a poly with at least one point.
- *
- * @param {Array<Number>|TypedArray<Number>} poly
- * @param {vtkPoints} points
- * @param {Bound} bounds
- */
-export function getBounds(poly, points, bounds) {
-  const n = poly.length;
-  const p = [];
-
-  points.getPoint(poly[0], p);
-  bounds[0] = p[0];
-  bounds[1] = p[0];
-  bounds[2] = p[1];
-  bounds[3] = p[1];
-  bounds[4] = p[2];
-  bounds[5] = p[2];
-
-  for (let j = 1; j < n; j++) {
-    points.getPoint(poly[j], p);
-    vtkBoundingBox.addPoint(bounds, p);
-  }
-  const length = vtkBoundingBox.getLength(bounds);
-  return vtkMath.dot(length, length);
+//------------------------------------------------------------------------------
+function getBounds(points) {
+  return points.getBounds();
 }
 
-// ---------------------------------------------------
-/**
- * Compute the normal of a polygon and return its norm.
- *
- * TBD: This does the same thing as vtkPolygon.computeNormal,
- * but in a more generic way. Maybe we can keep the public
- * static method somehow and have the private method use it instead.
- *
- * @param {Array<Number>|TypedArray<Number>} poly
- * @param {vtkPoints} points
- * @param {Vector3} normal
- * @returns {Number}
- */
-export function getNormal(poly, points, normal) {
-  normal.length = 3;
-  normal[0] = 0.0;
-  normal[1] = 0.0;
-  normal[2] = 0.0;
-
+//------------------------------------------------------------------------------
+function getNormal(points, normal) {
+  const n = points.getNumberOfPoints();
+  const nn = [0.0, 0.0, 0.0];
   const p0 = [];
   let p1 = [];
   let p2 = [];
-  const v1 = [];
-  const v2 = [];
 
-  points.getPoint(poly[0], p0);
-  points.getPoint(poly[1], p1);
+  points.getPoint(0, p0);
+  points.getPoint(1, p1);
 
-  for (let j = 2; j < poly.length; j++) {
-    points.getPoint(poly[j], p2);
+  const v1 = [];
+  const v2 = [];
+  for (let j = 2; j < n; j++) {
+    points.getPoint(j, p2);
     vtkMath.subtract(p2, p1, v1);
     vtkMath.subtract(p0, p1, v2);
 
-    const n = [0, 0, 0];
-    vtkMath.cross(v1, v2, n);
-    vtkMath.add(normal, n, normal);
+    nn[0] += v1[1] * v2[2] - v1[2] * v2[1];
+    nn[1] += v1[2] * v2[0] - v1[0] * v2[2];
+    nn[2] += v1[0] * v2[1] - v1[1] * v2[0];
 
     [p1, p2] = [p2, p1];
   }
 
-  return vtkMath.normalize(normal);
+  const norm = Math.sqrt(vtkMath.dot(nn, nn));
+  normal[0] = nn[0] / norm;
+  normal[1] = nn[1] / norm;
+  normal[2] = nn[2] / norm;
+  return norm;
+}
+
+//------------------------------------------------------------------------------
+function isConvex(points) {
+  let i;
+  const v = [
+    [0, 0, 0],
+    [0, 0, 0],
+    [0, 0, 0],
+  ];
+  let tmp;
+  const a = [0, 0, 0];
+  let aMag;
+  const b = [0, 0, 0];
+  let bMag;
+  const n = [0, 0, 0];
+  const ni = [0, 0, 0];
+  let nComputed = false;
+  const numPts = points.getNumberOfPoints();
+
+  if (numPts < 3) {
+    return false;
+  }
+
+  if (numPts === 3) {
+    return true;
+  }
+
+  points.getPoint(0, v[1]);
+  points.getPoint(1, v[2]);
+
+  for (i = 0; i <= numPts; i++) {
+    tmp = v[0];
+    v[0] = v[1];
+    v[1] = v[2];
+    v[2] = tmp;
+
+    points.getPoint((i + 2) % numPts, v[2]);
+
+    // order is important!!! to maintain consistency with polygon vertex order
+    a[0] = v[2][0] - v[1][0];
+    a[1] = v[2][1] - v[1][1];
+    a[2] = v[2][2] - v[1][2];
+    b[0] = v[0][0] - v[1][0];
+    b[1] = v[0][1] - v[1][1];
+    b[2] = v[0][2] - v[1][2];
+
+    if (!nComputed) {
+      aMag = vtkMath.norm(a);
+      bMag = vtkMath.norm(b);
+      if (aMag > VTK_DBL_EPSILON && bMag > VTK_DBL_EPSILON) {
+        vtkMath.cross(a, b, n);
+        nComputed =
+          vtkMath.norm(n) > VTK_DBL_EPSILON * (aMag < bMag ? bMag : aMag);
+      }
+    } else {
+      vtkMath.cross(a, b, ni);
+      if (vtkMath.norm(ni) > VTK_DBL_EPSILON && vtkMath.dot(n, ni) < 0) {
+        return false;
+      }
+    }
+  }
+
+  return true;
+}
+
+//------------------------------------------------------------------------------
+function computeCentroid(points, normal) {
+  const numPts = points.getNumberOfPoints();
+  // Strategy:
+  // - Compute centroid of projected polygon on (x,y) if polygon is projectible, (x,z) otherwise
+  // - Accumulate signed projected area as well, which is needed in the centroid's formula
+  // - Infer 3rd dimension using polygon's normal.
+  let i = 0;
+  if (numPts < 2) {
+    macro.vtkWarningMacro('Not enough points to compute centroid.');
+  }
+  // Handle to the doubled area of the projected polygon on (x,y) or (x,z) if the polygon
+  // projected on (x,y) is degenerate.
+  let a = 0.0;
+  const p0 = [...points.getPoint(0), 0, 0, 0];
+
+  let xOffset = 0;
+  let yOffset = 0;
+  let c = [];
+  // Checking if the polygon is colinear with z axis.
+  // If it is, the centroid computation axis shall be shifted
+  // because the projected polygon on (x,y) is degenerate.
+  const z = [0, 0, 1];
+  vtkMath.cross(normal, z, c);
+  // If the normal is orthogonal with z axis, the projected polygon is then a line...
+  if (Math.abs(c[0] * c[0] + c[1] * c[1] + c[2] * c[2] - 1.0) <= EPSILON) {
+    yOffset = 1;
+    const y = [0, 1, 0];
+    vtkMath.cross(normal, y, c);
+    // If the normal is orthogonal with y axis, the projected polygon is then a line...
+    if (Math.abs(c[0] * c[0] + c[1] * c[1] + c[2] * c[2] - 1.0) <= EPSILON) {
+      xOffset = 1;
+    }
+  }
+
+  c = [0, 0, 0];
+
+  let maxabsdet = 0;
+  for (i = 0; i < numPts; i++) {
+    const point = points.getPoint(i);
+    p0[3 * !(i % 2)] = point[0];
+    p0[3 * !(i % 2) + 1] = point[1];
+    p0[3 * !(i % 2) + 2] = point[2];
+    const det =
+      p0[3 * (i % 2) + xOffset] * p0[3 * !(i % 2) + 1 + yOffset] -
+      p0[3 * !(i % 2) + xOffset] * p0[3 * (i % 2) + 1 + yOffset];
+    c[xOffset] += (p0[xOffset] + p0[3 + xOffset]) * det;
+    c[1 + yOffset] += (p0[1 + yOffset] + p0[4 + yOffset]) * det;
+    a += det;
+    maxabsdet = Math.abs(det) > maxabsdet ? Math.abs(det) : maxabsdet;
+  }
+  if (Math.abs(a) < EPSILON * maxabsdet) {
+    // Polygon is degenerate
+    macro.vtkWarningMacro('Warning: polygon is degenerate, no centroid.');
+    return null;
+  }
+  c[xOffset] /= 3.0 * a;
+  c[1 + yOffset] /= 3.0 * a;
+  c[2 - xOffset - yOffset] =
+    (1.0 / normal[2 - xOffset - yOffset]) *
+    (-normal[xOffset] * c[xOffset] -
+      normal[1 + yOffset] * c[1 + yOffset] +
+      vtkMath.dot(normal, p0));
+  return c;
+}
+
+//------------------------------------------------------------------------------
+function computeArea(points, normal) {
+  const numPts = points.getNumberOfPoints();
+  if (numPts < 3) {
+    return 0.0;
+  }
+  if (normal.length === 0) {
+    return NaN;
+  }
+
+  let area = 0.0;
+
+  // Select the projection direction
+  const nx = normal[0] > 0.0 ? normal[0] : -normal[0]; // abs x-coord
+  const ny = normal[1] > 0.0 ? normal[1] : -normal[1]; // abs y-coord
+  const nz = normal[2] > 0.0 ? normal[2] : -normal[2]; // abs z-coord
+
+  // coord is the index of the axis with biggest normal coordinate
+  const a = nx > nz ? 0 : 2;
+  const b = ny > nz ? 1 : 2;
+  const coord = nx > ny ? a : b;
+
+  // compute area of the 2D projection
+  const x0 = [];
+  const x1 = [];
+  const x2 = [];
+
+  for (let i = 0; i < numPts; i++) {
+    points.getPoint(i, x0);
+    points.getPoint((i + 1) % numPts, x1);
+    points.getPoint((i + 2) % numPts, x2);
+    switch (coord) {
+      case 0:
+        area += x1[1] * (x2[2] - x0[2]);
+        break;
+      case 1:
+        area += x1[0] * (x2[2] - x0[2]);
+        break;
+      case 2:
+        area += x1[0] * (x2[1] - x0[1]);
+        break;
+      default:
+        area += 0;
+        break;
+    }
+  }
+
+  // scale to get area before projection
+  switch (coord) {
+    case 0:
+      area /= 2.0 * nx;
+      break;
+    case 1:
+      area /= 2.0 * ny;
+      break;
+    case 2:
+      area /= 2.0 * nz;
+      break;
+    default:
+      area /= 1;
+      break;
+  }
+  return Math.abs(area);
+}
+
+//------------------------------------------------------------------------------
+function distanceToPolygon(x, points, closest) {
+  const outObj = { t: Number.MIN_VALUE, distance: 0 };
+  const bounds = points.getBounds();
+  // First check to see if the point is inside the polygon
+  // do a quick bounds check
+  if (
+    x[0] >= bounds[0] &&
+    x[0] <= bounds[1] &&
+    x[1] >= bounds[2] &&
+    x[1] <= bounds[3] &&
+    x[2] >= bounds[4] &&
+    x[2] <= bounds[5]
+  ) {
+    const n = [0, 0, 0];
+    getNormal(points, n);
+    if (pointInPolygon(x, points.getData(), bounds, n)) {
+      closest[0] = x[0];
+      closest[1] = x[1];
+      closest[2] = x[2];
+      return outObj;
+    }
+  }
+
+  // Not inside, compute the distance of the point to the edges.
+  let minDist2 = Number.MAX_VALUE;
+  const p0 = [0, 0, 0];
+  const p1 = [0, 0, 0];
+  const c = [0, 0, 0];
+  let t = 0;
+  let dist2;
+  const numPts = points.getNumberOfPoints();
+  for (let i = 0; i < numPts; i++) {
+    points.getPoint(i, p0);
+    points.getPoint((i + 1) % numPts, p1);
+    const distToLine = vtkLine.distanceToLine(x, p0, p1, c);
+    t = distToLine.t;
+    dist2 = distToLine.distance;
+    if (dist2 < minDist2) {
+      minDist2 = dist2;
+      closest[0] = c[0];
+      closest[1] = c[1];
+      closest[2] = c[2];
+    }
+  }
+
+  outObj.t = t;
+  outObj.distance = Math.sqrt(minDist2);
+  return outObj;
+}
+
+//------------------------------------------------------------------------------
+function arePolygonsCoincident(points1, points2) {
+  // Check if the planes defined by the polygons points are coincident
+  const origin1 = [0, 0, 0];
+  points1.getPoint(0, origin1);
+  const origin2 = [0, 0, 0];
+  points2.getPoint(0, origin2);
+  const lineBetweenPlanes = vtkMath.subtract(origin1, origin2, []);
+  const normal1 = [0, 0, 0];
+  const normal2 = [0, 0, 0];
+  getNormal(points1, normal1);
+  getNormal(points2, normal2);
+  const dot1 = vtkMath.dot(lineBetweenPlanes, normal1);
+  const dot2 = vtkMath.dot(lineBetweenPlanes, normal2);
+
+  if (!(dot1 === 0 && dot2 === 0)) return null;
+
+  const bounds1 = points1.getBounds();
+  const bounds2 = points2.getBounds();
+  const corner11 = [0, 0, 0];
+  const corner12 = [0, 0, 0];
+  const corner21 = [0, 0, 0];
+  const corner22 = [0, 0, 0];
+  vtkBoundingBox.computeCornerPoints(bounds1, corner11, corner12);
+  vtkBoundingBox.computeCornerPoints(bounds2, corner21, corner22);
+
+  // Intersect the two bounding boxes. If the result is one of the bounding boxes,
+  // one is included in the other
+  const bounds1Copy = [...bounds1];
+  const intersection = vtkBoundingBox.intersect(bounds1Copy, bounds2);
+
+  // If the bboxes intersect and the polygons are coplanar, the polygons are
+  // either overlapping or included in one another
+  if (intersection) {
+    if (
+      vtkMath.areEquals(bounds1Copy, bounds1) ||
+      vtkMath.areEquals(bounds1Copy, bounds2)
+    )
+      return PolygonWithCellIntersectionState.INCLUDED;
+    return PolygonWithCellIntersectionState.OVERLAP;
+  }
+  return null;
+}
+
+//------------------------------------------------------------------------------
+// Create a local s-t coordinate system for a polygon. The point p0 is
+// the origin of the local system, p10 is s-axis vector, and p20 is the
+// t-axis vector. (These are expressed in the modelling coordinate system and
+// are vectors of dimension [3].) The values l10 and l20 are the lengths of
+// the vectors p10 and p20
+function parameterizePolygon(points, normal) {
+  let i;
+  let j;
+  let s;
+  let t;
+  const p = [0, 0, 0];
+  const p1 = [0, 0, 0];
+  const p2 = [0, 0, 0];
+  const sbounds = [0, 0];
+  const tbounds = [0, 0];
+  const x1 = [0, 0, 0];
+  const x2 = [0, 0, 0];
+  const numPts = points.getNumberOfPoints();
+
+  const p0 = [];
+  const p10 = [];
+  let l10 = 0;
+  const p20 = [];
+  let l20 = 0;
+
+  const outObj = {
+    p0,
+    p10,
+    l10,
+    p20,
+    l20,
+  };
+
+  if (numPts < 3) {
+    return macro.vtkWarningMacro(
+      'Cannot parameterize polygon with less than 2 points'
+    );
+  }
+
+  //  This is a two pass process: first create a p' coordinate system
+  //  that is then adjusted to ensure that the polygon points are all in
+  //  the range 0<=s,t<=1.  The p' system is defined by the polygon normal,
+  //  first vertex and the first edge.
+  //
+  points.getPoint(0, x1);
+  points.getPoint(1, x2);
+  for (i = 0; i < 3; i++) {
+    p0[i] = x1[i];
+    p10[i] = x2[i] - x1[i];
+  }
+  vtkMath.cross(normal, p10, p20);
+
+  // Determine lengths of edges
+  //
+  l10 = vtkMath.dot(p10, p10);
+  l20 = vtkMath.dot(p20, p20);
+  if (l10 === 0.0 || l20 === 0.0) {
+    return outObj;
+  }
+
+  //  Now evaluate all polygon points to determine min/max parametric
+  //  coordinate values.
+  //
+  // first vertex has (s,t) = (0,0)
+  sbounds[0] = 0.0;
+  sbounds[1] = 0.0;
+  tbounds[0] = 0.0;
+  tbounds[1] = 0.0;
+
+  for (i = 1; i < numPts; i++) {
+    points.getPoint(i, x1);
+    for (j = 0; j < 3; j++) {
+      p[j] = x1[j] - p0[j];
+    }
+    s = (p[0] * p10[0] + p[1] * p10[1] + p[2] * p10[2]) / l10;
+    t = (p[0] * p20[0] + p[1] * p20[1] + p[2] * p20[2]) / l20;
+    sbounds[0] = s < sbounds[0] ? s : sbounds[0];
+    sbounds[1] = s > sbounds[1] ? s : sbounds[1];
+    tbounds[0] = t < tbounds[0] ? t : tbounds[0];
+    tbounds[1] = t > tbounds[1] ? t : tbounds[1];
+  }
+
+  //  Re-evaluate coordinate system
+  //
+  for (i = 0; i < 3; i++) {
+    p1[i] = p0[i] + sbounds[1] * p10[i] + tbounds[0] * p20[i];
+    p2[i] = p0[i] + sbounds[0] * p10[i] + tbounds[1] * p20[i];
+    p0[i] = p0[i] + sbounds[0] * p10[i] + tbounds[0] * p20[i];
+    p10[i] = p1[i] - p0[i];
+    p20[i] = p2[i] - p0[i];
+  }
+  l10 = vtkMath.norm(p10);
+  l20 = vtkMath.norm(p20);
+  outObj.l10 = l10;
+  outObj.l20 = l20;
+
+  return outObj;
+}
+
+//------------------------------------------------------------------------------
+// Compute interpolation weights using mean value coordinate.
+function interpolateFunctionsUsingMVC(point, points) {
+  const numPts = points.getNumberOfPoints();
+  const weights = vtkMath.createArray(numPts);
+
+  // create local array for storing point-to-vertex vectors and distances
+  const uVec = vtkMath.createArray(numPts * 3);
+  const dist = vtkMath.createArray(numPts);
+  for (let i = 0; i < numPts; i++) {
+    const pt = points.getPoint(i);
+    // point-to-vertex vector
+    uVec[3 * i] = pt[0] - point[0];
+    uVec[3 * i + 1] = pt[1] - point[1];
+    uVec[3 * i + 2] = pt[2] - point[2];
+    // distance
+    dist[i] = vtkMath.norm([uVec[3 * i], uVec[3 * i + 1], uVec[3 * i + 2]]);
+
+    // handle special case when the point is really close to a vertex
+    if (dist[i] <= EPSILON) {
+      weights[i] = 1.0;
+      return weights;
+    }
+
+    uVec[3 * i] /= dist[i];
+    uVec[3 * i + 1] /= dist[i];
+    uVec[3 * i + 2] /= dist[i];
+  }
+
+  // Now loop over all vertices to compute weight
+  // w_i = ( tan(theta_i/2) + tan(theta_(i+1)/2) ) / dist_i
+  // To do consider the simplification of
+  // tan(alpha/2) = (1-cos(alpha))/sin(alpha)
+  //              = (d0*d1 - cross(u0, u1))/(2*dot(u0,u1))
+  const tanHalfTheta = vtkMath.createArray(numPts);
+  for (let i = 0; i < numPts; i++) {
+    const i1 = (i + 1) % numPts;
+
+    const l = Math.sqrt(
+      vtkMath.distance2BetweenPoints(
+        [uVec[3 * i], uVec[3 * i + 1], uVec[3 * i + 2]],
+        [uVec[3 * i1], uVec[3 * i1 + 1], uVec[3 * i1 + 2]]
+      )
+    );
+    const theta = 2.0 * Math.asin(l / 2.0);
+
+    // special case where x lies on an edge
+    if (Math.PI - theta < 0.001) {
+      weights[i] = dist[i1] / (dist[i] + dist[i1]);
+      weights[i1] = 1 - weights[i];
+      return weights;
+    }
+
+    tanHalfTheta[i] = Math.tan(theta / 2.0);
+  }
+
+  // Normal case
+  for (let i = 0; i < numPts; i++) {
+    let i1 = i - 1;
+    if (i1 === -1) {
+      i1 = numPts - 1;
+    }
+
+    weights[i] = (tanHalfTheta[i] + tanHalfTheta[i1]) / dist[i];
+  }
+
+  // normalize weight
+  let sum = 0.0;
+  for (let i = 0; i < numPts; i++) {
+    sum += weights[i];
+  }
+
+  if (Math.abs(sum) <= EPSILON) {
+    return weights;
+  }
+
+  for (let i = 0; i < numPts; i++) {
+    weights[i] /= sum;
+  }
+
+  return weights;
+}
+
+//------------------------------------------------------------------------------
+function interpolateFunctions(point, points, useMVCInterpolation) {
+  // Compute interpolation weights using mean value coordinate.
+  if (useMVCInterpolation) {
+    return interpolateFunctionsUsingMVC(point, points);
+  }
+
+  const weights = [];
+
+  // Compute interpolation weights using 1/r**2 normalized sum.
+  let i;
+  const numPts = points.getNumberOfPoints();
+  let sum;
+  const pt = [0, 0, 0];
+
+  for (sum = 0.0, i = 0; i < numPts; i++) {
+    points.getPoint(i, pt);
+    weights[i] = vtkMath.distance2BetweenPoints(point, pt);
+    if (weights[i] === 0.0) {
+      // exact hit
+      for (let j = 0; j < numPts; j++) {
+        weights[j] = 0.0;
+      }
+      weights[i] = 1.0;
+      return weights;
+    }
+    weights[i] = 1.0 / weights[i];
+    sum += weights[i];
+  }
+
+  for (i = 0; i < numPts; i++) {
+    weights[i] /= sum;
+  }
+  return weights;
+}
+
+// ------------------------------------------------------------------------------
+// Given a 3D point "point" return inside, outside cell, or failure
+// Evaluate parametric coordinates,
+// distance squared of point "point" to cell (in particular, the sub-cell indicated),
+// closest point on cell to "point" (unless closestPoint is null, in which case, the
+// closest point and dist2 are not found), and interpolation weights in cell.
+function evaluatePosition(point, points, normal) {
+  const cp = [0, 0, 0];
+  const ray = [0, 0, 0];
+
+  const res = {
+    intersection: false,
+    dist: Infinity,
+    pcoords: [],
+    weights: [],
+    closestPoint: [],
+  };
+
+  const { p0, p10, l10, p20, l20 } = parameterizePolygon(points, normal);
+  res.weights = interpolateFunctions(point, points, false);
+
+  vtkPlane.projectPoint(point, p0, normal, cp);
+
+  for (let i = 0; i < 3; i++) {
+    ray[i] = cp[i] - p0[i];
+  }
+  res.pcoords[0] = vtkMath.dot(ray, p10) / (l10 * l10);
+  res.pcoords[1] = vtkMath.dot(ray, p20) / (l20 * l20);
+  res.pcoords[2] = 0.0;
+
+  if (
+    res.pcoords[0] >= 0.0 &&
+    res.pcoords[0] <= 1.0 &&
+    res.pcoords[1] >= 0.0 &&
+    res.pcoords[1] <= 1.0 &&
+    pointInPolygon(cp, points.getData(), points.getBounds(), normal) ===
+      PolygonWithPointIntersectionState.INSIDE
+  ) {
+    res.closestPoint = [...cp];
+    res.dist = vtkMath.distance2BetweenPoints(point, res.closestPoint);
+    res.intersection = true;
+    return res;
+  }
+
+  // If here, point is outside of polygon, so need to find distance to boundary
+  //
+
+  let t;
+  let dist2;
+  const closest = [0, 0, 0];
+  const pt1 = [0, 0, 0];
+  const pt2 = [0, 0, 0];
+
+  const numPts = points.getNumberOfPoints();
+  for (let minDist2 = Number.MAX_VALUE, i = 0; i < numPts; i++) {
+    points.getPoint(i, pt1);
+    points.getPoint((i + 1) % numPts, pt2);
+    dist2 = vtkLine.distanceToLine(point, pt1, pt2, t, closest);
+    if (dist2 < minDist2) {
+      res.closestPoint = [...closest];
+      minDist2 = dist2;
+    }
+  }
+  res.dist = dist2;
+  return res;
+}
+
+//------------------------------------------------------------------------------
+function intersectWithLine(p1, p2, points, normal) {
+  let outObj = {
+    intersection: false,
+    betweenPoints: false,
+    t: Number.MAX_VALUE,
+    x: [],
+  };
+  const tol2 = TOLERANCE * TOLERANCE;
+
+  // Intersect plane of the polygon with line
+  //
+  const data = points.getData();
+  outObj = vtkPlane.intersectWithLine(
+    p1,
+    p2,
+    [data[0], data[1], data[2]],
+    normal
+  );
+  if (!outObj.intersection) {
+    return outObj;
+  }
+
+  // Evaluate position
+  const position = evaluatePosition(outObj.x, points, normal);
+  if (!(position.intersection || position.dist <= tol2)) {
+    outObj.intersection = false;
+  }
+
+  return outObj;
+}
+
+//------------------------------------------------------------------------------
+function intersectConvex2DCells(cell, points, normal) {
+  // Intersect the six total edges of the two triangles against each other. Two points are
+  // all that are required.
+  const numPts = points.getNumberOfPoints();
+  const x0 = [];
+  const x1 = [];
+  let lineIntersection;
+  const outObj = {
+    intersection: PolygonWithCellIntersectionState.NO_INTERSECTION,
+    x0: [],
+    x1: [],
+  };
+  let idx = 0;
+  const t2 = TOLERANCE * TOLERANCE;
+
+  const numPts2 = cell.getNumberOfPoints();
+
+  const poly = new Array(numPts);
+  for (let i = 0; i < poly.length; i++) poly[i] = i;
+  const coincidence = arePolygonsCoincident(cell.getPoints(), points);
+  if (coincidence !== null) {
+    outObj.intersection = coincidence;
+    return outObj;
+  }
+
+  // Loop over edges of second polygon and intersect against first polygon
+  let i;
+  for (i = 0; i < numPts2; i++) {
+    cell.getPoints().getPoint(i, x0);
+    cell.getPoints().getPoint((i + 1) % numPts2, x1);
+
+    lineIntersection = intersectWithLine(x0, x1, points, normal);
+    if (lineIntersection.intersection) {
+      if (idx === 0) {
+        outObj.x0 = lineIntersection.x;
+        idx++;
+      } else if (
+        (x1[0] - x0[0]) * (x1[0] - x0[0]) +
+          (x1[1] - x0[1]) * (x1[1] - x0[1]) +
+          (x1[2] - x0[2]) * (x1[2] - x0[2]) >
+          t2 &&
+        !vtkMath.areEquals(lineIntersection.x, outObj.x0)
+      ) {
+        outObj.intersection =
+          PolygonWithCellIntersectionState.LINE_INTERSECTION;
+        outObj.x1 = lineIntersection.x;
+        return outObj;
+      }
+    } // if edge intersection
+  } // over all edges
+
+  for (i = 0; i < numPts; i++) {
+    points.getPoint(i, x0);
+    points.getPoint((i + 1) % numPts, x1);
+
+    lineIntersection = cell.intersectWithLine(x0, x1);
+
+    if (lineIntersection.intersection) {
+      if (idx === 0) {
+        outObj.x0 = lineIntersection.x;
+        idx++;
+      } else if (
+        (x1[0] - x0[0]) * (x1[0] - x0[0]) +
+          (x1[1] - x0[1]) * (x1[1] - x0[1]) +
+          (x1[2] - x0[2]) * (x1[2] - x0[2]) >
+          t2 &&
+        !vtkMath.areEquals(lineIntersection.x, outObj.x0)
+      ) {
+        outObj.intersection =
+          PolygonWithCellIntersectionState.LINE_INTERSECTION;
+        outObj.x1 = lineIntersection.x;
+        return outObj;
+      }
+    } // if edge intersection
+  } // over all edges
+
+  // Evaluate what we got
+  if (idx === 1) {
+    outObj.intersection = PolygonWithCellIntersectionState.POINT_INTERSECTION; // everything intersecting at single point
+    return outObj;
+  }
+  outObj.intersection = PolygonWithCellIntersectionState.NO_INTERSECTION;
+  return outObj;
 }
 
 // ----------------------------------------------------------------------------
@@ -221,9 +898,18 @@ export function getNormal(poly, points, normal) {
 
 const STATIC = {
   PolygonWithPointIntersectionState,
+  PolygonWithCellIntersectionState,
   pointInPolygon,
   getBounds,
   getNormal,
+  isConvex,
+  computeCentroid,
+  computeArea,
+  evaluatePosition,
+  intersectWithLine,
+  intersectConvex2DCells,
+  interpolateFunctions,
+  parameterizePolygon,
 };
 
 // ----------------------------------------------------------------------------
@@ -234,36 +920,33 @@ function vtkPolygon(publicAPI, model) {
   // Set our classname
   model.classHierarchy.push('vtkPolygon');
 
-  function computeNormal() {
-    const v1 = [0, 0, 0];
-    const v2 = [0, 0, 0];
-    model.normal = [0, 0, 0];
-    const anchor = [...model.firstPoint.point];
-
-    let point = model.firstPoint;
-    for (let i = 0; i < model.pointCount; i++) {
-      vtkMath.subtract(point.point, anchor, v1);
-      vtkMath.subtract(point.next.point, anchor, v2);
-
-      const n = [0, 0, 0];
-      vtkMath.cross(v1, v2, n);
-      vtkMath.add(model.normal, n, model.normal);
-
-      point = point.next;
-    }
-
-    return vtkMath.normalize(model.normal);
+  function toCorrectId(i, nbPoints) {
+    return i < 0 ? nbPoints - 1 : i % nbPoints;
   }
 
-  function computeMeasure(point) {
+  function computeMeasure(pointId, idsToHandle) {
     const v1 = [0, 0, 0];
     const v2 = [0, 0, 0];
     const v3 = [0, 0, 0];
     const v4 = [0, 0, 0];
-
-    vtkMath.subtract(point.point, point.previous.point, v1);
-    vtkMath.subtract(point.next.point, point.point, v2);
-    vtkMath.subtract(point.previous.point, point.next.point, v3);
+    const i = idsToHandle.indexOf(pointId);
+    const numPts = model.points.getNumberOfPoints();
+
+    const current = [0, 0, 0];
+    model.points.getPoint(idsToHandle[i], current);
+    const next = [0, 0, 0];
+    model.points.getPoint(
+      toCorrectId(idsToHandle[toCorrectId(i + 1, idsToHandle.length)], numPts),
+      next
+    );
+    const previous = [0, 0, 0];
+    model.points.getPoint(
+      toCorrectId(idsToHandle[toCorrectId(i - 1, idsToHandle.length)], numPts),
+      previous
+    );
+    vtkMath.subtract(current, previous, v1);
+    vtkMath.subtract(next, current, v2);
+    vtkMath.subtract(previous, next, v3);
     vtkMath.cross(v1, v2, v4);
 
     const area = vtkMath.dot(v4, model.normal);
@@ -277,36 +960,54 @@ function vtkPolygon(publicAPI, model) {
     return (perimeter * perimeter) / area;
   }
 
-  function canRemoveVertex(point) {
-    if (model.pointCount <= 3) {
+  function canRemoveVertex(pointId, idsToRemove) {
+    if (model.points.getNumberOfPoints() <= 3) {
       return true;
     }
 
-    const previous = point.previous;
-    const next = point.next;
+    const i = idsToRemove.indexOf(pointId);
 
+    const previous = [0, 0, 0];
+    const next = [0, 0, 0];
     const v = [0, 0, 0];
-    vtkMath.subtract(next.point, previous.point, v);
-
     const sN = [0, 0, 0];
+    const nextNext = [0, 0, 0];
+
+    const prevId = idsToRemove[toCorrectId(i - 1, idsToRemove.length)];
+    const nextId = idsToRemove[toCorrectId(i + 1, idsToRemove.length)];
+
+    model.points.getPoint(prevId, previous);
+    model.points.getPoint(nextId, next);
+
+    vtkMath.subtract(next, previous, v);
+
     vtkMath.cross(v, model.normal, sN);
     vtkMath.normalize(sN);
     if (vtkMath.norm(sN) === 0) {
       return false;
     }
+    model.points.getPoint(
+      idsToRemove[toCorrectId(i + 2, idsToRemove.length)],
+      nextNext
+    );
 
-    let val = vtkPlane.evaluate(sN, previous.point, next.next.point);
+    let val = vtkPlane.evaluate(sN, previous, nextNext);
     // eslint-disable-next-line no-nested-ternary
     let currentSign = val > EPSILON ? 1 : val < -EPSILON ? -1 : 0;
     let oneNegative = currentSign < 0 ? 1 : 0;
 
     for (
-      let vertex = next.next.next;
-      vertex.id !== previous.id;
-      vertex = vertex.next
+      let vertexId = idsToRemove[toCorrectId(i + 3, idsToRemove.length)];
+      vertexId !== prevId;
+      vertexId = idsToRemove[toCorrectId(vertexId + 1, idsToRemove.length)]
     ) {
-      const previousVertex = vertex.previous;
-      val = vtkPlane.evaluate(sN, previous.point, vertex.point);
+      const prevId2 =
+        idsToRemove[toCorrectId(vertexId - 1, idsToRemove.length)];
+      const previousVertex = [0, 0, 0];
+      model.points.getPoint(prevId2, previousVertex);
+      const vertex = [0, 0, 0];
+      model.points.getPoint(vertexId, vertex);
+      val = vtkPlane.evaluate(sN, previousVertex, vertex);
       // eslint-disable-next-line no-nested-ternary
       const sign = val > EPSILON ? 1 : val < -EPSILON ? -1 : 0;
 
@@ -317,10 +1018,10 @@ function vtkPolygon(publicAPI, model) {
 
         if (
           vtkLine.intersection(
-            previous.point,
-            next.point,
-            vertex.point,
-            previousVertex.point,
+            previous,
+            next,
+            vertex,
+            previousVertex,
             [0],
             [0]
           ) === IntersectionState.YES_INTERSECTION
@@ -334,102 +1035,138 @@ function vtkPolygon(publicAPI, model) {
     return oneNegative === 1;
   }
 
-  function removePoint(point, queue) {
-    model.pointCount -= 1;
+  function removePoint(pointId, queue, idsToHandle) {
+    const i = idsToHandle.indexOf(pointId);
 
-    const previous = point.previous;
-    const next = point.next;
+    const prevId = idsToHandle[toCorrectId(i - 1, idsToHandle.length)];
+    const nextId = idsToHandle[toCorrectId(i + 1, idsToHandle.length)];
+    const currentId = idsToHandle[i];
 
-    model.tris = model.tris.concat(point.point);
-    model.tris = model.tris.concat(next.point);
-    model.tris = model.tris.concat(previous.point);
+    const previous = [0, 0, 0];
+    model.points.getPoint(prevId, previous);
+    const next = [0, 0, 0];
+    model.points.getPoint(nextId, next);
+    const toRemove = [0, 0, 0];
+    model.points.getPoint(currentId, toRemove);
 
-    previous.next = next;
-    next.previous = previous;
+    queue.deleteById(prevId);
+    queue.deleteById(nextId);
 
-    queue.deleteById(previous.id);
-    queue.deleteById(next.id);
+    idsToHandle.splice(i, 1);
 
-    const previousMeasure = computeMeasure(previous);
+    const previousMeasure = computeMeasure(prevId, idsToHandle);
     if (previousMeasure > 0) {
-      queue.push(previousMeasure, previous);
+      queue.push(previousMeasure, prevId);
+    } else {
+      idsToHandle.splice(idsToHandle.indexOf(prevId), 1);
     }
 
-    const nextMeasure = computeMeasure(next);
+    const nextMeasure = computeMeasure(nextId, idsToHandle);
     if (nextMeasure > 0) {
-      queue.push(nextMeasure, next);
-    }
-
-    if (point.id === model.firstPoint.id) {
-      model.firstPoint = next;
+      queue.push(nextMeasure, nextId);
+    } else {
+      idsToHandle.splice(idsToHandle.indexOf(nextId), 1);
     }
+    return [...toRemove, ...next, ...previous];
   }
 
   function earCutTriangulation() {
-    computeNormal();
-
+    if (model.points.getNumberOfPoints() < 3) return null;
     const vertexQueue = vtkPriorityQueue.newInstance();
-    let point = model.firstPoint;
-    for (let i = 0; i < model.pointCount; i++) {
-      const measure = computeMeasure(point);
+    // The algorithm needs to remove points from the polygon, so it is handled here by storing
+    // the ids of the points that remain to handle
+    const idsToHandle = [...model.pointsIds];
+    for (let i = 0; i < idsToHandle.length; i++) {
+      const measure = computeMeasure(i, model.pointsIds);
       if (measure > 0) {
-        vertexQueue.push(measure, point);
+        vertexQueue.push(measure, i);
       }
-
-      point = point.next;
     }
-
-    while (model.pointCount > 2 && vertexQueue.length() > 0) {
-      if (model.pointCount === vertexQueue.length()) {
+    const triangles = [];
+    while (idsToHandle.length > 2 && vertexQueue.length() > 0) {
+      if (idsToHandle.length === vertexQueue.length()) {
         // convex
-        const pointToRemove = vertexQueue.pop();
-        removePoint(pointToRemove, vertexQueue);
+        const idPointToRemove = vertexQueue.pop();
+        triangles.push(
+          ...removePoint(idPointToRemove, vertexQueue, idsToHandle)
+        );
       } else {
         // concave
-        const pointToRemove = vertexQueue.pop();
-        if (canRemoveVertex(pointToRemove)) {
-          removePoint(pointToRemove, vertexQueue);
+        const idPointToRemove = vertexQueue.pop();
+        if (canRemoveVertex(idPointToRemove, idsToHandle)) {
+          triangles.push(
+            ...removePoint(idPointToRemove, vertexQueue, idsToHandle)
+          );
         }
       }
     }
 
-    return model.pointCount <= 2;
+    return idsToHandle.length <= 2 ? triangles : null;
   }
 
+  publicAPI.computeNormal = () => getNormal(model.points, model.normal);
+
+  publicAPI.getBounds = () => getBounds(model.points);
+
   publicAPI.triangulate = () => {
-    if (!model.firstPoint) {
+    if (model.points.getNumberOfPoints() === 0) {
       return null;
     }
 
     return earCutTriangulation();
   };
 
-  publicAPI.setPoints = (points) => {
-    model.pointCount = points.length;
-
-    model.firstPoint = {
-      id: 0,
-      point: points[0],
-      next: null,
-      previous: null,
-    };
-
-    let currentPoint = model.firstPoint;
-    for (let i = 1; i < model.pointCount; i++) {
-      currentPoint.next = {
-        id: i,
-        point: points[i],
-        next: null,
-        previous: currentPoint,
-      };
-      currentPoint = currentPoint.next;
+  publicAPI.setPoints = (points, pointsIds) => {
+    let pointsData = points;
+    if (Array.isArray(pointsData[0])) {
+      pointsData = points.flat(2);
     }
-
-    model.firstPoint.previous = currentPoint;
-    currentPoint.next = model.firstPoint;
+    model.points = vtkPoints.newInstance({ values: pointsData });
+    if (!pointsIds) {
+      model.pointsIds = new Array(pointsData.length / 3);
+      for (let i = 0; i < points.length; i++) model.pointsIds[i] = i;
+    } else {
+      model.pointsIds = pointsIds;
+    }
+    publicAPI.computeNormal();
   };
 
-  publicAPI.getPointArray = () => model.tris;
+  publicAPI.computeCentroid = () => computeCentroid(model.points, model.normal);
+
+  //------------------------------------------------------------------------------
+  // Compute the area of the polygon (oriented in 3D space). It uses an
+  // efficient approach where the area is computed in 2D and then projected into
+  // 3D space.
+  publicAPI.computeArea = () => computeArea(model.points, model.normal);
+
+  publicAPI.distanceToPolygon = (x, closest) =>
+    distanceToPolygon(x, model.points, closest);
+
+  publicAPI.isConvex = () => isConvex(model.points);
+
+  publicAPI.pointInPolygon = (point) =>
+    pointInPolygon(
+      point,
+      model.points.getData(),
+      publicAPI.getBounds(),
+      model.normal
+    );
+
+  //------------------------------------------------------------------------------
+  // Compute interpolation weights using 1/r**2 normalized sum or mean value
+  // coordinate.
+  publicAPI.interpolateFunctions = (point, useMVCInterpolation) =>
+    interpolateFunctions(point, model.points, useMVCInterpolation);
+
+  //------------------------------------------------------------------------------
+  //
+  // Intersect this polygon with finite line defined by p1 & p2
+  //
+  publicAPI.intersectWithLine = (x1, x2) =>
+    intersectWithLine(x1, x2, model.points, model.normal);
+
+  publicAPI.intersectConvex2DCells = (cell) =>
+    intersectConvex2DCells(cell, model.points, model.normal);
 }
 
 // ----------------------------------------------------------------------------
@@ -437,9 +1174,7 @@ function vtkPolygon(publicAPI, model) {
 // ----------------------------------------------------------------------------
 
 const DEFAULT_VALUES = {
-  firstPoint: null,
-  pointCount: 0,
-  tris: [],
+  points: [],
 };
 
 // ----------------------------------------------------------------------------
@@ -447,9 +1182,15 @@ const DEFAULT_VALUES = {
 export function extend(publicAPI, model, initialValues = {}) {
   Object.assign(model, DEFAULT_VALUES, initialValues);
 
+  vtkCell.extend(publicAPI, model, initialValues);
+
   // Build VTK API
   macro.obj(publicAPI, model);
   vtkPolygon(publicAPI, model);
+  model.normal = [];
+  if (initialValues.points) {
+    publicAPI.computeNormal();
+  }
 }
 
 // ----------------------------------------------------------------------------
diff --git a/Sources/Common/DataModel/Polygon/test/testPolygon.js b/Sources/Common/DataModel/Polygon/test/testPolygon.js
new file mode 100644
index 00000000000..de347a86408
--- /dev/null
+++ b/Sources/Common/DataModel/Polygon/test/testPolygon.js
@@ -0,0 +1,545 @@
+import test from 'tape-catch';
+import vtkMath from 'vtk.js/Sources/Common/Core/Math';
+import vtkPolygon from 'vtk.js/Sources/Common/DataModel/Polygon';
+import {
+  PolygonWithPointIntersectionState,
+  PolygonWithCellIntersectionState,
+  EPSILON,
+} from 'vtk.js/Sources/Common/DataModel/Polygon/Constants';
+
+test('Test vtkPolygon instance', (t) => {
+  t.ok(vtkPolygon, 'Make sure the class definition exists');
+  const instance = vtkPolygon.newInstance();
+  t.ok(instance);
+  t.end();
+});
+
+test('Test vtkPolygon computeNormal', (t) => {
+  const polygon = vtkPolygon.newInstance();
+
+  const poly = [
+    [0, 0, 0],
+    [3, 0, 0],
+    [3, 3, 3],
+    [0, 3, 3],
+  ];
+  polygon.setPoints(poly);
+
+  t.ok(
+    polygon.computeNormal() - 25.455844 < EPSILON &&
+      vtkMath.areEquals(
+        polygon.get().normal,
+        [0, -0.7071067811865476, 0.7071067811865476]
+      )
+  );
+
+  t.end();
+});
+
+test('Test vtkPolygon computeCentroid', (t) => {
+  const polygon = vtkPolygon.newInstance();
+
+  let poly = [
+    [0, 0, 0],
+    [3, 0, 0],
+    [3, 3, 3],
+    [0, 3, 3],
+  ];
+  polygon.setPoints(poly);
+
+  t.ok(
+    vtkMath.areEquals([1.5, 1.5, 1.5], polygon.computeCentroid(), EPSILON),
+    'centroid should be [1.5, 1.5, 1.5]'
+  );
+
+  // Same points but in a different order to test degenerate polygon
+  poly = [
+    [0, 0, 0],
+    [3, 0, 0],
+    [0, 3, 3],
+    [3, 3, 3],
+  ];
+  polygon.setPoints(poly);
+
+  t.ok(
+    polygon.computeCentroid() === null,
+    'giving degenerate polygon returns null'
+  );
+
+  t.end();
+});
+
+test('Test vtkPolygon computeArea', (t) => {
+  const polygon = vtkPolygon.newInstance();
+
+  polygon.setPoints([
+    [0, 10, 0],
+    [0, 0, 1],
+  ]);
+
+  t.ok(
+    polygon.computeArea() < EPSILON,
+    'Area of polygon with less than 3 points should be 0'
+  );
+
+  let poly = [
+    [0, 0, 0],
+    [3, 0, 0],
+    [3, 3, 3],
+    [0, 3, 3],
+  ];
+  polygon.setPoints(poly);
+
+  t.ok(
+    polygon.computeArea() - 3 * 3 * 2 * Math.sqrt(2) < EPSILON,
+    'Area should be 3*3*2*sqrt(2)'
+  );
+
+  // Same points but in a different order to test degenerate polygon
+  poly = [
+    [0, 0, 0],
+    [3, 0, 0],
+    [0, 3, 3],
+    [3, 3, 3],
+  ];
+  polygon.setPoints(poly);
+  t.ok(Number.isNaN(polygon.computeArea()), 'Area of degenerate area is NaN');
+
+  t.end();
+});
+
+test('Test vtkPolygon distanceToPolygon', (t) => {
+  const polygon = vtkPolygon.newInstance();
+
+  const poly = [
+    [0, 0, 0],
+    [3, 0, 0],
+    [3, 3, 3],
+    [0, 3, 3],
+  ];
+  polygon.setPoints(poly);
+
+  const closest = [0, 0, 0];
+  let distToPolygon = polygon.distanceToPolygon([1.5, 1.5, 1.5], closest);
+  t.ok(
+    distToPolygon.t - Number.MIN_VALUE < EPSILON &&
+      distToPolygon.distance === 0 &&
+      vtkMath.areEquals([1.5, 1.5, 1.5], closest),
+    'Inside Polygon'
+  );
+
+  distToPolygon = polygon.distanceToPolygon([1.5, 0, -3], closest);
+  t.ok(
+    distToPolygon.t - 1.5 < EPSILON &&
+      distToPolygon.distance - 3 < EPSILON &&
+      vtkMath.areEquals([1.5, 0, 0], closest), // TBD
+    'Point outside polygon'
+  );
+
+  t.end();
+});
+
+test('Test vtkPolygon triangulate', (t) => {
+  const polygon = vtkPolygon.newInstance();
+
+  let poly = [
+    [0, 0, 0],
+    [3, 0, 0],
+    [3, 3, 3],
+    [0, 3, 3],
+  ];
+  polygon.setPoints(poly);
+  let triangulation = polygon.triangulate();
+  let expected = [0, 0, 0, 3, 0, 0, 0, 3, 3, 3, 3, 3, 0, 3, 3, 3, 0, 0];
+  t.ok(vtkMath.areEquals(triangulation, expected), 'Polygon with 4 points');
+  poly = [
+    [0, 0, 0],
+    [1, 2, 1],
+    [2, 0, 0],
+    [1, 3, 1],
+  ];
+  polygon.setPoints(poly);
+
+  triangulation = polygon.triangulate();
+  expected = [0, 0, 0, 1, 2, 1, 1, 3, 1, 2, 0, 0, 1, 3, 1, 1, 2, 1];
+
+  t.ok(
+    vtkMath.areEquals(triangulation, expected),
+    'Can triangulate non convex polygon'
+  );
+
+  poly = [
+    [1, 0, 0],
+    [2, 0, 0],
+    [3, 1, 0],
+    [3, 2, 0],
+    [2, 3, 0],
+    [1, 3, 0],
+    [0, 2, 0],
+    [0, 1, 0],
+  ];
+  polygon.setPoints(poly);
+
+  triangulation = polygon.triangulate();
+  // prettier-ignore
+  expected = [
+    1, 0, 0,
+    2, 0, 0,
+    0, 1, 0,
+    0, 2, 0,
+    0, 1, 0,
+    1, 3, 0,
+    2, 3, 0,
+    1, 3, 0,
+    3, 2, 0,
+    3, 1, 0,
+    3, 2, 0,
+    2, 0, 0,
+    1, 3, 0,
+    0, 1, 0,
+    3, 2, 0,
+    2, 0, 0,
+    3, 2, 0,
+    0, 1, 0,
+  ];
+  t.ok(vtkMath.areEquals(triangulation, expected), 'Can triangulate');
+
+  poly = [
+    [0, 0, 0],
+    [3, 0, 0],
+    [0, 3, 3],
+    [3, 3, 3],
+  ];
+  polygon.setPoints(poly);
+
+  triangulation = polygon.triangulate();
+  t.ok(triangulation === null, 'Cannot triangulate degenerate polygon');
+
+  poly = [
+    [0, 1, 2],
+    [3, 4, 5],
+  ];
+  polygon.setPoints(poly);
+  triangulation = polygon.triangulate();
+  t.ok(triangulation === null, 'Cannot triangulate polygon with 2 points');
+
+  t.end();
+});
+
+test('Test vtkPolygon pointInPolygon', (t) => {
+  const polygon = vtkPolygon.newInstance();
+  let poly = [
+    [0, 0, 0],
+    [3, 0, 0],
+    [3, 3, 3],
+    [0, 3, 3],
+  ];
+  polygon.setPoints(poly);
+
+  t.ok(
+    polygon.pointInPolygon([1.5, 1.5, 1.5]) ===
+      PolygonWithPointIntersectionState.INSIDE,
+    '[1.5,1.5,1.5] is in polygon'
+  );
+  t.ok(
+    polygon.pointInPolygon([5, 5, 5]) ===
+      PolygonWithPointIntersectionState.OUTSIDE,
+    '[5,5,5] is not in polygon'
+  );
+  t.ok(
+    polygon.pointInPolygon([1.5, 0, 0]) ===
+      PolygonWithPointIntersectionState.INSIDE,
+    '[1.5,0,0] is on edge of polygon (considered inside)'
+  );
+  // Same points but in a different order to test degenerate polygon
+  poly = [
+    [0, 0, 0],
+    [3, 0, 0],
+    [0, 3, 3],
+    [3, 3, 3],
+  ];
+  polygon.setPoints(poly);
+
+  t.end();
+});
+
+test('Test vtkPolygon isConvex', (t) => {
+  const polygon = vtkPolygon.newInstance();
+
+  let poly = [
+    [0, 0, 0],
+    [3, 0, 0],
+    [3, 3, 3],
+    [0, 3, 3],
+  ];
+  polygon.setPoints(poly);
+
+  t.ok(polygon.isConvex(), 'convex polygon');
+  poly = [
+    [0, 0, 0],
+    [1, 2, 1],
+    [2, 0, 0],
+    [1, 3, 1],
+  ];
+  polygon.setPoints(poly);
+
+  t.ok(!polygon.isConvex(), 'non convex polygon');
+
+  // Same points but in a different order to test degenerate polygon
+  poly = [
+    [0, 0, 0],
+    [3, 0, 0],
+    [0, 3, 3],
+    [3, 3, 3],
+  ];
+  polygon.setPoints(poly);
+  t.ok(!polygon.isConvex(), 'degenerate polygon');
+
+  t.end();
+});
+
+test('Test vtkPolygon interpolateFunctions', (t) => {
+  const polygon = vtkPolygon.newInstance();
+
+  let poly = [
+    [0, 0, 0],
+    [3, 0, 0],
+    [3, 3, 3],
+    [0, 3, 3],
+  ];
+  polygon.setPoints(poly);
+
+  t.ok(
+    vtkMath.areEquals(
+      polygon.interpolateFunctions([1.5, 1.5, 1.5], false),
+      [0.25, 0.25, 0.25, 0.25]
+    ),
+    'convex polygon'
+  );
+
+  // Test for MVC method
+  t.ok(
+    vtkMath.areEquals(
+      polygon.interpolateFunctions([1.5, 1.5, 1.5], true),
+      [0.25, 0.25, 0.25, 0.25]
+    ),
+    'convex polygon MVC'
+  );
+
+  t.ok(
+    vtkMath.areEquals(
+      polygon.interpolateFunctions(
+        [3 + EPSILON, 3 + EPSILON, 3 + EPSILON],
+        false
+      ),
+      [0, 0, 1, 0]
+    ),
+    'interpolate point close to one of the polygon point'
+  );
+
+  t.ok(
+    vtkMath.areEquals(
+      polygon.interpolateFunctions(
+        [
+          3 + (1 / Math.sqrt(3)) * EPSILON,
+          3 + (1 / Math.sqrt(3)) * EPSILON,
+          3 + (1 / Math.sqrt(3)) * EPSILON,
+        ],
+        true
+      ),
+      [0, 0, 1, 0]
+    ),
+    'interpolate point close to one of the point of the polygon, MVC interpolation'
+  );
+
+  // Same points but in a different order to test degenerate polygon
+  poly = [
+    [0, 0, 0],
+    [3, 0, 0],
+    [0, 3, 3],
+    [3, 3, 3],
+  ];
+  polygon.setPoints(poly);
+  t.ok(
+    vtkMath.areEquals(
+      polygon.interpolateFunctions([1.5, 1.5, 1.5], false),
+      [0.25, 0.25, 0.25, 0.25]
+    ),
+    'degenerate polygon'
+  );
+
+  t.end();
+});
+
+test('Test vtkPolygon intersectWithLine', (t) => {
+  const polygon = vtkPolygon.newInstance();
+  const poly = [
+    [0, 0, 0],
+    [3, 0, 0],
+    [3, 3, 3],
+    [0, 3, 3],
+  ];
+  polygon.setPoints(poly);
+
+  let p0 = [1.5, 1.5, 0];
+  let p1 = [1.5, 1.5, 5];
+  let intersect = polygon.intersectWithLine(p0, p1);
+  vtkMath.roundVector(intersect.x, intersect.x, 6);
+  t.ok(
+    intersect.intersection &&
+      intersect.betweenPoints &&
+      intersect.t - 0.3 <= EPSILON &&
+      vtkMath.areEquals(intersect.x, [1.5, 1.5, 1.5]),
+    'Point in polygon'
+  );
+
+  p0 = [3, 1.5, 0];
+  p1 = [3, 1.5, 3];
+  intersect = polygon.intersectWithLine(p0, p1);
+  t.ok(
+    intersect.intersection &&
+      intersect.betweenPoints &&
+      vtkMath.areEquals(intersect.x, [3, 1.5, 1.5]),
+    'Point on edge of polygon'
+  );
+
+  p0 = [10, 1, 1];
+  p1 = [10, 6, 6];
+  intersect = polygon.intersectWithLine(p0, p1);
+  t.ok(
+    !intersect.intersection &&
+      !intersect.betweenPoints &&
+      intersect.t >= Number.MAX_VALUE &&
+      intersect.x.length === 0,
+    'Line parallele to polygon'
+  );
+
+  p0 = [5, 5, 0];
+  p1 = [5, 5, 5];
+  intersect = polygon.intersectWithLine(p0, p1);
+  vtkMath.roundVector(intersect.x, intersect.x, 6);
+  console.log(intersect);
+  t.ok(
+    !intersect.intersection &&
+      intersect.betweenPoints &&
+      intersect.t - 1 < EPSILON &&
+      vtkMath.areEquals(intersect.x, [5, 5, 5]),
+    'Intersect plane but not polygon'
+  );
+
+  t.end();
+});
+
+test('Test vtkPolygon intersectConvex2DCells', (t) => {
+  const polygon = vtkPolygon.newInstance();
+  const poly = [
+    [0, 0, 0],
+    [3, 0, 0],
+    [3, 3, 3],
+    [0, 3, 3],
+  ];
+  polygon.setPoints(poly);
+
+  const polygonToIntersect = vtkPolygon.newInstance();
+  let polyToIntersect = [
+    [0, 1.5, 0],
+    [3, 1.5, 0],
+    [3, 1.5, 1.5],
+    [0, 1.5, 3],
+  ];
+  polygonToIntersect.setPoints(polyToIntersect);
+
+  let intersection = polygon.intersectConvex2DCells(polygonToIntersect);
+  t.ok(
+    intersection.intersection ===
+      PolygonWithCellIntersectionState.LINE_INTERSECTION &&
+      ((vtkMath.areEquals(intersection.x0, [0, 1.5, 1.5]) &&
+        vtkMath.areEquals(intersection.x1, [3, 1.5, 1.5])) ||
+        (vtkMath.areEquals(intersection.x0, [3, 1.5, 1.5]) &&
+          vtkMath.areEquals(intersection.x1, [0, 1.5, 1.5]))),
+    'line intersection'
+  );
+
+  polyToIntersect = [
+    [0, 1.5, 1.5],
+    [3, 1.5, 1.5],
+    [3, 1.5, 1.5],
+    [0, 1.5, 3],
+  ];
+  polygonToIntersect.setPoints(polyToIntersect);
+
+  intersection = polygon.intersectConvex2DCells(polygonToIntersect);
+  console.log(intersection);
+  t.ok(
+    intersection.intersection ===
+      PolygonWithCellIntersectionState.LINE_INTERSECTION &&
+      ((vtkMath.areEquals(intersection.x0, [0, 1.5, 1.5]) &&
+        vtkMath.areEquals(intersection.x1, [3, 1.5, 1.5])) ||
+        (vtkMath.areEquals(intersection.x0, [3, 1.5, 1.5]) &&
+          vtkMath.areEquals(intersection.x1, [0, 1.5, 1.5]))),
+    'line intersection on edge of one polygon'
+  );
+
+  polyToIntersect = [
+    [1, 1.5, 1.5],
+    [3, 1.5, 3],
+    [0, 1.5, 3],
+  ];
+  polygonToIntersect.setPoints(polyToIntersect);
+
+  intersection = polygon.intersectConvex2DCells(polygonToIntersect);
+  t.ok(
+    intersection.intersection ===
+      PolygonWithCellIntersectionState.POINT_INTERSECTION &&
+      vtkMath.areEquals(intersection.x0, [1, 1.5, 1.5]) &&
+      vtkMath.areEquals(intersection.x1, []),
+    'point intersection'
+  );
+
+  polyToIntersect = [
+    [7, 7, 7],
+    [4, 4, 4],
+    [10, 10, 10],
+  ];
+  polygonToIntersect.setPoints(polyToIntersect);
+
+  intersection = polygon.intersectConvex2DCells(polygonToIntersect);
+  t.ok(
+    intersection.intersection ===
+      PolygonWithCellIntersectionState.NO_INTERSECTION &&
+      vtkMath.areEquals(intersection.x0, []) &&
+      vtkMath.areEquals(intersection.x1, []),
+    'no intersection'
+  );
+
+  polyToIntersect = [
+    [1, 1, 1],
+    [2, 1, 1],
+    [2, 2, 2],
+    [1, 2, 2],
+  ];
+  polygonToIntersect.setPoints(polyToIntersect);
+  intersection = polygon.intersectConvex2DCells(polygonToIntersect);
+  t.ok(
+    intersection.intersection === PolygonWithCellIntersectionState.INCLUDED &&
+      vtkMath.areEquals(intersection.x0, []) &&
+      vtkMath.areEquals(intersection.x1, []),
+    'coincident polygons'
+  );
+
+  polyToIntersect = [
+    [2, 1, 1],
+    [4, 1, 1],
+    [4, 2, 2],
+    [2, 2, 2],
+  ];
+  polygonToIntersect.setPoints(polyToIntersect);
+  intersection = polygon.intersectConvex2DCells(polygonToIntersect);
+  t.ok(
+    intersection.intersection === PolygonWithCellIntersectionState.OVERLAP &&
+      vtkMath.areEquals(intersection.x0, []) &&
+      vtkMath.areEquals(intersection.x1, []),
+    'overlaping polygons'
+  );
+  t.end();
+});
diff --git a/Sources/Filters/General/ContourTriangulator/helper.js b/Sources/Filters/General/ContourTriangulator/helper.js
index 91a9d8f971b..5e0f07b9bd7 100644
--- a/Sources/Filters/General/ContourTriangulator/helper.js
+++ b/Sources/Filters/General/ContourTriangulator/helper.js
@@ -775,9 +775,7 @@ export function vtkCCSSplitAtPinchPoints(
     n = poly.length;
 
     bounds = [];
-    tol =
-      CCS_POLYGON_TOLERANCE *
-      Math.sqrt(vtkPolygon.getBounds(poly, points, bounds));
+    tol = CCS_POLYGON_TOLERANCE * Math.sqrt(vtkPolygon.getBounds(points));
 
     if (tol === 0) {
       // eslint-disable-next-line no-continue
@@ -929,8 +927,7 @@ export function vtkCCSFindTrueEdges(polys, points, polyEdges, originalEdges) {
     let m = n;
 
     // Compute bounds for tolerance
-    const bounds = [];
-    const tol2 = vtkPolygon.getBounds(oldPoly, points, bounds) * atol2;
+    const tol2 = vtkPolygon.getBounds(points) * atol2;
 
     // The new poly
     const newPoly = [];
@@ -1230,7 +1227,7 @@ export function vtkCCSPrepareForPolyInPoly(outerPoly, points, pp, bounds) {
 
   // Find the bounding box and tolerance for the polygon
   return (
-    vtkPolygon.getBounds(outerPoly, points, bounds) *
+    vtkPolygon.getBounds(points) *
     (CCS_POLYGON_TOLERANCE * CCS_POLYGON_TOLERANCE)
   );
 }
diff --git a/Sources/Filters/General/ContourTriangulator/index.js b/Sources/Filters/General/ContourTriangulator/index.js
index 8414ef17a9e..cdddc55412b 100644
--- a/Sources/Filters/General/ContourTriangulator/index.js
+++ b/Sources/Filters/General/ContourTriangulator/index.js
@@ -60,7 +60,7 @@ function triangulateContours(
     let maxnorm = 0;
     const n = [];
     for (let i = 0; i < newPolys.length; i++) {
-      const norm = vtkPolygon.getNormal(newPolys[i], points, n);
+      const norm = vtkPolygon.getNormal(newPolys[i], points, n) ** 2;
       if (norm > maxnorm) {
         maxnorm = norm;
         computedNormal[0] = n[0];
@@ -186,7 +186,7 @@ function triangulatePolygon(polygon, points, triangles) {
 
   let success = true;
   const normal = [];
-  const norm = vtkPolygon.getNormal(poly, points, normal);
+  const norm = vtkPolygon.getNormal(poly, points, normal) ** 2;
   if (norm !== 0) {
     success = vtkCCSTriangulate(
       poly,
diff --git a/Sources/Filters/General/PaintFilter/index.js b/Sources/Filters/General/PaintFilter/index.js
index 9c772df3f3e..6b766e55fe1 100644
--- a/Sources/Filters/General/PaintFilter/index.js
+++ b/Sources/Filters/General/PaintFilter/index.js
@@ -203,12 +203,16 @@ function vtkPaintFilter(publicAPI, model) {
       }
       polygon.setPoints(poly);
 
-      if (!polygon.triangulate()) {
+      const points = polygon.triangulate();
+      const triangulationSuceeded = points !== null;
+      if (triangulationSuceeded) {
         console.log('triangulation failed!');
       }
 
-      const points = polygon.getPointArray();
-      const triangleList = new Float32Array(points.length);
+      const triangleList = triangulationSuceeded
+        ? new Float32Array(points.length)
+        : new Float32Array(0);
+
       const numPoints = Math.floor(triangleList.length / 3);
       for (let i = 0; i < numPoints; i++) {
         const point = points.slice(3 * i, 3 * i + 3);
diff --git a/Sources/Filters/General/TriangleFilter/index.js b/Sources/Filters/General/TriangleFilter/index.js
index 0be51816142..e9a8834a119 100644
--- a/Sources/Filters/General/TriangleFilter/index.js
+++ b/Sources/Filters/General/TriangleFilter/index.js
@@ -56,13 +56,16 @@ function vtkTriangleFilter(publicAPI, model) {
           const polygon = vtkPolygon.newInstance();
           polygon.setPoints(cellPoints);
 
-          if (!polygon.triangulate()) {
+          const newCellPoints = polygon.triangulate();
+          const triangulationSucceeded = newCellPoints !== null;
+          if (!triangulationSucceeded) {
             vtkWarningMacro(`Triangulation failed at cellOffset ${c}`);
             ++model.errorCount;
           }
 
-          const newCellPoints = polygon.getPointArray();
-          const numSimplices = Math.floor(newCellPoints.length / 9);
+          const numSimplices = triangulationSucceeded
+            ? Math.floor(newCellPoints.length / 9)
+            : 0;
           const triPts = [];
           triPts.length = 9;
           for (let i = 0; i < numSimplices; i++) {

From bb6c1994482d999ec1408cfa192d0dc2fcc7c954 Mon Sep 17 00:00:00 2001
From: joannacirillo <cirillojoanna@gmail.com>
Date: Wed, 31 Aug 2022 17:14:56 +0200
Subject: [PATCH 2/2] fix(polygon): implement polygonWithPolygonIntersection

WORK IN PROGRESS
---
 Sources/Common/DataModel/Polygon/index.js  | 94 +++++++++++++++++++++-
 Sources/Common/DataModel/Triangle/index.js | 59 ++++++++++++++
 2 files changed, 152 insertions(+), 1 deletion(-)

diff --git a/Sources/Common/DataModel/Polygon/index.js b/Sources/Common/DataModel/Polygon/index.js
index f01de7aaf83..2267edff146 100644
--- a/Sources/Common/DataModel/Polygon/index.js
+++ b/Sources/Common/DataModel/Polygon/index.js
@@ -5,6 +5,7 @@ import vtkPlane from 'vtk.js/Sources/Common/DataModel/Plane';
 import vtkPriorityQueue from 'vtk.js/Sources/Common/Core/PriorityQueue';
 import { IntersectionState } from 'vtk.js/Sources/Common/DataModel/Line/Constants';
 import vtkBoundingBox from 'vtk.js/Sources/Common/DataModel/BoundingBox';
+import vtkTriangle from 'vtk.js/Sources/Common/DataModel/Triangle';
 import vtkCell from 'vtk.js/Sources/Common/DataModel/Cell';
 import vtkPoints from 'vtk.js/Sources/Common/Core/Points';
 
@@ -394,7 +395,7 @@ function computeArea(points, normal) {
 }
 
 //------------------------------------------------------------------------------
-function distanceToPolygon(x, points, closest) {
+export function distanceToPolygon(x, points, closest) {
   const outObj = { t: Number.MIN_VALUE, distance: 0 };
   const bounds = points.getBounds();
   // First check to see if the point is inside the polygon
@@ -444,6 +445,97 @@ function distanceToPolygon(x, points, closest) {
   return outObj;
 }
 
+//------------------------------------------------------------------------------
+// Method intersects two polygons. You must supply the number of points and
+// point coordinates (npts, *pts) and the bounding box (bounds) of the two
+// polygons. Also supply a tolerance squared for controlling
+// error. The method returns 1 if there is an intersection, and 0 if
+// not. A single point of intersection x[3] is also returned if there
+// is an intersection.
+function intersectPolygonWithPolygon(points1, points2, x) {
+  const n = [0, 0, 0];
+  const coords = [0, 0, 0];
+  const p1 = [0, 0, 0];
+  const p2 = [0, 0, 0];
+  const ray = [0, 0, 0];
+  let t;
+
+  //  Intersect each edge of first polygon against second
+  //
+  getNormal(points2, n);
+
+  const npts = points1.getNumberOfPoints();
+  const npts2 = points2.getNumberOfPoints();
+
+  const bounds2 = points2.getBounds();
+  for (let i = 0; i < npts; i++) {
+    points1.getPoint(3 * i, p1);
+    points1.getPoint((i + 1) % npts, p2);
+
+    for (let j = 0; j < 3; j++) {
+      ray[j] = p2[j] - p1[j];
+    }
+    if (!vtkBoundingBox.intersectBox(bounds2, p1, ray, coords, t)) {
+      // eslint-disable-next-line no-continue
+      continue;
+    }
+
+    if (vtkPlane.intersectWithLine(p1, p2, x, n) === 1) {
+      if (
+        vtkTriangle.pointInTriangle(
+          x,
+          [...points2.getPoint(0)],
+          [...points2.getPoint(3)],
+          [...points2.getPoint(6)]
+        ) ||
+        pointInPolygon(x, points2, bounds2, n) ===
+          PolygonWithPointIntersectionState.INSIDE
+      ) {
+        return true;
+      }
+    } else {
+      return false;
+    }
+  }
+
+  //  Intersect each edge of second polygon against first
+  //
+  getNormal(points1, n);
+
+  const bounds = points1.getBounds();
+  for (i = 0; i < npts2; i++) {
+    points2.getPoint(3 * i, p1);
+    points2.getPoint((i + 1) % npts, p2);
+
+    for (j = 0; j < 3; j++) {
+      ray[j] = p2[j] - p1[j];
+    }
+
+    if (!vtkBoundingBoxBox.intersectBox(bounds, p1, ray, coords, t)) {
+      continue;
+    }
+
+    if (vtkPlane.intersectWithLine(p1, p2, x, n).intersection) {
+      if (
+        vtkTriangle.pointInTriangle(
+          x,
+          [...points1.getPoint(0)],
+          [...points1.getPoint(3)],
+          [...points1.getPoint(6)]
+        ) ||
+        pointInPolygon(x, points1, bounds, n) ===
+          PolygonWithPointIntersectionState.INSIDE
+      ) {
+        return 1;
+      }
+    } else {
+      return 0;
+    }
+  }
+
+  return 0;
+}
+
 //------------------------------------------------------------------------------
 function arePolygonsCoincident(points1, points2) {
   // Check if the planes defined by the polygons points are coincident
diff --git a/Sources/Common/DataModel/Triangle/index.js b/Sources/Common/DataModel/Triangle/index.js
index fb6e0b3080e..ea926f12553 100644
--- a/Sources/Common/DataModel/Triangle/index.js
+++ b/Sources/Common/DataModel/Triangle/index.js
@@ -4,6 +4,8 @@ import * as vtkMath from 'vtk.js/Sources/Common/Core/Math';
 import vtkLine from 'vtk.js/Sources/Common/DataModel/Line';
 import vtkPlane from 'vtk.js/Sources/Common/DataModel/Plane';
 
+const TOLERANCE = 1e-6;
+
 // ----------------------------------------------------------------------------
 // Global methods
 // ----------------------------------------------------------------------------
@@ -221,6 +223,62 @@ function intersectWithTriangle(p1, q1, r1, p2, q2, r2, tolerance = 1e-6) {
   return { intersect: true, coplanar, pt1, pt2, surfaceId };
 }
 
+//------------------------------------------------------------------------------
+// Given a point x, determine whether it is inside (within the
+// a tolerance squared) the triangle defined by the three
+// coordinate values p1, p2, p3. Method is via comparing dot products.
+// (Note: in current implementation the tolerance only works in the
+// neighborhood of the three vertices of the triangle.
+function pointInTriangle(x, p1, p2, p3) {
+  const tol2 = TOLERANCE * TOLERANCE;
+  const x1 = [0, 0, 0];
+  const x2 = [0, 0, 0];
+  const x3 = [0, 0, 0];
+  const v13 = [0, 0, 0];
+  const v21 = [0, 0, 0];
+  const v32 = [0, 0, 0];
+  const n1 = [0, 0, 0];
+  const n2 = [0, 0, 0];
+  const n3 = [0, 0, 0];
+
+  //  Compute appropriate vectors
+  //
+  for (let i = 0; i < 3; i++) {
+    x1[i] = x[i] - p1[i];
+    x2[i] = x[i] - p2[i];
+    x3[i] = x[i] - p3[i];
+    v13[i] = p1[i] - p3[i];
+    v21[i] = p2[i] - p1[i];
+    v32[i] = p3[i] - p2[i];
+  }
+
+  //  See whether intersection point is within tolerance of a vertex.
+  //
+  if (
+    x1[0] * x1[0] + x1[1] * x1[1] + x1[2] * x1[2] <= tol2 ||
+    x2[0] * x2[0] + x2[1] * x2[1] + x2[2] * x2[2] <= tol2 ||
+    x3[0] * x3[0] + x3[1] * x3[1] + x3[2] * x3[2] <= tol2
+  ) {
+    return true;
+  }
+
+  //  If not near a vertex, check whether point is inside of triangular face.
+  //
+  //  Obtain normal off of triangular face
+  //
+  vtkMath.cross(x1, v13, n1);
+  vtkMath.cross(x2, v21, n2);
+  vtkMath.cross(x3, v32, n3);
+
+  //  Check whether ALL the three normals go in same direction
+  //
+  return (
+    vtkMath.dot(n1, n2) >= 0.0 &&
+    vtkMath.dot(n2, n3) >= 0.0 &&
+    vtkMath.dot(n1, n3) >= 0.0
+  );
+}
+
 // ----------------------------------------------------------------------------
 // Static API
 // ----------------------------------------------------------------------------
@@ -229,6 +287,7 @@ export const STATIC = {
   computeNormalDirection,
   computeNormal,
   intersectWithTriangle,
+  pointInTriangle,
 };
 
 // ----------------------------------------------------------------------------