ContourTriangulator

Introduction

vtkContourTriangulator

Methods

extend

Method use to decorate a given object (publicAPI+model) with vtkContourTriangulator characteristics.

Argument Type Required Description
publicAPI Yes object on which methods will be bounds (public)
model Yes object on which data structure will be bounds (protected)
initialValues IContourTriangulatorInitialValues No (default: {})

getTriangulatePolys

Returns the behavior of the filter regarding polys.

newInstance

Method use to create a new instance of vtkContourTriangulator

Argument Type Required Description
initialValues IContourTriangulatorInitialValues No for pre-setting some of its content

requestData

Argument Type Required Description
inData any Yes
outData any Yes

setTriangulatePolys

Sets the behavior of the filter regarding polys.

Argument Type Required Description
triangulate boolean Yes whether the filter should triangulate polys or leave them untouched. True by default

triangulateContours

This is a complex subroutine that takes a collection of lines that
were formed by cutting a polydata with a plane, and generates
a face that has those lines as its edges. The lines must form one
or more closed contours, but they need not be sorted.

Only “numLine” lines starting from “firstLine” are used to create new
polygons, and the new polygons are appended to “polys”. The normal of
the cut plane must be provided so that polys will be correctly oriented.

Given some closed contour lines, create a triangle mesh that fills
those lines. The input lines do not have to be in tail-to-tip order.
Only numLines starting from firstLine will be used. Note that holes
can be indicated by contour loops whose normals are in the opposite
direction to the provided normal.

Argument Type Required Description
polyData vtkPolyData Yes
firstLine Number Yes
numLines Number Yes
polys vtkCellArray Yes
normal Nullable. Yes If null, the function will compute the normal of the polys.
triangulatePolys Boolean No (default: true) If set to true the resulting polygons will be triangulated, otherwise the polygons
themselves will be added to the output.
diagnoseOnTriangulationError Boolean No (default: false) If this option is set to true and there was a triangulation error
this will add the polys as outlines to the output.

Returns

Type Description
Boolean Returns true if triangulation was successful, false otherwise.

triangulatePolygon

A robust method for triangulating a polygon. It cleans up the polygon
and then applies the ear-cut triangulation. A zero return value
indicates that triangulation failed.

Argument Type Required Description
polygon Array. or TypedArray Yes Array of point indices defining the polygon
points vtkPoints Yes The point coordinates of the polygon
triangles vtkCellArray Yes The cell array that is going to be filled with the triangulation

Returns

Type Description
Boolean Returns true if triangulation was successful, false otherwise.

Source

Constants.js
export const CCS_POLYGON_TOLERANCE = 1e-5;

export default { CCS_POLYGON_TOLERANCE };
helper.js
import macro from 'vtk.js/Sources/macros';
import vtkPoints from 'vtk.js/Sources/Common/Core/Points';
import * as vtkMath from 'vtk.js/Sources/Common/Core/Math';
import vtkLine from 'vtk.js/Sources/Common/DataModel/Line';
import vtkPolygon from 'vtk.js/Sources/Common/DataModel/Polygon';
import vtkIncrementalOctreePointLocator from 'vtk.js/Sources/Common/DataModel/IncrementalOctreePointLocator';
import { VtkDataTypes } from 'vtk.js/Sources/Common/Core/DataArray/Constants';
import { CCS_POLYGON_TOLERANCE } from './Constants';
import { PolygonWithPointIntersectionState } from '../../../Common/DataModel/Polygon/Constants';

const { vtkErrorMacro } = macro;

/**
* Reverse the elements between the indices firstIdx and lastIdx of the given array arr.
*
* @param {Array|TypedArray} arr
* @param {Number} firstIdx
* @param {Number} lastIdx
*/
export function reverseElements(
arr,
firstIdx = undefined,
lastIdx = undefined
) {
const first = firstIdx ?? 0;
const last = lastIdx ?? arr.length - 1;
const mid = first + Math.floor((last - first) / 2);
for (let i = first; i <= mid; ++i) {
[arr[i], arr[last - (i - first)]] = [arr[last - (i - first)], arr[i]];
}
}

// ---------------------------------------------------
/**
* Compute the quality of a triangle.
*
* @param {Vector3} p0
* @param {Vector3} p1
* @param {Vector3} p2
* @param {Vector3} normal
* @returns {Number}
*/
export function vtkCCSTriangleQuality(p0, p1, p2, normal) {
const u = [];
const v = [];
const w = [];

vtkMath.subtract(p1, p0, u);
vtkMath.subtract(p2, p1, v);
vtkMath.subtract(p0, p2, w);

const area2 =
(u[1] * v[2] - u[2] * v[1]) * normal[0] +
(u[2] * v[0] - u[0] * v[2]) * normal[1] +
(u[0] * v[1] - u[1] * v[0]) * normal[2];

let perim =
Math.sqrt(u[0] * u[0] + u[1] * u[1] + u[2] * u[2]) +
Math.sqrt(v[0] * v[0] + v[1] * v[1] + v[2] * v[2]) +
Math.sqrt(w[0] * w[0] + w[1] * w[1] + w[2] * w[2]);

perim *= perim; // square the perimeter
perim = perim !== 0 ? perim : 1.0;

// use a normalization factor so equilateral quality is 1.0
return (area2 / perim) * 10.392304845413264;
}

// ---------------------------------------------------
/**
* Insert a triangle, and subdivide that triangle if one of
* its edges originally had more than two points before
* vtkCCSFindTrueEdges was called. Is called by vtkCCSTriangulate.
*
* @param {vtkCellArray} polys
* @param {Array|TypedArray} poly
* @param {Vector3} trids
* @param {Array|TypedArray} polyEdges
* @param {Array|TypedArray} originalEdges
*/
export function vtkCCSInsertTriangle(
polys,
poly,
trids,
polyEdges,
originalEdges
) {
const nextVert = [1, 2, 0];

// To store how many of originalEdges match
let edgeCount = 0;
const edgeLocs = [-1, -1, -1];

// Check for original edge matches
for (let vert = 0; vert < 3; vert++) {
const currId = trids[vert];
const edgeLoc = polyEdges[currId];
if (edgeLoc >= 0) {
let nextId = currId + 1;
if (nextId === poly.length) {
nextId = 0;
}

// Is the triangle edge a polygon edge?
if (nextId === trids[nextVert[vert]]) {
edgeLocs[vert] = edgeLoc;
edgeCount++;
}
}
}

if (edgeCount === 0) {
// No special edge handling, so just do one triangle
polys.insertNextCell([poly[trids[0]], poly[trids[1]], poly[trids[2]]]);
} else {
// Make triangle fans for edges with extra points
const edgePts = [
[poly[trids[0]], poly[trids[1]]],
[poly[trids[1]], poly[trids[2]]],
[poly[trids[2]], poly[trids[0]]],
];

// Find out which edge has the most extra points
let maxPoints = 0;
let currSide = 0;
for (let i = 0; i < 3; i++) {
if (edgeLocs[i] >= 0) {
const edgeLoc = edgeLocs[i];
const npts = originalEdges[edgeLoc];
const pts = originalEdges.slice(edgeLoc + 1, edgeLoc + 1 + npts);
if (!(edgePts[i][0] === pts[0] || edgePts[i][1] === pts[npts - 1])) {
vtkErrorMacro('assertion error in vtkCCSInsertTriangle');
}
if (npts > maxPoints) {
maxPoints = npts;
currSide = i;
}
edgePts[i] = pts;
}
}

// Find the edges before/after the edge with most points
const prevSide = (currSide + 2) % 3;
const nextSide = (currSide + 1) % 3;

// If other edges have only 2 points, nothing to do with them
const prevNeeded = edgePts[prevSide].length > 2;
const nextNeeded = edgePts[nextSide].length > 2;

// The tail is the common point in the triangle fan
const tailPtIds = [];
tailPtIds[prevSide] = edgePts[currSide][1];
tailPtIds[currSide] = edgePts[prevSide][0];
tailPtIds[nextSide] = edgePts[currSide][edgePts[currSide].length - 2];

// Go through the sides and make the fans
for (let side = 0; side < 3; side++) {
if (
(side !== prevSide || prevNeeded) &&
(side !== nextSide || nextNeeded)
) {
let m = 0;
let n = edgePts[side].length - 1;

if (side === currSide) {
m += prevNeeded;
n -= nextNeeded;
}

for (let k = m; k < n; k++) {
polys.insertNextCell([
edgePts[side][k],
edgePts[side][k + 1],
tailPtIds[side],
]);
}
}
}
}
}

