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..9fcd07cbb88 100644 --- a/Sources/Common/DataModel/Polygon/index.d.ts +++ b/Sources/Common/DataModel/Polygon/index.d.ts @@ -1,66 +1,240 @@ 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; } export interface vtkPolygon extends vtkObject { + /** + * Set the polygon's points + * Points must be ordered in counterclockwise order + * @param {Vector3[]|Array} points The polygon's points. + * @param {Array} pointIds pointIds + */ + setPoints(points: Vector3[]|Array, pointIds?: Array): 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; + + /** + * 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; + /** * 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|TypedArray} vertices Vertices of the polygon + * @param {Array|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|TypedArray, - bounds: Bounds, - normal: Vector3 -): PolygonIntersectionState; + point: Vector3, + vertices: Array | 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; /** * Method used to decorate a given object (publicAPI+model) with vtkPolygon characteristics. @@ -69,7 +243,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 +257,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..3a211325aff 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,684 @@ 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|TypedArray} 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|TypedArray} 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 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 +847,18 @@ export function getNormal(poly, points, normal) { const STATIC = { PolygonWithPointIntersectionState, + PolygonWithCellIntersectionState, pointInPolygon, getBounds, getNormal, + isConvex, + computeCentroid, + computeArea, + evaluatePosition, + intersectWithLine, + intersectConvex2DCells, + interpolateFunctions, + parameterizePolygon, }; // ---------------------------------------------------------------------------- @@ -234,36 +869,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 +909,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 +967,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 +984,135 @@ 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.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 +1120,7 @@ function vtkPolygon(publicAPI, model) { // ---------------------------------------------------------------------------- const DEFAULT_VALUES = { - firstPoint: null, - pointCount: 0, - tris: [], + points: [], }; // ---------------------------------------------------------------------------- @@ -447,9 +1128,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); + if (initialValues.points) { + model.normal = []; + 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..cc66034297c --- /dev/null +++ b/Sources/Common/DataModel/Polygon/test/testPolygon.js @@ -0,0 +1,514 @@ +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 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++) {