// ---------------------------------------------------
/**
* Triangulate a polygon that has been simplified by FindTrueEdges.
* This will re-insert the original edges. The output triangles are
* appended to "polys" and, for each stored triangle, "color" will
* be added to "scalars". The final two arguments (polygon and
* triangles) are only for temporary storage.
* The return value is true if triangulation was successful.
*
* @param {Array} poly
* @param {vtkPoints} points
* @param {Array} polyEdges
* @param {Array} originalEdges
* @param {vtkCellArray} polys
* @param {Vector3} normal
* @returns {boolean}
*/
export function vtkCCSTriangulate(
poly,
points,
polyEdges,
originalEdges,
polys,
normal
) {
let n = poly.length;

// If the poly is a line, then skip it
if (n < 3) {
return true;
}

// If the poly is a triangle, then pass it
if (n === 3) {
const trids = [0, 1, 2];
vtkCCSInsertTriangle(polys, poly, trids, polyEdges, originalEdges);
return true;
}

// If the poly has 4 or more points, triangulate it
let triangulationFailure = false;
let ppoint = [];
let point = [];
let npoint = [];
let i = 0;
let j = 0;
let k = 0;

const verts = [];
verts.length = n;

for (i = 0; i < n; i++) {
verts[i] = [i, 0];
}

// compute the triangle quality for each vert
k = n - 2;
points.getPoint(poly[verts[k][0]], point);
i = n - 1;
points.getPoint(poly[verts[i][0]], npoint);

let concave = 0;
let maxq = 0;
let maxi = 0;
for (j = 0; j < n; j++) {
[ppoint, point, npoint] = [point, npoint, ppoint];
points.getPoint(poly[verts[j][0]], npoint);

const q = vtkCCSTriangleQuality(ppoint, point, npoint, normal);
if (q > maxq) {
maxi = i;
maxq = q;
}
concave += q < 0;
verts[i][1] = q;
i = j;
}

let foundEar;
// perform the ear-cut triangulation
for (;;) {
// if no potential ears were found, then fail
if (maxq <= Number.MIN_VALUE) {
triangulationFailure = true;
break;
}

i = maxi;
j = i + 1 !== n ? i + 1 : 0;
k = i !== 0 ? i - 1 : n - 1;

if (verts[i][1] > 0) {
foundEar = true;
points.getPoint(poly[verts[j][0]], npoint);
points.getPoint(poly[verts[k][0]], ppoint);

// only do ear check if there are concave vertices
if (concave) {
// get the normal of the split plane
const v = [];
const u = [];

vtkMath.subtract(npoint, ppoint, v);
vtkMath.cross(v, normal, u);
const d = vtkMath.dot(ppoint, u);

let jj = j + 1 !== n ? j + 1 : 0;
let x = [];
points.getPoint(poly[verts[jj][0]], x);
let side = vtkMath.dot(x, u) < d;
let foundNegative = side;

// check for crossings of the split plane
jj = jj + 1 !== n ? jj + 1 : 0;
let y = [];
const s = [];
const t = [];
for (; foundEar && jj !== k; jj = jj + 1 !== n ? jj + 1 : 0) {
[x, y] = [y, x];
points.getPoint(poly[verts[jj][0]], x);
const sside = vtkMath.dot(x, u) < d;
// XOR
if (side ? !sside : sside) {
side = !side;
foundNegative = true;
foundEar =
vtkLine.intersection(ppoint, npoint, x, y, s, t) ===
vtkLine.IntersectionState.NO_INTERSECTION;
}
}
foundEar &&= foundNegative;
}

if (!foundEar) {
// don't try again until it is split
verts[i][1] = Number.MIN_VALUE;
} else {
// create a triangle from vertex and neighbors
const trids = [verts[i][0], verts[j][0], verts[k][0]];
vtkCCSInsertTriangle(polys, poly, trids, polyEdges, originalEdges);

// remove the vertex i
verts.splice(i, 1);
k -= i === 0;
j -= j !== 0;

// break if this was final triangle
if (--n < 3) {
break;
}

// re-compute quality of previous point
const kk = k !== 0 ? k - 1 : n - 1;
points.getPoint(poly[verts[kk][0]], point);
const kq = vtkCCSTriangleQuality(point, ppoint, npoint, normal);
concave -= verts[k][1] < 0 && kq >= 0;
verts[k][1] = kq;

// re-compute quality of next point
const jj = j + 1 !== n ? j + 1 : 0;
points.getPoint(poly[verts[jj][0]], point);
const jq = vtkCCSTriangleQuality(ppoint, npoint, point, normal);
concave -= verts[j][1] < 0 && jq >= 0;
verts[j][1] = jq;
}
}

// find the highest-quality ear candidate
maxi = 0;
maxq = verts[0][1];
for (i = 1; i < n; i++) {
const q = verts[i][1];
if (q > maxq) {
maxi = i;
maxq = q;
}
}
}

return !triangulationFailure;
}

// ---------------------------------------------------
/**
* Create polygons from line segments.
*
* @param {vtkPolyData} polyData
* @param {Number} firstLine
* @param {Number} endLine
* @param {Boolean} oriented
* @param {Array} newPolys
* @param {Array} incompletePolys
*/
export function vtkCCSMakePolysFromLines(
polyData,
firstLine,
endLine,
oriented,
newPolys,
incompletePolys
) {
let npts = 0;
let pts = [];

// Bitfield for marking lines as used
const usedLines = new Uint8Array(endLine - firstLine); // defaults to 0

// Require cell links to get lines from pointIds
polyData.buildLinks(polyData.getPoints().getNumberOfPoints());

let numNewPolys = 0;
let remainingLines = endLine - firstLine;

while (remainingLines > 0) {
// Create a new poly
const polyId = numNewPolys++;
const poly = [];
newPolys.push(poly);

let lineId = 0;
let completePoly = false;

// start the poly
for (lineId = firstLine; lineId < endLine; lineId++) {
if (!usedLines[lineId - firstLine]) {
pts = polyData.getCellPoints(lineId).cellPointIds;
npts = pts.length;

let n = npts;
if (npts > 2 && pts[0] === pts[npts - 1]) {
n = npts - 1;
completePoly = true;
}
poly.length = n;
for (let i = 0; i < n; i++) {
poly[i] = pts[i];
}
break;
}
}

usedLines[lineId - firstLine] = 1;
remainingLines--;

let noLinesMatch = remainingLines === 0 && !completePoly;

while (!completePoly && !noLinesMatch && remainingLines > 0) {
// This is cleared if a match is found
noLinesMatch = true;

// Number of points in the poly
const npoly = poly.length;

const lineEndPts = [];
const endPts = [poly[npoly - 1], poly[0]];

// For both open ends of the polygon
for (let endIdx = 0; endIdx < 2; endIdx++) {
const matches = [];
const cells = polyData.getPointCells(endPts[endIdx]);

// Go through all lines that contain this endpoint
for (let icell = 0; icell < cells.length; icell++) {
lineId = cells[icell];
if (
lineId >= firstLine &&
lineId < endLine &&
!usedLines[lineId - firstLine]
) {
pts = polyData.getCellPoints(lineId).cellPointIds;
npts = pts.length;
lineEndPts[0] = pts[0];
lineEndPts[1] = pts[npts - 1];

// Check that poly end matches line end
if (
endPts[endIdx] === lineEndPts[endIdx] ||
(!oriented && endPts[endIdx] === lineEndPts[1 - endIdx])
) {
matches.push(lineId);
}
}
}

if (matches.length > 0) {
// Multiple matches mean we need to decide which path to take
if (matches.length > 1) {
// Remove double-backs
let k = matches.length;
do {
lineId = matches[--k];
pts = polyData.getCellPoints(lineId).cellPointIds;
npts = pts.length;
lineEndPts[0] = pts[0];
lineEndPts[1] = pts[npts - 1];
// check if line is reversed
const r = endPts[endIdx] !== lineEndPts[endIdx];

if (
(!r &&
((endIdx === 0 && poly[npoly - 2] === pts[1]) ||
(endIdx === 1 && poly[1] === pts[npts - 2]))) ||
(r &&
((endIdx === 0 && poly[npoly - 2] === pts[npts - 2]) ||
(endIdx === 1 && poly[1] === pts[1])))
) {
matches.splice(k, 1);
}
} while (k > 0 && matches.length > 1);

// If there are multiple matches due to intersections,
// they should be dealt with here.
}

lineId = matches[0];
pts = polyData.getCellPoints(lineId).cellPointIds;
npts = pts.length;
lineEndPts[0] = pts[0];
lineEndPts[1] = pts[npts - 1];

// Do both ends match?
if (endPts[endIdx] === lineEndPts[endIdx]) {
completePoly = endPts[1 - endIdx] === lineEndPts[1 - endIdx];
} else {
completePoly = endPts[1 - endIdx] === lineEndPts[endIdx];
}

if (endIdx === 0) {
for (let i = 1; i < npts - (completePoly ? 1 : 0); i++) {
poly.push(pts[i]);
}
} else {
for (let i = completePoly ? 1 : 0; i < npts - 1; i++) {
poly.unshift(pts[i]);
}
}

if (endPts[endIdx] !== lineEndPts[endIdx]) {
// reverse the ids in the added line
let pit = poly.length;
let ptsIt = completePoly ? 1 : 0;
let ptsEnd = npts - 1;
if (endIdx === 1) {
pit = npts - 1 - (completePoly ? 1 : 0);
ptsIt = 1;
ptsEnd = npts - (completePoly ? 1 : 0);
}
while (ptsIt !== ptsEnd) {
poly[--pit] = poly[ptsIt++];
}
}

usedLines[lineId - firstLine] = 1;
remainingLines--;
noLinesMatch = false;
}
}
}

// Check for incomplete polygons
if (noLinesMatch) {
incompletePolys.push(polyId);
}
}
}

// ---------------------------------------------------
/**
* Join polys that have loose ends, as indicated by incompletePolys.
* Any polys created will have a normal opposite to the supplied normal,
* and any new edges that are created will be on the hull of the point set.
* Shorter edges will be preferred over long edges.
*
* @param {Array[]} polys
* @param {Array} incompletePolys
* @param {vtkPoints} points
* @param {Vector3} normal
*/
export function vtkCCSJoinLooseEnds(polys, incompletePolys, points, normal) {
// Relative tolerance for checking whether an edge is on the hull
const tol = CCS_POLYGON_TOLERANCE;

// A list of polys to remove when everything is done
const removePolys = [];

const p1 = [];
const p2 = [];
let poly1;
let poly2;
let pt1;
let pt2;
let dMin;
let iMin;
let v;
let d;

let n = incompletePolys.length;
while (n !== 0) {
poly1 = polys[incompletePolys[n - 1]];
pt1 = poly1[poly1.length - 1];
points.getPoint(pt1, p1);

dMin = Number.MAX_VALUE;
iMin = 0;

for (let i = 0; i < n; i++) {
poly2 = polys[incompletePolys[i]];
pt2 = poly2[0];
points.getPoint(pt2, p2);

// The next few steps verify that edge [p1, p2] is on the hull
v = [p2[0] - p1[0], p2[1] - p1[1], p2[2] - p1[2]];
d = vtkMath.norm(v);
if (d !== 0) {
v[0] /= d;
v[1] /= d;
v[2] /= d;
}

// Compute the midpoint of the edge
const pm = [
0.5 * (p1[0] + p2[0]),
0.5 * (p1[1] + p2[1]),
0.5 * (p1[2] + p2[2]),
];

// Create a plane equation
const pc = [];
vtkMath.cross(normal, v, pc);
pc[3] = -vtkMath.dot(pc, pm);

// Check that all points are inside the plane. If they aren't, then
// the edge is not on the hull of the pointset.
let badPoint = false;
const m = polys.length;
const p = [];
for (let j = 0; j < m && !badPoint; j++) {
const poly = polys[j];
const npts = poly.length;
for (let k = 0; k < npts; k++) {
const ptId = poly[k];
if (ptId !== pt1 && ptId !== pt2) {
points.getPoint(ptId, p);
const val = p[0] * pc[0] + p[1] * pc[1] + p[2] * pc[2] + pc[3];
const r2 = vtkMath.distance2BetweenPoints(p, pm);

// Check distance from plane against the tolerance
if (val < 0 && val * val > tol * tol * r2) {
badPoint = true;
break;
}
}
}

// If no bad points, then this edge is a candidate
if (!badPoint && d < dMin) {
dMin = d;
iMin = i;
}
}
}

// If a match was found, append the polys
if (dMin < Number.MAX_VALUE) {
// Did the poly match with itself?
if (iMin === n - 1) {
// Mark the poly as closed
incompletePolys.pop();
} else {
const id2 = incompletePolys[iMin];

// Combine the polys
// for (let i = 1; i < polys[id2].length; i++) {
// poly1.push(polys[id2][i]);
// }
poly1.push(...polys[id2]);

// Erase the second poly
removePolys.push(id2);
incompletePolys.splice(iMin, 1);
}
} else {
// If no match, erase this poly from consideration
removePolys.push(incompletePolys[n - 1]);
incompletePolys.pop();
}
n = incompletePolys.length;
}

// Remove polys that couldn't be completed
removePolys.sort((a, b) => a - b);
let i = removePolys.length;
while (i > 0) {
// Remove items in reverse order
polys.splice(removePolys[--i], 1);
}

// Clear the incompletePolys vector, it's indices are no longer valid
incompletePolys.length = 0;
}

// ---------------------------------------------------
/**
* Given three vectors p.p1, p.p2, and p.p3, this routine
* checks to see if progressing from p1 to p2 to p3 is a clockwise
* or counterclockwise progression with respect to the normal.
* The return value is -1 for clockwise, +1 for counterclockwise,
* and 0 if any two of the vectors are coincident.
*
* @param {Vector3} p
* @param {Vector3} p1
* @param {Vector3} p2
* @param {Vector3} p3
* @param {Vector3} normal
* @returns {Number}
*/
export function vtkCCSVectorProgression(p, p1, p2, p3, normal) {
const v1 = [p1[0] - p[0], p1[1] - p[1], p1[2] - p[2]];
const v2 = [p2[0] - p[0], p2[1] - p[1], p2[2] - p[2]];
const v3 = [p3[0] - p[0], p3[1] - p[1], p3[2] - p[2]];
const w1 = [];
const w2 = [];

vtkMath.cross(v2, v1, w1);
vtkMath.cross(v2, v3, w2);
const s1 = vtkMath.dot(w1, normal);
const s2 = vtkMath.dot(w2, normal);

if (s1 !== 0 && s2 !== 0) {
const sb1 = s1 < 0;
const sb2 = s2 < 0;

// if sines have different signs
// XOR
if (sb1 ? !sb2 : sb2) {
// return -1 if s2 is -ve
return 1 - 2 * sb2;
}

const c1 = vtkMath.dot(v2, v1);
const l1 = vtkMath.norm(v1);
const c2 = vtkMath.dot(v2, v3);
const l2 = vtkMath.norm(v3);

// ck is the difference of the cosines, flipped in sign if sines are +ve
const ck = (c2 * l2 - c1 * l1) * (1 - sb1 * 2);

if (ck !== 0) {
// return the sign of ck
return 1 - 2 * (ck < 0);
}
}

return 0;
}

// ---------------------------------------------------
/**
* Check for self-intersection. Split the figure-eights.
* This assumes that all intersections occur at existing
* vertices, i.e. no new vertices will be created. Returns
* the number of splits made.
*
* @param {Array[]} polys
* @param {vtkPoints} points
* @param {Array} polyGroups
* @param {Array} polyEdges
* @param {Vector3} normal
* @param {Boolean} oriented
*/
export function vtkCCSSplitAtPinchPoints(
polys,
points,
polyGroups,
polyEdges,
normal,
oriented
) {
const tryPoints = vtkPoints.newInstance({
dataType: VtkDataTypes.DOUBLE,
empty: true,
});

const locator = vtkIncrementalOctreePointLocator.newInstance();

let splitCount = 0;
let poly;
let n;
let bounds;
let tol;

for (let i = 0; i < polys.length; i++) {
poly = polys[i];
n = poly.length;

bounds = [];
tol =
CCS_POLYGON_TOLERANCE *
Math.sqrt(vtkPolygon.getBounds(poly, points, bounds));

if (tol === 0) {
// eslint-disable-next-line no-continue
continue;
}

tryPoints.initialize();
locator.setTolerance(tol);
locator.initPointInsertion(tryPoints, bounds);

let foundMatch = false;
let idx1 = 0;
let idx2 = 0;
let unique = 0;
const point = [];
const p1 = [];
const p2 = [];
const p3 = [];

for (idx2 = 0; idx2 < n; idx2++) {
points.getPoint(poly[idx2], point);

const { success, pointIdx } = locator.insertUniquePoint(point, 0);
if (!success) {
// Need vertIdx to match poly indices, so force point insertion
locator.insertNextPoint(point);

// Do the points have different pointIds?
idx1 = pointIdx;
unique = poly[idx2] !== poly[idx1];

if (idx2 > idx1 + 2 - unique && n + idx1 > idx2 + 2 - unique) {
if (oriented) {
// Make sure that splitting this poly won't create a hole poly
let prevIdx = n + idx1 - 1;
let midIdx = idx1 + 1;
let nextIdx = idx2 + 1;
if (prevIdx >= n) {
prevIdx -= n;
}
if (midIdx >= n) {
midIdx -= n;
}
if (nextIdx >= n) {
nextIdx -= n;
}

points.getPoint(poly[prevIdx], p1);
points.getPoint(poly[midIdx], p2);
points.getPoint(poly[nextIdx], p3);

if (vtkCCSVectorProgression(point, p1, p2, p3, normal) > 0) {
foundMatch = true;
break;
}
} else {
foundMatch = true;
break;
}
}
}
}

if (foundMatch) {
splitCount++;

// Split off a new poly
const m = idx2 - idx1;
const oldPoly = polys[i];
const oldEdges = polyEdges[i];
const newPoly1 = oldPoly.slice(idx1, idx1 + m + unique);
const newEdges1 = oldEdges.slice(idx1, idx1 + m + unique);
const newPoly2 = new Array(n - m + unique);
const newEdges2 = new Array(n - m + unique);

if (unique) {
newEdges1[m] = -1;
}

// The poly that is split off, which might have more intersections
for (let j = 0; j < idx1 + unique; j++) {
newPoly2[j] = oldPoly[j];
newEdges2[j] = oldEdges[j];
}
if (unique) {
newEdges2[idx1] = -1;
}
for (let k = idx2; k < n; k++) {
newPoly2[k - m + unique] = oldPoly[k];
newEdges2[k - m + unique] = oldEdges[k];
}

polys[i] = newPoly1;
polyEdges[i] = newEdges1;
polys.push(newPoly2);
polyEdges.push(newEdges2);

// Unless polygroup was clear (because poly was reversed),
// make a group with one entry for the new poly
polyGroups.length = polys.length;
if (polyGroups[i].length > 0) {
polyGroups[polys.length - 1].push(polys.length - 1);
}
}
}
return splitCount;
}

// ---------------------------------------------------
/**
* The polygons might have a lot of extra points, i.e. points
* in the middle of the edges. Remove those points, but keep
* the original edges as polylines in the originalEdges array.
* Only original edges with more than two points will be kept.
*
* @param {Array[]} polys
* @param {vtkPoints} points
* @param {Array} polyEdges
* @param {Array} originalEdges
*/
export function vtkCCSFindTrueEdges(polys, points, polyEdges, originalEdges) {
// Tolerance^2 for angle to see if line segments are parallel
const atol2 = CCS_POLYGON_TOLERANCE * CCS_POLYGON_TOLERANCE;

const p0 = [];
const p1 = [];
const p2 = [];
const v1 = [];
const v2 = [];
let l1;
let l2;

for (let polyId = 0; polyId < polys.length; polyId++) {
const oldPoly = polys[polyId];
const n = oldPoly.length;
const newEdges = [];
polyEdges.push(newEdges);

// Only useful if poly has more than three sides
if (n < 4) {
newEdges[0] = -1;
newEdges[1] = -1;
newEdges[2] = -1;
// eslint-disable-next-line no-continue
continue;
}

// While we remove points, m keeps track of how many points are left
let m = n;

// Compute bounds for tolerance
const bounds = [];
const tol2 = vtkPolygon.getBounds(oldPoly, points, bounds) * atol2;

// The new poly
const newPoly = [];
let cornerPointId = 0;
let oldOriginalId = -1;

// Keep the partial edge from before the first corner is found
const partialEdge = [];
let cellCount = 0;

points.getPoint(oldPoly[n - 1], p0);
points.getPoint(oldPoly[0], p1);
vtkMath.subtract(p1, p0, v1);
l1 = vtkMath.dot(v1, v1);

for (let j = 0; j < n; j++) {
let k = j + 1;
if (k >= n) {
k -= n;
}

points.getPoint(oldPoly[k], p2);
vtkMath.subtract(p2, p1, v2);
l2 = vtkMath.dot(v2, v2);

// Dot product is |v1||v2|cos(theta)
const c = vtkMath.dot(v1, v2);
// sin^2(theta) = (1 - cos^2(theta))
// and c*c = l1*l2*cos^2(theta)
const s2 = l1 * l2 - c * c;

// In the small angle approximation, sin(theta) == theta, so
// s2/(l1*l2) is the angle that we want to check, but it's not
// a valid check if l1 or l2 is very close to zero.

const pointId = oldPoly[j];

// Keep the point if:
// 1) removing it would create a 2-point poly OR
// 2) it's more than "tol" distance from the prev point AND
// 3) the angle is greater than atol:
if (
m <= 3 ||
(l1 > tol2 && (c < 0 || l1 < tol2 || l2 < tol2 || s2 > l1 * l2 * atol2))
) {
// Complete the previous edge only if the final point count
// will be greater than two
if (cellCount > 1) {
if (pointId !== oldOriginalId) {
originalEdges.push(pointId);
cellCount++;
}
// Update the number of segments in the edge
const countLocation = originalEdges.length - cellCount - 1;
originalEdges[countLocation] = cellCount;
newEdges.push(countLocation);
} else if (cellCount === 0) {
partialEdge.push(pointId);
} else {
newEdges.push(-1);
}

newPoly.push(pointId);

// Start a new edge with cornerPointId as a "virtual" point
cornerPointId = pointId;
oldOriginalId = pointId;
cellCount = 1;

// Rotate to the next point
p0[0] = p2[0];
p0[1] = p2[1];
p0[2] = p2[2];
p1[0] = p2[0];
p1[1] = p2[1];
p1[2] = p2[2];
v1[0] = v2[0];
v1[1] = v2[1];
v1[2] = v2[2];
l1 = l2;
} else {
if (cellCount > 0 && pointId !== oldOriginalId) {
// First check to see if we have to add cornerPointId
if (cellCount === 1) {
originalEdges.push(1); // new edge
originalEdges.push(cornerPointId);
}
// Then add the new point
originalEdges.push(pointId);
oldOriginalId = pointId;
cellCount++;
} else {
// No corner yet, so save the point
partialEdge.push(pointId);
}

// Reduce the count
m--;

// Join the previous two segments, since the point was removed
p1[0] = p2[0];
p1[1] = p2[1];
p1[2] = p2[2];
vtkMath.subtract(p2, p0, v1);
l1 = vtkMath.dot(v1, v1);
}
}

// Add the partial edge to the end
for (let ii = 0; ii < partialEdge.length; ii++) {
const pointId = partialEdge[ii];
if (pointId !== oldOriginalId) {
if (cellCount === 1) {
originalEdges.push(1); // new edge
originalEdges.push(cornerPointId);
}
originalEdges.push(pointId);
oldOriginalId = pointId;
cellCount++;
}
}

// Finalize
if (cellCount > 1) {
// Update the number of segments in the edge
const countLocation = originalEdges.length - cellCount - 1;
originalEdges[countLocation] = cellCount;
newEdges.push(countLocation);
}
polys[polyId] = newPoly;
}
}

// ---------------------------------------------------
/**
* Reverse a cleaned-up polygon along with the info about
* all of its original vertices.
*
* @param {Array} poly
* @param {Array} edges
* @param {Array} originalEdges
*/
export function vtkCCSReversePoly(poly, edges, originalEdges) {
reverseElements(poly, 1, poly.length - 1);
edges.reverse();
for (let i = 0; i < edges.length; i++) {
if (edges[i] >= 0) {
const firstPtsIdx = edges[i] + 1;
const npts = originalEdges[edges[i]];
reverseElements(originalEdges, firstPtsIdx, firstPtsIdx + npts - 1);
}
}
}

// ---------------------------------------------------
/**
* Check the sense of the polygon against the given normal. Returns
* zero if the normal is zero.
*
* @param {Array} poly
* @param {vtkPoints} points
* @param {Vector3} normal
*/
export function vtkCCSCheckPolygonSense(poly, points, normal) {
// Compute the normal
const pnormal = [0.0, 0.0, 0.0];
const p0 = [];
const p1 = [];
const p2 = [];
const v1 = [];
const v2 = [];
const v = [];

points.getPoint(poly[0], p0);
points.getPoint(poly[1], p1);
vtkMath.subtract(p1, p0, v1);

for (let jj = 2; jj < poly.length; jj++) {
points.getPoint(poly[jj], p2);
vtkMath.subtract(p2, p0, v2);
vtkMath.cross(v1, v2, v);
vtkMath.add(pnormal, v, pnormal);
p1[0] = p2[0];
p1[1] = p2[1];
p1[2] = p2[2];
v1[0] = v2[0];
v1[1] = v2[1];
v1[2] = v2[2];
}

// Check the normal
const d = vtkMath.dot(pnormal, normal);

return { isNormalNotZero: d !== 0, sense: d > 0 };
}

// ---------------------------------------------------
/**
* Check whether innerPoly is inside outerPoly.
* The normal is needed to verify the polygon orientation.
* The values of pp, bounds, and tol2 must be precomputed
* by calling vtkCCSPrepareForPolyInPoly() on outerPoly.
*
* @param {Array} outerPoly
* @param {Array} innerPoly
* @param {vtkPoints} points
* @param {Vector3} normal
* @param {Float64Array} pp
* @param {Bounds} bounds
* @param {Number} tol2
*/
export function vtkCCSPolyInPoly(
outerPoly,
innerPoly,
points,
normal,
pp,
bounds,
tol2
) {
// Find a vertex of poly "j" that isn't on the edge of poly "i".
// This is necessary or the PointInPolygon might return "true"
// based only on roundoff error.
const n = outerPoly.length;
const m = innerPoly.length;
const p = [];
const q1 = [];
const q2 = [];

for (let jj = 0; jj < m; jj++) {
// Semi-randomize the point order
// eslint-disable-next-line no-bitwise
const kk = (jj >> 1) + (jj & 1) * ((m + 1) >> 1);
points.getPoint(innerPoly[kk], p);
const intersectionState = vtkPolygon.pointInPolygon(p, pp, bounds, normal);
if (intersectionState === PolygonWithPointIntersectionState.FAILURE) {
vtkErrorMacro('Error finding point in polygon in vtkCCSPolyInPoly');
}
if (intersectionState !== PolygonWithPointIntersectionState.OUTSIDE) {
let pointOnEdge = 0;
points.getPoint(outerPoly[n - 1], q1);

for (let ii = 0; ii < n; ii++) {
points.getPoint(outerPoly[ii], q2);
// This method returns distance squared
const { distance } = vtkLine.distanceToLine(p, q1, q2);
if (distance < tol2) {
pointOnEdge = 1;
break;
}
q1[0] = q2[0];
q1[1] = q2[1];
q1[2] = q2[2];
}

if (!pointOnEdge) {
// Good result, point is in polygon
return true;
}
}
}

// No matches found
return false;
}

// ---------------------------------------------------
/**
* Precompute values needed for the PolyInPoly check.
* The values that are returned are as follows:
* pp: an array of the polygon vertices
* bounds: the polygon bounds
* tol2: a tolerance value based on the size of the polygon
* (note: pp must be pre-allocated to the 3*outerPoly.length)
*
* @param {Array} outerPoly
* @param {vtkPoints} points
* @param {Float64Array} pp
* @param {Bounds} bounds
*/
export function vtkCCSPrepareForPolyInPoly(outerPoly, points, pp, bounds) {
const n = outerPoly.length;

if (n === 0) {
return 0.0; // to avoid false positive warning about uninitialized value
}

// Pull out the points
const point = [];
let j = 0;
for (let i = 0; i < n; i++) {
points.getPoint(outerPoly[i], point);
pp[j++] = point[0];
pp[j++] = point[1];
pp[j++] = point[2];
}

// Find the bounding box and tolerance for the polygon
return (
vtkPolygon.getBounds(outerPoly, points, bounds) *
(CCS_POLYGON_TOLERANCE * CCS_POLYGON_TOLERANCE)
);
}

// ---------------------------------------------------
/**
* Check for polygons within polygons. Group the polygons
* if they are within each other. Reverse the sense of
* the interior "hole" polygons. A hole within a hole
* will be reversed twice and will become its own group.
*
* @param {Array} newPolys
* @param {vtkPoints} points
* @param {Array} polyGroups
* @param {Array} polyEdges
* @param {Array} originalEdges
* @param {Vector3} normal
* @param {Boolean} oriented
*/
export function vtkCCSMakeHoleyPolys(
newPolys,
points,
polyGroups,
polyEdges,
originalEdges,
normal,
oriented
) {
const numNewPolys = newPolys.length;
if (numNewPolys <= 1) {
return;
}

// Use bit arrays to keep track of inner polys
const polyReversed = [];
const innerPolys = [];

// GroupCount is an array only needed for unoriented polys
let groupCount;
if (!oriented) {
groupCount = new Int32Array(numNewPolys);
}

// Find the maximum poly size
let nmax = 1;
for (let kk = 0; kk < numNewPolys; kk++) {
nmax = Math.max(nmax, newPolys[kk].length);
}

// These are some values needed for poly-in-poly checks
const pp = new Float64Array(3 * nmax);
const bounds = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0];
let tol2;

// Go through all polys
for (let i = 0; i < numNewPolys; i++) {
const n = newPolys[i].length;

if (n < 3) {
// eslint-disable-next-line no-continue
continue;
}

// Check if poly is reversed
const { isNormalNotZero, sense } = vtkCCSCheckPolygonSense(
newPolys[i],
points,
normal
);
if (isNormalNotZero) {
polyReversed[i] = !sense;
}

// Precompute some values needed for poly-in-poly checks
tol2 = vtkCCSPrepareForPolyInPoly(newPolys[i], points, pp, bounds);

// Look for polygons inside of this one
for (let j = 0; j < numNewPolys; j++) {
if (j !== i && newPolys[j].length >= 3) {
// Make sure polygon i is not in polygon j
const pg = polyGroups[j];
if (!pg.includes(i)) {
if (
vtkCCSPolyInPoly(
newPolys[i],
newPolys[j],
points,
normal,
pp.subarray(3 * n),
bounds,
tol2
)
) {
// Add to group
polyGroups[i].push(j);
if (groupCount) {
groupCount[j] += 1;
}
}
}
}
}
}

if (!oriented) {
// build a stack of polys that aren't inside other polys=
const outerPolyStack = [];
for (let ll = 0; ll < numNewPolys; ll++) {
if (groupCount[ll] === 0) {
outerPolyStack.push(ll);
}
}

let j;
while (outerPolyStack.length > 0) {
j = outerPolyStack.length - 1;
outerPolyStack.pop();

if (polyReversed[j]) {
vtkCCSReversePoly(newPolys[j], polyEdges[j], originalEdges);
polyReversed[j] = false;
}

if (polyGroups[j].length > 1) {
// Convert the group into a bit array, to make manipulation easier
innerPolys.length = 0;
for (let k = 1; k < polyGroups[j].length; k++) {
const jj = polyGroups[j][k];
if (groupCount[jj] > 1) {
groupCount[jj] -= 2;
if (groupCount[jj] === 0) {
outerPolyStack.push(jj);
}
} else {
innerPolys[jj] = 1;
polyGroups[jj].length = 0;
if (!polyReversed[jj]) {
vtkCCSReversePoly(newPolys[jj], polyEdges[jj], originalEdges);
polyReversed[jj] = false;
}
}
}

// Use the bit array to recreate the polyGroup
polyGroups[j].length = 0;
polyGroups[j].push(j);
for (let jj = 0; jj < numNewPolys; jj++) {
if (innerPolys[jj]) {
polyGroups[j].push(jj);
}
}
}
}
} else {
// oriented
for (let j = 0; j < numNewPolys; j++) {
// Remove the groups for reversed polys
if (polyReversed[j]) {
polyGroups[j].length = 0;
}
// Polys inside the interior polys have their own groups, so remove
// them from this group
else if (polyGroups[j].length > 1) {
// Convert the group into a bit array, to make manipulation easier
innerPolys.length = 0;
for (let k = 1; k < polyGroups[j].length; k++) {
innerPolys[polyGroups[j][k]] = true;
}

// Look for non-reversed polys inside this one
for (let kk = 1; kk < polyGroups[j].length; kk++) {
// jj is the index of the inner poly
const jj = polyGroups[j][kk];
// If inner poly is not reversed then
if (!polyReversed[jj]) {
// Remove that poly and all polys inside of it from the group
for (let ii = 0; ii < polyGroups[jj].length; ii++) {
innerPolys[polyGroups[jj][ii]] = false;
}
}
}

// Use the bit array to recreate the polyGroup
polyGroups[j].length = 0;
polyGroups[j].push(j);
for (let jj = 0; jj < numNewPolys; jj++) {
if (innerPolys[jj]) {
polyGroups[j].push(jj);
}
}
}
}
}

// delete[] groupCount;
}

// ---------------------------------------------------
/**
* Check line segment with point Ids (i, j) to make sure that it
* doesn't cut through the edges of any polys in the group.
* Return value of zero means check failed and the cut is not
* usable.
*
* @param {Array[]} polys
* @param {vtkPoints} points
* @param {Vector3} normal
* @param {Array} polyGroup
* @param {Number} outerPolyId
* @param {Number} innerPolyId
* @param {Number} outerIdx
* @param {Number} innerIdx
*/
export function vtkCCSCheckCut(
polys,
points,
normal,
polyGroup,
outerPolyId,
innerPolyId,
outerIdx,
innerIdx
) {
const ptId1 = polys[outerPolyId][outerIdx];
const ptId2 = polys[innerPolyId][innerIdx];

const tol = CCS_POLYGON_TOLERANCE;

const p1 = [];
const p2 = [];
points.getPoint(ptId1, p1);
points.getPoint(ptId2, p2);

const w = [];
vtkMath.subtract(p2, p1, w);
const l = vtkMath.normalize(w);

// Cuts between coincident points are good
if (l === 0) {
return true;
}

// Define a tolerance with units of distance squared
const tol2 = l * l * tol * tol;

// Check the sense of the cut: it must be pointing "in" for both polys.
let polyId = outerPolyId;
let polyIdx = outerIdx;

let r = p1;
const r1 = [];
let r2 = p2;
const r3 = [];

for (let ii = 0; ii < 2; ii++) {
const poly = polys[polyId];
const n = poly.length;
let prevIdx = n - polyIdx - 1;
let nextIdx = polyIdx + 1;
if (prevIdx >= n) {
prevIdx -= n;
}
if (nextIdx >= n) {
nextIdx -= n;
}

points.getPoint(poly[prevIdx], r1);
points.getPoint(poly[nextIdx], r3);

if (vtkCCSVectorProgression(r, r1, r2, r3, normal) > 0) {
return false;
}

polyId = innerPolyId;
polyIdx = innerIdx;
r = p2;
r2 = p1;
}

// Check for intersections of the cut with polygon edges.
// First, create a cut plane that divides space at the cut line.
const pc = [];
vtkMath.cross(normal, w, pc);
pc[3] = -vtkMath.dot(pc, p1);

for (let i = 0; i < polyGroup.length; i++) {
const poly = polys[polyGroup[i]];
const n = poly.length;

const q1 = [];
const q2 = [];
let qtId1 = poly[n - 1];
points.getPoint(qtId1, q1);
let v1 = pc[0] * q1[0] + pc[1] * q1[1] + pc[2] * q1[2] + pc[3];
let c1 = v1 > 0;

for (let j = 0; j < n; j++) {
const qtId2 = poly[j];
points.getPoint(qtId2, q2);
const v2 = pc[0] * q2[0] + pc[1] * q2[1] + pc[2] * q2[2] + pc[3];
const c2 = v2 > 0;

// If lines share an endpoint, they can't intersect,
// so don't bother with the check.
if (
ptId1 !== qtId1 &&
ptId1 !== qtId2 &&
ptId2 !== qtId1 &&
ptId2 !== qtId2
) {
// Check for intersection
if ((c1 ? !c2 : c2) || v1 * v1 < tol2 || v2 * v2 < tol2) {
vtkMath.subtract(q2, q1, w);
if (vtkMath.dot(w, w) > 0) {
const qc = [];
vtkMath.cross(w, normal, qc);
qc[3] = -vtkMath.dot(qc, q1);

const u1 = qc[0] * p1[0] + qc[1] * p1[1] + qc[2] * p1[2] + qc[3];
const u2 = qc[0] * p2[0] + qc[1] * p2[1] + qc[2] * p2[2] + qc[3];
const d1 = u1 > 0;
const d2 = u2 > 0;

if (d1 ? !d2 : d2) {
// One final check to make sure endpoints aren't coincident
let p = p1;
let q = q1;
if (v2 * v2 < v1 * v1) {
p = p2;
}
if (u2 * u2 < u1 * u1) {
q = q2;
}
if (vtkMath.distance2BetweenPoints(p, q) > tol2) {
return false;
}
}
}
}
}

qtId1 = qtId2;
q1[0] = q2[0];
q1[1] = q2[1];
q1[2] = q2[2];
v1 = v2;
c1 = c2;
}
}

return true;
}

// ---------------------------------------------------
/**
* Check the quality of a cut between an outer and inner polygon.
* An ideal cut is one that forms a 90 degree angle with each
* line segment that it joins to. Smaller values indicate a
* higher quality cut.
*
* @param {Array} outerPoly
* @param {Array} innerPoly
* @param {Number} i
* @param {Number} j
* @param {vtkPoints} points
*/
export function vtkCCSCutQuality(outerPoly, innerPoly, i, j, points) {
const n = outerPoly.length;
const m = innerPoly.length;

const a = i > 0 ? i - 1 : n - 1;
const b = i < n - 1 ? i + 1 : 0;

const c = j > 0 ? j - 1 : m - 1;
const d = j < m - 1 ? j + 1 : 0;

const p0 = [];
const p1 = [];
const p2 = [];
points.getPoint(outerPoly[i], p1);
points.getPoint(innerPoly[j], p2);

const v1 = [];
const v2 = [];
vtkMath.subtract(p2, p1, v1);

const l1 = vtkMath.dot(v1, v1);
let l2;
let qmax = 0;
let q;

points.getPoint(outerPoly[a], p0);
vtkMath.subtract(p0, p1, v2);
l2 = vtkMath.dot(v2, v2);
if (l2 > 0) {
q = vtkMath.dot(v1, v2);
q *= q / l2;
if (q > qmax) {
qmax = q;
}
}

points.getPoint(outerPoly[b], p0);
vtkMath.subtract(p0, p1, v2);
l2 = vtkMath.dot(v2, v2);
if (l2 > 0) {
q = vtkMath.dot(v1, v2);
q *= q / l2;
if (q > qmax) {
qmax = q;
}
}

points.getPoint(innerPoly[c], p0);
vtkMath.subtract(p2, p0, v2);
l2 = vtkMath.dot(v2, v2);
if (l2 > 0) {
q = vtkMath.dot(v1, v2);
q *= q / l2;
if (q > qmax) {
qmax = q;
}
}

points.getPoint(innerPoly[d], p0);
vtkMath.subtract(p2, p0, v2);
l2 = vtkMath.dot(v2, v2);
if (l2 > 0) {
q = vtkMath.dot(v1, v2);
q *= q / l2;
if (q > qmax) {
qmax = q;
}
}

if (l1 > 0) {
return qmax / l1; // also l1 + qmax, incorporates distance;
}

return Number.MAX_VALUE;
}

// ---------------------------------------------------
/**
* Find the two sharpest verts on an inner (i.e. inside-out) poly.
*
* @param {Array} poly
* @param {vtkPoints} points
* @param {Vector3} normal
* @param {[Number, Number]} verts
*/
export function vtkCCSFindSharpestVerts(poly, points, normal, verts) {
const p1 = [];
const p2 = [];
const v1 = [];
const v2 = [];
const v = [];
let l1;
let l2;

const minVal = [0, 0];

verts[0] = 0;
verts[1] = 0;

const n = poly.length;
points.getPoint(poly[n - 1], p2);
points.getPoint(poly[0], p1);

vtkMath.subtract(p1, p2, v1);
l1 = Math.sqrt(vtkMath.dot(v1, v1));

for (let j = 0; j < n; j++) {
let k = j + 1;
if (k === n) {
k = 0;
}

points.getPoint(poly[k], p2);
vtkMath.subtract(p2, p1, v2);
l2 = Math.sqrt(vtkMath.dot(v2, v2));

vtkMath.cross(v1, v2, v);
const b = vtkMath.dot(v, normal);

if (b < 0 && l1 * l2 > 0) {
// Dot product is |v1||v2|cos(theta), range [-1, +1]
const val = vtkMath.dot(v1, v2) / (l1 * l2);
if (val < minVal[0]) {
minVal[1] = minVal[0];
minVal[0] = val;
verts[1] = verts[0];
verts[0] = j;
}
}

// Rotate to the next point
p1[0] = p2[0];
p1[1] = p2[1];
p1[2] = p2[2];
v1[0] = v2[0];
v1[1] = v2[1];
v1[2] = v2[2];
l1 = l2;
}
}

// ---------------------------------------------------
/**
* Find two valid cuts between outerPoly and innerPoly.
* Used by vtkCCSCutHoleyPolys.
*
* @param {Array} polys
* @param {Array} polyGroup
* @param {Number} outerPolyId
* @param {Number} innerPolyId
* @param {vtkPoints} points
* @param {Vector3} normal
* @param {Array[]} cuts
* @param {Boolean} exhaustive
*/
export function vtkCCSFindCuts(
polys,
polyGroup,
outerPolyId,
innerPolyId,
points,
normal,
cuts,
exhaustive
) {
const outerPoly = polys[outerPolyId];
const innerPoly = polys[innerPolyId];
const innerSize = innerPoly.length;
// Find the two sharpest points on the inner poly
const verts = [];
vtkCCSFindSharpestVerts(innerPoly, points, normal, verts);

// A list of cut locations according to quality
const cutlist = [];
cutlist.length = outerPoly.length;

// Search for potential cuts (need to find two cuts)
let cutId = 0;
cuts[0][0] = 0;
cuts[0][1] = 0;
cuts[1][0] = 0;
cuts[1][1] = 0;

let foundCut = false;
for (cutId = 0; cutId < 2; cutId++) {
const count = exhaustive ? innerSize : 3;

for (let i = 0; i < count && !foundCut; i++) {
// Semi-randomize the search order
// TODO: Does this do the same as in C++?
// eslint-disable-next-line no-bitwise
let j = (i >> 1) + (i & 1) * ((innerSize + 1) >> 1);
// Start at the best first point
j = (j + verts[cutId]) % innerSize;

for (let kk = 0; kk < outerPoly.length; kk++) {
const q = vtkCCSCutQuality(outerPoly, innerPoly, kk, j, points);
cutlist[kk] = [q, kk];
}

cutlist.sort((a, b) => a[0] - b[0]);

for (let lid = 0; lid < cutlist.length; lid++) {
const k = cutlist[lid][1];

// If this is the second cut, do extra checks
if (cutId > 0) {
// Make sure cuts don't share an endpoint
if (j === cuts[0][1] || k === cuts[0][0]) {
// eslint-disable-next-line no-continue
continue;
}

// Make sure cuts don't intersect
const p1 = [];
const p2 = [];
points.getPoint(outerPoly[cuts[0][0]], p1);
points.getPoint(innerPoly[cuts[0][1]], p2);

const q1 = [];
const q2 = [];
points.getPoint(outerPoly[k], q1);
points.getPoint(innerPoly[j], q2);

let u;
let v;
if (
vtkLine.intersection(p1, p2, q1, q2, u, v) ===
vtkLine.IntersectionState.YES_INTERSECTION
) {
// eslint-disable-next-line no-continue
continue;
}
}

// This check is done for both cuts
if (
vtkCCSCheckCut(
polys,
points,
normal,
polyGroup,
outerPolyId,
innerPolyId,
k,
j
)
) {
cuts[cutId][0] = k;
cuts[cutId][1] = j;
foundCut = true;
break;
}
}
}

if (!foundCut) {
return false;
}
}

return true;
}

// ---------------------------------------------------
/**
* Helper for vtkCCSCutHoleyPolys. Change a polygon and a hole
* into two separate polygons by making two cuts between them.
*
* @param {Array[]} polys
* @param {Array} polyEdges
* @param {Number} outerPolyId
* @param {Number} innerPolyId
* @param {vtkPoints} points
* @param {Array[]} cuts
*/
export function vtkCCSMakeCuts(
polys,
polyEdges,
outerPolyId,
innerPolyId,
points,
cuts
) {
const q = [];
const r = [];
for (let bb = 0; bb < 2; bb++) {
const ptId1 = polys[outerPolyId][cuts[bb][0]];
const ptId2 = polys[innerPolyId][cuts[bb][1]];
points.getPoint(ptId1, q);
points.getPoint(ptId2, r);
}

const outerPoly = polys[outerPolyId];
const innerPoly = polys[innerPolyId];

const outerEdges = polyEdges[outerPolyId];
const innerEdges = polyEdges[innerPolyId];

// Generate new polys from the cuts
const n = outerPoly.length;
const m = innerPoly.length;
let idx;

// Generate poly1
const n1 = n * (cuts[1][0] < cuts[0][0]) + cuts[1][0] - cuts[0][0] + 1;
const n2 = n1 + m * (cuts[0][1] < cuts[1][1]) + cuts[0][1] - cuts[1][1] + 1;

const poly1 = [];
poly1.length = n2;
const edges1 = new Array(n2);

idx = cuts[0][0];
for (let i1 = 0; i1 < n1; i1++) {
const k = idx++;
poly1[i1] = outerPoly[k];
edges1[i1] = outerEdges[k];
idx *= idx !== n;
}
edges1[n1 - 1] = -1;

idx = cuts[1][1];
for (let i2 = n1; i2 < n2; i2++) {
const k = idx++;
poly1[i2] = innerPoly[k];
edges1[i2] = innerEdges[k];
idx *= idx !== m;
}
edges1[n2 - 1] = -1;

// Generate poly2
const m1 = n * (cuts[0][0] < cuts[1][0]) + cuts[0][0] - cuts[1][0] + 1;
const m2 = m1 + m * (cuts[1][1] < cuts[0][1]) + cuts[1][1] - cuts[0][1] + 1;

const poly2 = [];
poly2.length = m2;
const edges2 = new Array(m2);

idx = cuts[1][0];
for (let j1 = 0; j1 < m1; j1++) {
const k = idx++;
poly2[j1] = outerPoly[k];
edges2[j1] = outerEdges[k];
idx *= idx !== n;
}
edges2[m1 - 1] = -1;

idx = cuts[0][1];
for (let j2 = m1; j2 < m2; j2++) {
const k = idx++;
poly2[j2] = innerPoly[k];
edges2[j2] = innerEdges[k];
idx *= idx !== m;
}
edges2[m2 - 1] = -1;

// Replace outerPoly and innerPoly with these new polys
polys[outerPolyId] = poly1;
polys[innerPolyId] = poly2;
polyEdges[outerPolyId] = edges1;
polyEdges[innerPolyId] = edges2;
}

// ---------------------------------------------------
/**
* After the holes have been identified, make cuts between the
* outer poly and each hole. Make two cuts per hole. The only
* strict requirement is that the cut must not intersect any
* edges, but it's best to make sure that no really sharp angles
* are created.
*
* @param {Array[]} polys
* @param {vtkPoints} points
* @param {Array[]} polyGroups
* @param {Array} polyEdges
* @param {Vector3} normal
* @returns {boolean}
*/
export function vtkCCSCutHoleyPolys(
polys,
points,
polyGroups,
polyEdges,
normal
) {
let cutFailure = 0;

// Go through all groups and cut out the first inner poly that is
// found. Every time an inner poly is cut out, the groupId counter
// is reset because cutting a poly creates a new group.
let groupId = 0;
while (groupId < polyGroups.length) {
const polyGroup = polyGroups[groupId];

// Only need to make a cut if the group size is greater than 1
if (polyGroup.length > 1) {
// The first member of the group is the outer poly
const outerPolyId = polyGroup[0];

// The second member of the group is the first inner poly
let innerPolyId = polyGroup[1];

// Sort the group by size, do largest holes first
let innerBySize = new Array(polyGroup.length);

for (let i = 1; i < polyGroup.length; i++) {
innerBySize[i] = [polys[polyGroup[i]].length, i];
}

innerBySize = [
innerBySize[0],
...innerBySize.splice(1).sort((a, b) => a[0] - b[0]),
];
reverseElements(innerBySize, 1, innerBySize.length - 1);

// Need to check all inner polys in sequence, until one succeeds.
// Do a quick search first, then do an exhaustive search.
let madeCut = 0;
let inner = 0;
for (let exhaustive = 0; exhaustive < 2 && !madeCut; exhaustive++) {
for (let j = 1; j < polyGroup.length; j++) {
inner = innerBySize[j][1];
innerPolyId = polyGroup[inner];

const cuts = [];
if (
vtkCCSFindCuts(
polys,
polyGroup,
outerPolyId,
innerPolyId,
points,
normal,
cuts,
exhaustive
)
) {
vtkCCSMakeCuts(
polys,
polyEdges,
outerPolyId,
innerPolyId,
points,
cuts
);
madeCut = 1;
break;
}
}
}

if (madeCut) {
// Move successfully cut innerPolyId into its own group
polyGroup.splice(inner, 1);
// Only add if innerPolyId hasn't been set already.
// Having the same poly occur as both polyGroup and
// innerPoly would cause an infinite loop.
if (polyGroups[innerPolyId].length === 0) {
polyGroups[innerPolyId].push(innerPolyId);
}
} else {
// Remove all failed inner polys from the group
for (let k = 1; k < polyGroup.length; k++) {
innerPolyId = polyGroup[k];
// Only add if innerPolyId hasn't been set already.
// Having the same poly occur as both polyGroup and
// innerPoly would cause an infinite loop.
if (polyGroups[innerPolyId].length === 0) {
polyGroups[innerPolyId].push(innerPolyId);
}
}
polyGroup.length = 1;
cutFailure = 1;
}

// If there are other interior polys in the group, find out whether
// they are in poly1 or poly2
if (polyGroup.length > 1) {
const poly1 = polys[outerPolyId];
const pp = new Float64Array(3 * poly1.length);
const bounds = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0];
const tol2 = vtkCCSPrepareForPolyInPoly(poly1, points, pp, bounds);

let nextGroupId = groupId;
let ii = 1;
while (ii < polyGroup.length) {
if (
vtkCCSPolyInPoly(
poly1,
polys[polyGroup[ii]],
points,
normal,
pp,
bounds,
tol2
)
) {
// Keep this poly in polyGroup
ii++;
} else {
// Move this poly to poly2 group
polyGroups[innerPolyId].push(polyGroup[ii]);
polyGroup.splice(ii, 1);

// Reduce the groupId to ensure that this new group will get cut
if (innerPolyId < nextGroupId) {
nextGroupId = innerPolyId;
}
}
}

// Set the groupId for the next iteration
groupId = nextGroupId;
// eslint-disable-next-line no-continue
continue;
}
}
// Increment to the next group
groupId++;
}
return !cutFailure;
}
index.d.ts
import { vtkAlgorithm, vtkObject } from "../../../interfaces";
import { Nullable, TypedArray, Vector3 } from "../../../types";
import vtkCellArray from "../../../Common/Core/CellArray";
import vtkPolyData from "../../../Common/DataModel/PolyData";
import vtkPoints from "../../../Common/Core/Points";

/**
*
*/
export interface IContourTriangulatorInitialValues {
triangulatePolys?: boolean;
}

type vtkContourTriangulatorBase = vtkObject & vtkAlgorithm;

export interface vtkContourTriangulator extends vtkContourTriangulatorBase {
/**
*
* @param {any} inData
* @param {any} outData
*/
requestData(inData: any, outData: any): void;

/**
* Sets the behavior of the filter regarding polys.
* @param {boolean} triangulate whether the filter should triangulate polys
* or leave them untouched. True by default
* @return {boolean} true if it changes
*/
setTriangulatePolys(triangulate: boolean): boolean;

/**
* Returns the behavior of the filter regarding polys.
* @return {boolean} True if the filter triangulates polys, false otherwise.
*/
getTriangulatePolys(): boolean;
}

// ----------------------------------------------------------------------------
// Static API
// ----------------------------------------------------------------------------

/**
* This is a complex subroutine that takes a collection of lines that
* were formed by cutting a polydata with a plane, and generates
* a face that has those lines as its edges. The lines must form one
* or more closed contours, but they need not be sorted.
*
* Only "numLine" lines starting from "firstLine" are used to create new
* polygons, and the new polygons are appended to "polys". The normal of
* the cut plane must be provided so that polys will be correctly oriented.
*
* Given some closed contour lines, create a triangle mesh that fills
* those lines. The input lines do not have to be in tail-to-tip order.
* Only numLines starting from firstLine will be used. Note that holes
* can be indicated by contour loops whose normals are in the opposite
* direction to the provided normal.
*
* @param {vtkPolyData} polyData
* @param {Number} firstLine
* @param {Number} numLines
* @param {vtkCellArray} polys
* @param {Nullable<Vector3>} normal If null, the function will compute
* the normal of the polys.
* @param {Boolean} [triangulatePolys] (default: true) If set to true
* the resulting polygons will be triangulated, otherwise the polygons
* themselves will be added to the output.
* @param {Boolean} [diagnoseOnTriangulationError] (default: false)
* If this option is set to true and there was a triangulation error
* this will add the polys as outlines to the output.
* @returns {Boolean} Returns true if triangulation was successful,
* false otherwise.
*/
export function triangulateContours(
polyData: vtkPolyData,
firstLine: number,
numLines: number,
polys: vtkCellArray,
normal: Nullable<Vector3>,
triangulatePolys?: boolean,
diagnoseOnTriangulationError?: boolean
): boolean;

/**
* A robust method for triangulating a polygon. It cleans up the polygon
* and then applies the ear-cut triangulation. A zero return value
* indicates that triangulation failed.
*
* @param {Array<Number>|TypedArray} polygon Array of point indices defining the polygon
* @param {vtkPoints} points The point coordinates of the polygon
* @param {vtkCellArray} triangles The cell array that is going to be
* filled with the triangulation
* @returns {Boolean} Returns true if triangulation was successful,
* false otherwise.
*/
export function triangulatePolygon(
polygon: Array<number> | TypedArray,
points: vtkPoints,
triangles: vtkCellArray
): boolean;

/**
* Method use to decorate a given object (publicAPI+model) with vtkContourTriangulator characteristics.
*
* @param publicAPI object on which methods will be bounds (public)
* @param model object on which data structure will be bounds (protected)
* @param {IContourTriangulatorInitialValues} [initialValues] (default: {})
*/
export function extend(
publicAPI: object,
model: object,
initialValues?: IContourTriangulatorInitialValues
): void;

// ----------------------------------------------------------------------------

/**
* Method use to create a new instance of vtkContourTriangulator
* @param {IContourTriangulatorInitialValues} [initialValues] for pre-setting some of its content
*/
export function newInstance(
initialValues?: IContourTriangulatorInitialValues
): vtkContourTriangulator;

/**
* vtkContourTriangulator
*/
export declare const vtkContourTriangulator: {
newInstance: typeof newInstance;
extend: typeof extend;
// static
triangulateContours: typeof triangulateContours;
triangulatePolygon: typeof triangulatePolygon;
};

export default vtkContourTriangulator;
index.js
import macro from 'vtk.js/Sources/macros';
import vtkCellArray from 'vtk.js/Sources/Common/Core/CellArray';
import vtkPolygon from 'vtk.js/Sources/Common/DataModel/Polygon';
import vtkPolyData from 'vtk.js/Sources/Common/DataModel/PolyData';
import { VtkDataTypes } from 'vtk.js/Sources/Common/Core/DataArray/Constants';

import {
vtkCCSCutHoleyPolys,
vtkCCSFindTrueEdges,
vtkCCSJoinLooseEnds,
vtkCCSMakeHoleyPolys,
vtkCCSMakePolysFromLines,
vtkCCSSplitAtPinchPoints,
vtkCCSTriangulate,
} from './helper';

const { vtkErrorMacro } = macro;

const TRIANGULATION_ERROR_DISPLAY = false;
const DIAGNOSE_ON_TRIANGULATION_ERROR = false;

//------------------------------------------------------------------------------
function triangulateContours(
polyData,
firstLine,
numLines,
polys,
normal,
triangulatePolys = true
) {
let triangulationFailure = false;

// If no cut lines were generated, there's nothing to do
if (numLines <= 0) {
return false;
}

const points = polyData.getPoints();

// Join all the new lines into connected groups, i.e. polygons.
// If we are lucky these will be simple, convex polygons. But
// we can't count on that.

const newPolys = [];
const incompletePolys = [];

const oriented = normal?.length < 3;
vtkCCSMakePolysFromLines(
polyData,
firstLine,
firstLine + numLines,
oriented,
newPolys,
incompletePolys
);

// if no normal specified, then compute one from largest contour
let computedNormal = normal;
if (!oriented) {
computedNormal = [0, 0, 1];
let maxnorm = 0;
const n = [];
for (let i = 0; i < newPolys.length; i++) {
const norm = vtkPolygon.getNormal(newPolys[i], points, n);
if (norm > maxnorm) {
maxnorm = norm;
computedNormal[0] = n[0];
computedNormal[1] = n[1];
computedNormal[2] = n[2];
}
}
}

// Join any loose ends. If the input was a closed surface then there
// will not be any loose ends, so this is provided as a service to users
// who want to clip a non-closed surface.
vtkCCSJoinLooseEnds(newPolys, incompletePolys, points, computedNormal);

// Some points might be in the middle of straight line segments.
// These points can be removed without changing the shape of the
// polys, and removing them makes triangulation more stable.
// Unfortunately removing these points also means that the polys
// will no longer form a watertight cap over the cut.

const polyEdges = [];
const originalEdges = [];
vtkCCSFindTrueEdges(newPolys, points, polyEdges, originalEdges);

// Next we have to check for polygons with holes, i.e. polygons that
// have other polygons inside. Each polygon is "grouped" with the
// polygons that make up its holes.

// Initialize each group to hold just one polygon.
const numNewPolys = newPolys.length;
const polyGroups = new Array(numNewPolys);
for (let i = 0; i < numNewPolys; i++) {
polyGroups[i] = [i];
}

// Find out which polys are holes in larger polys. Create a group
// for each poly where the first member of the group is the larger
// poly, and all other members are the holes. The number of polyGroups
// will be the same as the number of polys, and any polys that are
// holes will have a matching empty group.

vtkCCSMakeHoleyPolys(
newPolys,
points,
polyGroups,
polyEdges,
originalEdges,
computedNormal,
oriented
);

// Make cuts to create simple polygons out of the holey polys.
// After this is done, each polyGroup will have exactly 1 polygon,
// and no polys will be holes. This is currently the most computationally
// expensive part of the process.

if (
!vtkCCSCutHoleyPolys(
newPolys,
points,
polyGroups,
polyEdges,
computedNormal
)
) {
triangulationFailure = true;
}

// Some polys might be self-intersecting. Split the polys at each intersection point.
vtkCCSSplitAtPinchPoints(
newPolys,
points,
polyGroups,
polyEdges,
computedNormal,
oriented
);

// ------ Triangulation code ------

// Go through all polys and triangulate them
for (let polyId = 0; polyId < polyGroups.length; polyId++) {
// If group is empty, then poly was a hole without a containing poly
if (polyGroups[polyId].length === 0) {
// eslint-disable-next-line no-continue
continue;
}

if (!triangulatePolys) {
polys.insertNextCell(originalEdges.slice(1, originalEdges.length));
} else if (
!vtkCCSTriangulate(
newPolys[polyId],
points,
polyEdges[polyId],
originalEdges,
polys,
computedNormal
)
) {
triangulationFailure = false;
// Diagnostic code: show the polys as outlines
if (DIAGNOSE_ON_TRIANGULATION_ERROR) {
const lines = polyData.getLines();
const poly = newPolys[polyId];
lines.insertNextCell([poly.length + 1, ...poly, poly[0]]);
}
}
}

return !triangulationFailure;
}

// ---------------------------------------------------
function triangulatePolygon(polygon, points, triangles) {
const poly = [...polygon];
const polys = [poly];

const originalEdges = [];
const polyEdges = [];
vtkCCSFindTrueEdges(polys, points, polyEdges, originalEdges);
const edges = polyEdges[0];

let success = true;
const normal = [];
const norm = vtkPolygon.getNormal(poly, points, normal);
if (norm !== 0) {
success = vtkCCSTriangulate(
poly,
points,
edges,
originalEdges,
triangles,
normal
);
}
return success;
}

export const STATIC = { triangulateContours, triangulatePolygon };

function vtkContourTriangulator(publicAPI, model) {
// Set our className
model.classHierarchy.push('vtkContourTriangulator');

publicAPI.requestData = (inData, outData) => {
// implement requestData
const input = inData[0];
// FIXME: do not instantiate a new polydata each time the filter is executed.
const output = vtkPolyData.newInstance();
outData[0] = output;

if (!input) {
vtkErrorMacro('Invalid or missing input');
return false;
}

let triangulationError = false;

const lines = input.getLines();
if (lines == null || lines.getNumberOfCells() === 0) {
return true;
}

input.buildCells();

const polysArray = vtkCellArray.newInstance({
dataType: VtkDataTypes.DOUBLE,
empty: true,
});
output.setPolys(polysArray);
output.setPoints(input.getPoints());
output.getPointData().passData(input.getPointData());

triangulationError = !triangulateContours(
input,
input.getNumberOfVerts(),
lines.getNumberOfCells(),
polysArray,
null,
model.triangulatePolys
);

if (triangulationError && TRIANGULATION_ERROR_DISPLAY) {
vtkErrorMacro('Triangulation failed, output might have holes.');
}

return true;
};
}

// ----------------------------------------------------------------------------
// Object factory
// ----------------------------------------------------------------------------

const DEFAULT_VALUES = {
triangulatePolys: true,
};

// ----------------------------------------------------------------------------

export function extend(publicAPI, model, initialValues = {}) {
Object.assign(model, DEFAULT_VALUES, initialValues);

// Make this a VTK object
macro.obj(publicAPI, model);

// Also make it an algorithm with one input and one output
macro.algo(publicAPI, model, 1, 1);

macro.setGet(publicAPI, model, ['triangulatePolys']);

// Object specific methods
vtkContourTriangulator(publicAPI, model);
}

// ----------------------------------------------------------------------------

export const newInstance = macro.newInstance(extend, 'vtkContourTriangulator');

// ----------------------------------------------------------------------------

export default { newInstance, extend, ...STATIC };