diff --git a/Sources/Common/DataModel/AbstractPointLocator/index.d.ts b/Sources/Common/DataModel/AbstractPointLocator/index.d.ts new file mode 100644 index 00000000000..8b879633bf8 --- /dev/null +++ b/Sources/Common/DataModel/AbstractPointLocator/index.d.ts @@ -0,0 +1,54 @@ +import { vtkObject } from "../../../interfaces"; +import { Bounds } from "../../../types"; +import { ILocatorInitialValues } from "../Locator"; + +/** + * + */ +export interface IAbstractPointLocatorInitialValues + extends ILocatorInitialValues { + bounds?: Bounds; + numberOfBuckets: number; +} + +export interface vtkAbstractPointLocator extends vtkObject { + /** + * Set the bounds of this object. + * @param {Bounds} input + */ + setBounds(input: Bounds): void; + + /** + * Get the bounds of this object. + * @returns {Bounds} + */ + getBounds(): Bounds; +} + +// ---------------------------------------------------------------------------- +// Static API +// ---------------------------------------------------------------------------- + +/** + * Method use to decorate a given object (publicAPI+model) with vtkAbstractPointLocator characteristics. + * + * @param publicAPI object on which methods will be bounds (public) + * @param model object on which data structure will be bounds (protected) + * @param {IAbstractPointLocatorInitialValues} [initialValues] (default: {}) + */ +export function extend( + publicAPI: object, + model: object, + initialValues?: IAbstractPointLocatorInitialValues +): void; + +// ---------------------------------------------------------------------------- + +/** + * vtkAbstractPointLocator + */ +export declare const vtkAbstractPointLocator: { + extend: typeof extend; +}; + +export default vtkAbstractPointLocator; diff --git a/Sources/Common/DataModel/AbstractPointLocator/index.js b/Sources/Common/DataModel/AbstractPointLocator/index.js new file mode 100644 index 00000000000..342e5e4b94d --- /dev/null +++ b/Sources/Common/DataModel/AbstractPointLocator/index.js @@ -0,0 +1,39 @@ +import macro from 'vtk.js/Sources/macros'; +import vtkLocator from 'vtk.js/Sources/Common/DataModel/Locator'; + +function vtkAbstractPointLocator(publicAPI, model) { + // Set our className + model.classHierarchy.push('vtkAbstractPointLocator'); +} + +// ---------------------------------------------------------------------------- +// Object factory +// ---------------------------------------------------------------------------- + +function defaultValues(initialValues) { + return { + bounds: [], + numberOfBuckets: 0, + ...initialValues, + }; +} + +// ---------------------------------------------------------------------------- + +export function extend(publicAPI, model, initialValues = {}) { + vtkLocator.extend(publicAPI, model, defaultValues(initialValues)); + + // Make this a VTK object + macro.obj(publicAPI, model); + + macro.get(publicAPI, model, ['numberOfBuckets']); + + macro.setGetArray(publicAPI, model, ['bounds']); + + // Object specific methods + vtkAbstractPointLocator(publicAPI, model); +} + +// ---------------------------------------------------------------------------- + +export default { extend }; diff --git a/Sources/Common/DataModel/AbstractPointLocator/test/testLocator.js b/Sources/Common/DataModel/AbstractPointLocator/test/testLocator.js new file mode 100644 index 00000000000..9ba2d00ce07 --- /dev/null +++ b/Sources/Common/DataModel/AbstractPointLocator/test/testLocator.js @@ -0,0 +1,11 @@ +import test from 'tape-catch'; +import vtkAbstractPointLocator from 'vtk.js/Sources/Common/DataModel/AbstractPointLocator'; + +test('Test vtkAbstractPointLocator instance', (t) => { + t.ok(vtkAbstractPointLocator, 'Make sure the class definition exists'); + t.ok( + vtkAbstractPointLocator.newInstance === undefined, + 'Make sure class is abstract' + ); + t.end(); +}); diff --git a/Sources/Common/DataModel/IncrementalOctreeNode/index.d.ts b/Sources/Common/DataModel/IncrementalOctreeNode/index.d.ts new file mode 100644 index 00000000000..23f3bca87f3 --- /dev/null +++ b/Sources/Common/DataModel/IncrementalOctreeNode/index.d.ts @@ -0,0 +1,285 @@ +import { vtkAlgorithm, vtkObject } from "../../../interfaces"; +import { Bounds, Vector3 } from "../../../types"; +import vtkCellArray from "vtk.js/Sources/Common/DataModel/CellArray"; +import vtkPolyData from "vtk.js/Sources/Common/DataModel/PolyData"; +import vtkPoints from "vtk.js/Sources/Common/Core/Points"; + +/** + * + */ +export interface IIncrementalOctreeNodeInitialValues { + pointIdSet?: number[]; + minBounds?: Bounds; + maxBounds?: Bounds; + minDataBounds?: Bounds; + maxDataBounds?: Bounds; + parent?: vtkIncrementalOctreeNode; + children?: vtkIncrementalOctreeNode[]; +} + +export interface vtkIncrementalOctreeNode extends vtkObject { + /** + * Create a list object for storing point indices. + */ + createPointIdSet(): void; + + /** + * Set the spatial bounding box of the node. This function sets a default + * data bounding box. + * + * @param {Number} x1 + * @param {Number} x2 + * @param {Number} y1 + * @param {Number} y2 + * @param {Number} z1 + * @param {Number} z2 + */ + setBounds( + x1: number, + x2: number, + y1: number, + y2: number, + z1: number, + z2: number + ): void; + + /** + * Get the spatial bounding box of the node. The values are returned via + * an array in order of: x_min, x_max, y_min, y_max, z_min, z_max. + * + * @param {Bounds} bounds + */ + getBounds(bounds: Bounds): void; + + /** + * Given a point inserted to either this node (a leaf node) or a descendant + * leaf (of this node --- when this node is a non-leaf node), update the + * counter and the data bounding box for this node only. The data bounding box + * is considered only if updateData is non-zero. The returned value indicates + * whether (1) or not (0) the data bounding box is actually updated. Note that + * argument nHits must be 1 unless this node is updated with a number (nHits) + * of exactly duplicate points as a whole via a single call to this function. + * + * @param {Vector3} point + * @param {Number} nHits + * @param {Boolean} updateData + */ + updateCounterAndDataBounds( + point: Vector3, + nHits: number, + updateData: boolean + ): boolean; + + /** + * Given a point inserted to either this node (a leaf node) or a descendant + * leaf (of this node --- when this node is a non-leaf node), update the + * counter and the data bounding box recursively bottom-up until a specified + * node. The data bounding box is considered only if updateData is non-zero. + * The returned value indicates whether (true) or not (false) the data bounding box + * is actually updated. Note that argument nHits must be 1 unless this node + * is updated with a number (nHits) of exactly duplicate points as a whole + * via a single call to this function. + * + * @param {Vector3} point + * @param {Number} nHits + * @param {Boolean} updateData + * @param {vtkIncrementalOctreeNode} endNode + */ + updateCounterAndDataBoundsRecursively( + point: Vector3, + nHits: number, + updateData: boolean, + endNode: vtkIncrementalOctreeNode + ): boolean; + + /** + * Given a point, determine whether (true) or not (false) it is an exact duplicate + * of all the points, if any, maintained in this node. In other words, to + * check if this given point and all existing points, if any, of this node + * are exactly duplicate with one another. + * + * @param {Vector3} point + */ + containsDuplicatePointsOnly(point: Vector3): boolean; + + /** + * Determine whether or not this node is a leaf. + */ + isLeaf(): boolean; + + /** + * Get the child at the given index i. + * i must be an int between 0 and 7. + * + * @param {Number} i + */ + getChild(i: number): vtkIncrementalOctreeNode; + + /** + * Given a number (>= threshold) of all exactly duplicate points (accessible + * via points and pntIds, but with exactly the same 3D coordinate) maintained + * in this leaf node and a point (absolutely not a duplicate any more, with + * pntIdx storing the index in points)) to be inserted to this node, separate + * all the duplicate points from this new point by means of usually recursive + * node sub-division such that the former points are inserted to a descendant + * leaf while the new point is inserted to a sibling of this descendant leaf. + * Argument ptMode specifies whether the point is not inserted at all but only + * the point index is provided upon 0, the point is inserted via vtkPoints:: + * InsertPoint() upon 1, or this point is instead inserted through vtkPoints:: + * InsertNextPoint() upon 2. + * + * @param {vtkPoints} points + * @param {Number[]} pntIds + * @param {Vector3} newPnt + * @param {Number} pntIdx + * @param {Number} maxPts + * @param {Number} ptMode + */ + separateExactlyDuplicatePointsFromNewInsertion( + points: vtkPoints, + pntIds: number[], + newPnt: Vector3, + pntIdx: number, + maxPts: number, + ptMode: number + ): void; + + /** + * Divide this LEAF node into eight child nodes as the number of points + * maintained by this leaf node has reached the threshold maxPts while + * another point newPnt is just going to be inserted to it. The available + * point-indices pntIds are distributed to the child nodes based on the + * point coordinates (available through points). Note that this function + * can incur recursive node-division to determine the specific leaf node + * for accepting the new point (with pntIdx storing the index in points) + * because the existing maxPts points may fall within only one of the eight + * child nodes to make a radically imbalanced layout within the node (to + * be divided). Argument ptMode specifies whether the point is not inserted + * at all but instead only the point index is provided upon 0, the point is + * inserted via vtkPoints.InsertPoint() upon 1, or the point is inserted by + * vtkPoints.InsertNextPoint() upon 2. The returned value of this function + * indicates whether pntIds needs to be destroyed (1) or just unregistered + * from this node as it has been attached to another node (0). + * numberOfNodes in the tree is updated with new created nodes + * + * @param {vtkPoints} points + * @param {Number[]} pntIds + * @param {Vector3} newPnt + * @param {Number} pntIdx + * @param {Number} maxPts + * @param {Number} ptMode + * @param {Number} numberOfNodes + */ + createChildNodes( + points: vtkPoints, + pntIds: number[], + newPnt: Vector3, + pntIdx: number, + maxPts: number, + ptMode: number, + numberOfNodes: number + ): { success: boolean; numberOfNodes: number; pointIdx: number }; + + /** + * This function is called after a successful point-insertion check and + * only applies to a leaf node. Prior to a call to this function, the + * octree should have been retrieved top-down to find the specific leaf + * node in which this new point (newPt) will be inserted. The actual index + * of the new point (to be inserted to points) is stored in pntId. Argument + * ptMode specifies whether the point is not inserted at all but instead only + * the point index is provided upon 0, the point is inserted via vtkPoints. + * insertPoint() upon 1, or it is inserted via vtkPoints.insertNextPoint() + * upon 2. For case 0, pntId needs to be specified. For cases 1 and 2, the + * actual point index is returned via pntId. Note that this function always + * returns 1 to indicate the success of point insertion. + * numberOfNodes is the number of nodes present in the tree at this time. + * it is used to assign an ID to each node which can be used to associate + * application specific information with each node. It is updated if new nodes + * are added to the tree. + * + * @param {Number} points + * @param {Number} newPnt + * @param {Number} maxPts + * @param {Number} pntId + * @param {Number} ptMode + * @param {Number} numberOfNodes + */ + insertPoint( + points: number, + newPnt: number, + maxPts: number, + pntId: number, + ptMode: number, + numberOfNodes: number + ): { numberOfNodes: number; pointIdx: number }; + + /** + * Compute the minimum squared distance from a point to this node, with all + * six boundaries considered. The data bounding box is checked if checkData + * is non-zero. The closest on-boundary point is returned via closest. + * + * @param {Vector3} point + * @param {Vector3} closest + * @param {Boolean} innerOnly + * @param {vtkIncrementalOctreeNode} rootNode + * @param {Boolean} checkData + * @returns {Number} + */ + getDistance2ToBoundary( + point: Vector3, + closest: Vector3, + innerOnly: boolean, + rootNode: vtkIncrementalOctreeNode, + checkData: boolean + ): number; + + /** + * Given a point inside this node, get the minimum squared distance to all + * inner boundaries. An inner boundary is a node's face that is shared by + * another non-root node. + * + * @param {Vector3} point + * @param {vtkIncrementalOctreeNode} rootNode + */ + getDistance2ToInnerBoundary( + point: Vector3, + rootNode: vtkIncrementalOctreeNode + ): number; +} + +// ---------------------------------------------------------------------------- +// Static API +// ---------------------------------------------------------------------------- + +/** + * Method use to decorate a given object (publicAPI+model) with vtkIncrementalOctreeNode characteristics. + * + * @param publicAPI object on which methods will be bounds (public) + * @param model object on which data structure will be bounds (protected) + * @param {object} [initialValues] (default: {}) + */ +export function extend( + publicAPI: object, + model: object, + initialValues?: IIncrementalOctreeNodeInitialValues +): void; + +// ---------------------------------------------------------------------------- + +/** + * Method use to create a new instance of vtkIncrementalOctreeNode + * @param {IIncrementalOctreeNodeInitialValues} [initialValues] for pre-setting some of its content + */ +export function newInstance( + initialValues?: IIncrementalOctreeNodeInitialValues +): vtkIncrementalOctreeNode; + +/** + * vtkIncrementalOctreeNode + */ +export declare const vtkIncrementalOctreeNode: { + newInstance: typeof newInstance; + extend: typeof extend; +}; + +export default vtkIncrementalOctreeNode; diff --git a/Sources/Common/DataModel/IncrementalOctreeNode/index.js b/Sources/Common/DataModel/IncrementalOctreeNode/index.js new file mode 100644 index 00000000000..3566e70783b --- /dev/null +++ b/Sources/Common/DataModel/IncrementalOctreeNode/index.js @@ -0,0 +1,804 @@ +import macro from 'vtk.js/Sources/macros'; +import * as vtkMath from 'vtk.js/Sources/Common/Core/Math'; + +const { vtkErrorMacro } = macro; + +const OCTREENODE_INSERTPOINT = [ + (points, pointIdx, coords) => pointIdx, + (points, pointIdx, coords) => { + points.setTuple(pointIdx, coords); + return pointIdx; + }, + (points, pointIdx, coords) => points.insertNextTuple(coords), +]; + +// Given the index (0 ~ 7) of a child node, the spatial bounding axis (0 ~ 2 +// for x, y, and z), and the value (0 ~ 1 for min and max) to access, this LUT +// allows for rapid assignment of its spatial bounding box --- MinBounds[3] +// and MaxBounds[3], with each specific value or entry of this LUT pointing to +// MinBounds[3] for 0, center point for 1, or MaxBounds[3] for 2. +const OCTREE_CHILD_BOUNDS_LUT = [ + [ + [0, 1], + [0, 1], + [0, 1], + ], + [ + [1, 2], + [0, 1], + [0, 1], + ], + [ + [0, 1], + [1, 2], + [0, 1], + ], + [ + [1, 2], + [1, 2], + [0, 1], + ], + + [ + [0, 1], + [0, 1], + [1, 2], + ], + [ + [1, 2], + [0, 1], + [1, 2], + ], + [ + [0, 1], + [1, 2], + [1, 2], + ], + [ + [1, 2], + [1, 2], + [1, 2], + ], +]; + +function vtkIncrementalOctreeNode(publicAPI, model) { + // Set our className + model.classHierarchy.push('vtkIncrementalOctreeNode'); + + //------------------------------------------------------------------------------ + publicAPI.createPointIdSet = () => { + if (model.pointIdSet == null) { + model.pointIdSet = []; + // model.pointIdSet.allocate(initSize, growSize); + } + }; + + //------------------------------------------------------------------------------ + publicAPI.setBounds = (x1, x2, y1, y2, z1, z2) => { + model.minBounds[0] = x1; + model.maxBounds[0] = x2; + model.minBounds[1] = y1; + model.maxBounds[1] = y2; + model.minBounds[2] = z1; + model.maxBounds[2] = z2; + + model.minDataBounds[0] = x2; + model.maxDataBounds[0] = x1; + model.minDataBounds[1] = y2; + model.maxDataBounds[1] = y1; + model.minDataBounds[2] = z2; + model.maxDataBounds[2] = z1; + }; + + //------------------------------------------------------------------------------ + publicAPI.getBounds = (bounds) => { + bounds[0] = model.minBounds[0]; + bounds[1] = model.maxBounds[0]; + bounds[2] = model.minBounds[1]; + bounds[3] = model.maxBounds[1]; + bounds[4] = model.minBounds[2]; + bounds[5] = model.maxBounds[2]; + }; + + //------------------------------------------------------------------------------ + publicAPI.updateCounterAndDataBounds = (point, nHits, updateData) => { + model.numberOfPoints += nHits; + + if (!updateData) return false; + + let updated = false; + + if (point[0] < model.minDataBounds[0]) { + updated = true; + model.minDataBounds[0] = point[0]; + } + if (point[0] > model.maxDataBounds[0]) { + updated = true; + model.maxDataBounds[0] = point[0]; + } + + if (point[1] < model.minDataBounds[1]) { + updated = true; + model.minDataBounds[1] = point[1]; + } + if (point[1] > model.maxDataBounds[1]) { + updated = true; + model.maxDataBounds[1] = point[1]; + } + + if (point[2] < model.minDataBounds[2]) { + updated = true; + model.minDataBounds[2] = point[2]; + } + if (point[2] > model.maxDataBounds[2]) { + updated = true; + model.maxDataBounds[2] = point[2]; + } + + return updated; + }; + + //------------------------------------------------------------------------------ + publicAPI.updateCounterAndDataBoundsRecursively = ( + point, + nHits, + updateData, + endNode + ) => { + const updated = publicAPI.updateCounterAndDataBounds( + point, + nHits, + updateData + ); + + return model.parent === endNode + ? updated + : model.parent.updateCounterAndDataBoundsRecursively( + point, + nHits, + updated, + endNode + ); + }; + + //------------------------------------------------------------------------------ + publicAPI.containsDuplicatePointsOnly = (point) => + model.minDataBounds[0] === point[0] && + point[0] === model.maxDataBounds[0] && + model.minDataBounds[1] === point[1] && + point[1] === model.maxDataBounds[1] && + model.minDataBounds[2] === point[2] && + point[2] === model.maxDataBounds[2]; + + //------------------------------------------------------------------------------ + publicAPI.isLeaf = () => model.children == null; + + //------------------------------------------------------------------------------ + publicAPI.getChild = (i) => model.children[i]; + + //------------------------------------------------------------------------------ + publicAPI.separateExactlyDuplicatePointsFromNewInsertion = ( + points, + pntIds, + newPnt, + pntIdx, + maxPts, + ptMode + ) => { + // the number of points already maintained in this leaf node + // >= maxPts AND all of them are exactly duplicate with one another + // BUT the new point is not a duplicate of them any more + let pointIdx = pntIdx; + let i; + const dupPnt = [0.0, 0.0, 0.0]; + const octMin = [0.0, 0.0, 0.0]; + const octMid = [0.0, 0.0, 0.0]; + const octMax = [0.0, 0.0, 0.0]; + // TODO: Check + // double* boxPtr[3] = { null, null, null }; + // vtkIncrementalOctreeNode* ocNode = null; + // vtkIncrementalOctreeNode* duplic = this; + // vtkIncrementalOctreeNode* single = this; + const boxPtr = [null, null, null]; + let ocNode = null; + let duplic = publicAPI; + let single = publicAPI; + + // the coordinate of the duplicate points: note pntIds == model.pointIdSet + points.getPoint(pntIds[0], dupPnt); + + while (duplic === single) { + // as long as separation has not been achieved + // update the current (in recursion) node and access the bounding box info + ocNode = duplic; + octMid[0] = (ocNode.minBounds[0] + ocNode.maxBounds[0]) * 0.5; + octMid[1] = (ocNode.minBounds[1] + ocNode.maxBounds[1]) * 0.5; + octMid[2] = (ocNode.minBounds[2] + ocNode.maxBounds[2]) * 0.5; + boxPtr[0] = ocNode.minBounds; + boxPtr[1] = octMid; + boxPtr[2] = ocNode.maxBounds; + + // create eight child nodes + ocNode.children = [ + vtkIncrementalOctreeNode.newInstance(), + vtkIncrementalOctreeNode.newInstance(), + vtkIncrementalOctreeNode.newInstance(), + vtkIncrementalOctreeNode.newInstance(), + vtkIncrementalOctreeNode.newInstance(), + vtkIncrementalOctreeNode.newInstance(), + vtkIncrementalOctreeNode.newInstance(), + vtkIncrementalOctreeNode.newInstance(), + ]; + for (i = 0; i < 8; i++) { + // x-bound: axis 0 + octMin[0] = boxPtr[OCTREE_CHILD_BOUNDS_LUT[i][0][0]][0]; + octMax[0] = boxPtr[OCTREE_CHILD_BOUNDS_LUT[i][0][1]][0]; + + // y-bound: axis 1 + octMin[1] = boxPtr[OCTREE_CHILD_BOUNDS_LUT[i][1][0]][1]; + octMax[1] = boxPtr[OCTREE_CHILD_BOUNDS_LUT[i][1][1]][1]; + + // z-bound: axis 2 + octMin[2] = boxPtr[OCTREE_CHILD_BOUNDS_LUT[i][2][0]][2]; + octMax[2] = boxPtr[OCTREE_CHILD_BOUNDS_LUT[i][2][1]][2]; + + ocNode.children[i] = vtkIncrementalOctreeNode.newInstance(); + ocNode.children[i].setParent(ocNode); + ocNode.children[i].setBounds( + octMin[0], + octMax[0], + octMin[1], + octMax[1], + octMin[2], + octMax[2] + ); + } + + // determine the leaf node of the duplicate points & that of the new point + duplic = ocNode.children[ocNode.getChildIndex(dupPnt)]; + single = ocNode.children[ocNode.getChildIndex(newPnt)]; + } + // boxPtr[0] = null; + // boxPtr[1] = null; + // boxPtr[2] = null; + + // Now the duplicate points have been separated from the new point // + + // create a vtkIdList object for the new point + // update the counter and the data bounding box until the root node + // (including the root node) + pointIdx = OCTREENODE_INSERTPOINT[ptMode](points, pointIdx, newPnt); + // eslint-disable-next-line no-bitwise + // TODO: createPointIdSet is called with initSize and growSize in C++. + // single.createPointIdSet(maxPts >> 2, maxPts >> 1); + single.createPointIdSet(); + single.getPointIdSet().push(pointIdx); + single.updateCounterAndDataBoundsRecursively(newPnt, 1, 1, null); + + // We just need to reference pntIds while un-registering it from 'this'. + // This avoids deep-copying point ids from pntIds to duplic's PointIdSet. + // update the counter and the data bounding box, but until 'this' node + // (excluding 'this' node) + duplic.setPointIdSet(pntIds); + duplic.updateCounterAndDataBoundsRecursively( + dupPnt, + pntIds.length, + 1, + publicAPI + ); + return pointIdx; + }; + + //------------------------------------------------------------------------------ + publicAPI.createChildNodes = ( + points, + pntIds, + newPnt, + pntIdx, + maxPts, + ptMode, + numberOfNodes + ) => { + // There are two scenarios for which this function is invoked. + // + // (1) the number of points already maintained in this leaf node + // == maxPts AND not all of them are exactly duplicate + // AND the new point is not a duplicate of them all + // (2) the number of points already maintained in this leaf node + // >= maxPts AND all of them are exactly duplicate with one another + // BUT the new point is not a duplicate of them any more + + // address case (2) first if necessary + let nbNodes = numberOfNodes; + let pointIdx = pntIdx; + const sample = []; + points.getPoint(pntIds[0], sample); + if (publicAPI.containsDuplicatePointsOnly(sample)) { + pointIdx = publicAPI.separateExactlyDuplicatePointsFromNewInsertion( + points, + pntIds, + newPnt, + pointIdx, + maxPts, + ptMode + ); + return { success: false, nbNodes, pointIdx }; + } + + // then address case (1) below + + let i; + let target; + let dvidId = -1; // index of the sub-dividing octant, if any + let fullId = -1; // index of the full octant, if any + const numIds = [0, 0, 0, 0, 0, 0, 0, 0]; + const octMin = []; + const octMax = []; + const tempPt = []; + let tempId; + + const octMid = [ + (model.minBounds[0] + model.maxBounds[0]) * 0.5, + (model.minBounds[1] + model.maxBounds[1]) * 0.5, + (model.minBounds[2] + model.maxBounds[2]) * 0.5, + ]; + const boxPtr = [model.minBounds, octMid, model.maxBounds]; + + // create eight child nodes + model.children = []; + for (i = 0; i < 8; i++) { + // x-bound: axis 0 + octMin[0] = boxPtr[OCTREE_CHILD_BOUNDS_LUT[i][0][0]][0]; + octMax[0] = boxPtr[OCTREE_CHILD_BOUNDS_LUT[i][0][1]][0]; + + // y-bound: axis 1 + octMin[1] = boxPtr[OCTREE_CHILD_BOUNDS_LUT[i][1][0]][1]; + octMax[1] = boxPtr[OCTREE_CHILD_BOUNDS_LUT[i][1][1]][1]; + + // z-bound: axis 2 + octMin[2] = boxPtr[OCTREE_CHILD_BOUNDS_LUT[i][2][0]][2]; + octMax[2] = boxPtr[OCTREE_CHILD_BOUNDS_LUT[i][2][1]][2]; + + // This call internally sets the cener and default data bounding box, too. + model.children[i] = vtkIncrementalOctreeNode.newInstance(); + model.children[i].iD = nbNodes++; + // TODO: Check + // model.children[i].setParent(this); + model.children[i].setParent(model); + model.children[i].setBounds( + octMin[0], + octMax[0], + octMin[1], + octMax[1], + octMin[2], + octMax[2] + ); + + // allocate a list of point-indices (size = 2^n) for index registration + // eslint-disable-next-line no-bitwise + // model.children[i].createPointIdSet(maxPts >> 2, maxPts >> 1); + model.children[i].createPointIdSet(); + } + boxPtr[0] = null; + boxPtr[1] = null; + boxPtr[2] = null; + + // distribute the available point-indices to the eight child nodes + for (i = 0; i < maxPts; i++) { + tempId = pntIds[i]; + points.getPoint(tempId, tempPt); + target = publicAPI.getChildIndex(tempPt); + model.children[target].getPointIdSet().insertNextId(tempId); + model.children[target].updateCounterAndDataBounds(tempPt); + numIds[target]++; + } + + // locate the full child, just if any + for (i = 0; i < 8; i++) { + if (numIds[i] === maxPts) { + fullId = i; + break; + } + } + + target = publicAPI.getChildIndex(newPnt); + if (fullId === target) { + // The fact is that we are going to insert the new point to an already + // full octant (child node). Thus we need to further divide this child + // to avoid the overflow problem. + // TODO + ({ numberOfNodes: nbNodes, pointIdx } = model.children[ + target + ].createChildNodes( + points, + pntIds, + newPnt, + pointIdx, + maxPts, + ptMode, + nbNodes + )); + dvidId = fullId; + } else { + // the initial division is a success + pointIdx = OCTREENODE_INSERTPOINT[ptMode](points, pointIdx, newPnt); + model.children[target].getPointIdSet().push(pointIdx); + model.children[target].updateCounterAndDataBoundsRecursively( + newPnt, + 1, + 1, + null + ); + + // NOTE: The counter below might reach the threshold, though we delay the + // sub-division of this child node until the next point insertion occurs. + numIds[target]++; + } + + // Now it is time to reclaim those un-used vtkIdList objects, of which each + // either is empty or still needs to be deleted due to further division of + // the child node. This post-deallocation of the un-used vtkIdList objects + // (of some child nodes) is based on the assumption that retrieving the + // 'maxPts' points from vtkPoints and the associated 'maxPts' point-indices + // from vtkIdList is more expensive than reclaiming at most 8 vtkIdList + // objects at hand. + for (i = 0; i < 8; i++) { + if (numIds[i] === 0 || i === dvidId) { + model.children[i].deletePointIdSet(); + } + } + + // notify vtkIncrementalOctreeNode::InsertPoint() to destroy pntIds + return { success: true, numberOfNodes: nbNodes, pointIdx }; + }; + + //------------------------------------------------------------------------------ + publicAPI.insertPoint = ( + points, + newPnt, + maxPts, + pntId, + ptMode, + numberOfNodes + ) => { + let success; + let nbNodes = 0; + let pointIdx = pntId; + + if (model.pointIdSet) { + // there has been at least one point index + if ( + model.pointIdSet.length < maxPts || + publicAPI.containsDuplicatePointsOnly(newPnt) + ) { + // this leaf node is not full or + // this leaf node is full, but of all exactly duplicate points + // and the point under check is another duplicate of these points + pointIdx = OCTREENODE_INSERTPOINT[ptMode](points, pointIdx, newPnt); + // model.pointIdSet.insertNextId(*pntId); + model.pointIdSet.push(pointIdx); + publicAPI.updateCounterAndDataBoundsRecursively(newPnt, 1, 1, null); + } else { + // overflow: divide this node and delete the list of point-indices. + // Note that the number of exactly duplicate points might be greater + // than or equal to maxPts. + // TODO: pointIdx + ({ success, numberOfNodes: nbNodes } = publicAPI.createChildNodes( + points, + model.pointIdSet, + newPnt, + pointIdx, + maxPts, + ptMode, + numberOfNodes + )); + if (success) { + // TODO: Check + // model.pointIdSet.delete(); + model.pointIdSet = null; + } + model.pointIdSet = null; + } + } else { + // there has been no any point index registered in this leaf node + pointIdx = OCTREENODE_INSERTPOINT[ptMode](points, pointIdx, newPnt); + // model.pointIdSet = vtkIdList::New(); + // model.pointIdSet.allocate((maxPts >> 2), (maxPts >> 1)); + // model.pointIdSet.insertNextId(*pntId); + model.pointIdSet = []; + model.pointIdSet.push(pointIdx); + publicAPI.updateCounterAndDataBoundsRecursively(newPnt, 1, 1, null); + } + + return { numberOfNodes: numberOfNodes + nbNodes, pointIdx }; + }; + + //------------------------------------------------------------------------------ + publicAPI.getDistance2ToBoundary = ( + point, + closest, + innerOnly, + rootNode, + checkData + ) => { + // It is mandatory that GetMinDataBounds() and GetMaxDataBounds() be used. + // Direct access to MinDataBounds and MaxDataBounds might incur problems. + let thisMin = null; + let thisMax = null; + let rootMin = null; + let rootMax = null; + // TODO: Check + // let minDist = VTK_DOUBLE_MAX; + let minDist = Number.MAX_VALUE; // minimum distance to the boundaries + if (checkData) { + thisMin = publicAPI.getMinDataBounds(); + thisMax = publicAPI.getMaxDataBounds(); + rootMin = rootNode.getMinDataBounds(); + rootMax = rootNode.getMaxDataBounds(); + } else { + thisMin = model.minBounds; + thisMax = model.maxBounds; + rootMin = rootNode.getMinBounds(); + rootMax = rootNode.getMaxBounds(); + } + + let minFace = 0; // index of the face with min distance to the point + const beXless = Number(point[0] < thisMin[0]); + const beXmore = Number(point[0] > thisMax[0]); + const beYless = Number(point[1] < thisMin[1]); + const beYmore = Number(point[1] > thisMax[1]); + const beZless = Number(point[2] < thisMin[2]); + const beZmore = Number(point[2] > thisMax[2]); + const withinX = Number(!beXless && !beXmore); + const withinY = Number(!beYless && !beYmore); + const withinZ = Number(!beZless && !beZmore); + // eslint-disable-next-line no-bitwise + const xyzFlag = (withinZ << 2) + (withinY << 1) + withinX; + + switch (xyzFlag) { + case 0: { + // withinZ = 0; withinY = 0; withinX = 0 + // closest to a corner + + closest[0] = beXless ? thisMin[0] : thisMax[0]; + closest[1] = beYless ? thisMin[1] : thisMax[1]; + closest[2] = beZless ? thisMin[2] : thisMax[2]; + minDist = vtkMath.distance2BetweenPoints(point, closest); + break; + } + + case 1: { + // withinZ = 0; withinY = 0; withinX = 1 + // closest to an x-aligned edge + + closest[0] = point[0]; + closest[1] = beYless ? thisMin[1] : thisMax[1]; + closest[2] = beZless ? thisMin[2] : thisMax[2]; + minDist = vtkMath.distance2BetweenPoints(point, closest); + break; + } + + case 2: { + // withinZ = 0; withinY = 1; withinX = 0 + // closest to a y-aligned edge + + closest[0] = beXless ? thisMin[0] : thisMax[0]; + closest[1] = point[1]; + closest[2] = beZless ? thisMin[2] : thisMax[2]; + minDist = vtkMath.distance2BetweenPoints(point, closest); + break; + } + + case 3: { + // withinZ = 0; withinY = 1; withinX = 1 + // closest to a z-face + + if (beZless) { + minDist = thisMin[2] - point[2]; + closest[2] = thisMin[2]; + } else { + minDist = point[2] - thisMax[2]; + closest[2] = thisMax[2]; + } + + minDist *= minDist; + closest[0] = point[0]; + closest[1] = point[1]; + break; + } + + case 4: { + // withinZ = 1; withinY = 0; withinX = 0 + // cloest to a z-aligned edge + + closest[0] = beXless ? thisMin[0] : thisMax[0]; + closest[1] = beYless ? thisMin[1] : thisMax[1]; + closest[2] = point[2]; + minDist = vtkMath.distance2BetweenPoints(point, closest); + break; + } + + case 5: { + // withinZ = 1; withinY = 0; withinX = 1 + // closest to a y-face + + if (beYless) { + minDist = thisMin[1] - point[1]; + closest[1] = thisMin[1]; + } else { + minDist = point[1] - thisMax[1]; + closest[1] = thisMax[1]; + } + + minDist *= minDist; + closest[0] = point[0]; + closest[2] = point[2]; + break; + } + + case 6: { + // withinZ = 1; withinY = 1; withinX = 0 + // closest to an x-face + + if (beXless) { + minDist = thisMin[0] - point[0]; + closest[0] = thisMin[0]; + } else { + minDist = point[0] - thisMax[0]; + closest[0] = thisMax[0]; + } + + minDist *= minDist; + closest[1] = point[1]; + closest[2] = point[2]; + break; + } + + case 7: { + // withinZ = 1; withinY = 1; withinZ = 1 + // point is inside the box + + if (innerOnly) { + // check only inner boundaries + let faceDst; + + faceDst = point[0] - thisMin[0]; // x-min face + if (thisMin[0] !== rootMin[0] && faceDst < minDist) { + minFace = 0; + minDist = faceDst; + } + + faceDst = thisMax[0] - point[0]; // x-max face + if (thisMax[0] !== rootMax[0] && faceDst < minDist) { + minFace = 1; + minDist = faceDst; + } + + faceDst = point[1] - thisMin[1]; // y-min face + if (thisMin[1] !== rootMin[1] && faceDst < minDist) { + minFace = 2; + minDist = faceDst; + } + + faceDst = thisMax[1] - point[1]; // y-max face + if (thisMax[1] !== rootMax[1] && faceDst < minDist) { + minFace = 3; + minDist = faceDst; + } + + faceDst = point[2] - thisMin[2]; // z-min face + if (thisMin[2] !== rootMin[2] && faceDst < minDist) { + minFace = 4; + minDist = faceDst; + } + + faceDst = thisMax[2] - point[2]; // z-max face + if (thisMax[2] !== rootMax[2] && faceDst < minDist) { + minFace = 5; + minDist = faceDst; + } + } else { + // check all boundaries + const tmpDist = [ + point[0] - thisMin[0], + thisMax[0] - point[0], + point[1] - thisMin[1], + thisMax[1] - point[1], + point[2] - thisMin[2], + thisMax[2] - point[2], + ]; + + for (let i = 0; i < 6; i++) { + if (tmpDist[i] < minDist) { + minFace = i; + minDist = tmpDist[i]; + } + } + } + + // no square operation if no any inner boundary + if (minDist !== Number.MAX_VALUE) { + minDist *= minDist; + } + + closest[0] = point[0]; + closest[1] = point[1]; + closest[2] = point[2]; + + // minFace: the quad with the min distance to the point + // 0: x-min face ===> xyzIndx = 0: x and minFace & 1 = 0: thisMin + // 1: x-max face ===> xyzIndx = 0: x and minFace & 1 = 1: thisMax + // 2: y-min face ===> xyzIndx = 1: y and minFace & 1 = 0: thisMin + // 3: y-max face ===> xyzIndx = 1: y and minFace & 1 = 1: thisMax + // 4: z-min face ===> xyzIndx = 2: z and minFace & 1 = 0: thisMin + // 5: z-max face ===> xyzIndx = 2: z and minFace & 1 = 1: thisMax + const pMinMax = [thisMin, thisMax]; + // eslint-disable-next-line no-bitwise + const xyzIndx = minFace >> 1; + // eslint-disable-next-line no-bitwise + closest[xyzIndx] = pMinMax[minFace & 1][xyzIndx]; + + break; + } + default: + vtkErrorMacro('unexpected case in getDistance2ToBoundary'); + } + + return minDist; + }; + + //------------------------------------------------------------------------------ + publicAPI.getDistance2ToInnerBoundary = (point, rootNode) => { + const dummy = []; + return publicAPI.getDistance2ToBoundary(point, dummy, 0, rootNode, 0); + }; +} + +// ---------------------------------------------------------------------------- +// Object factory +// ---------------------------------------------------------------------------- + +const DEFAULT_VALUES = { + pointIdSet: null, + minBounds: [], + maxBounds: [], + minDataBounds: [], + maxDataBounds: [], + parent: null, + children: null, +}; + +// ---------------------------------------------------------------------------- + +export function extend(publicAPI, model, initialValues = {}) { + Object.assign(model, DEFAULT_VALUES, initialValues); + + // Make this a VTK object + macro.obj(publicAPI, model); + + macro.setGet(publicAPI, model, [ + 'minBounds', + 'maxBounds', + 'minDataBounds', + 'maxDataBounds', + ]); + + macro.get(publicAPI, model, ['pointIdSet']); + + // TODO: No get? + macro.set(publicAPI, model, ['parent']); + + // Object specific methods + vtkIncrementalOctreeNode(publicAPI, model); +} + +// ---------------------------------------------------------------------------- + +export const newInstance = macro.newInstance( + extend, + 'vtkIncrementalOctreeNode' +); + +// ---------------------------------------------------------------------------- + +export default { newInstance, extend }; diff --git a/Sources/Common/DataModel/IncrementalOctreeNode/test/testIncrementalOctreeNode.js b/Sources/Common/DataModel/IncrementalOctreeNode/test/testIncrementalOctreeNode.js new file mode 100644 index 00000000000..d244cb6538f --- /dev/null +++ b/Sources/Common/DataModel/IncrementalOctreeNode/test/testIncrementalOctreeNode.js @@ -0,0 +1,9 @@ +import test from 'tape-catch'; +import vtkIncrementalOctreeNode from 'vtk.js/Sources/Common/DataModel/IncrementalOctreeNode'; + +test('Test vtkIncrementalOctreeNode instance', (t) => { + t.ok(vtkIncrementalOctreeNode, 'Make sure the class definition exists'); + const instance = vtkIncrementalOctreeNode.newInstance(); + t.ok(instance); + t.end(); +}); diff --git a/Sources/Common/DataModel/IncrementalOctreePointLocator/index.d.ts b/Sources/Common/DataModel/IncrementalOctreePointLocator/index.d.ts new file mode 100644 index 00000000000..ef54954bb82 --- /dev/null +++ b/Sources/Common/DataModel/IncrementalOctreePointLocator/index.d.ts @@ -0,0 +1,64 @@ +import { vtkAlgorithm, vtkObject } from "../../../interfaces"; +import { Vector3 } from "../../../types"; +import vtkCellArray from "vtk.js/Sources/Common/DataModel/CellArray"; +import vtkPolyData from "vtk.js/Sources/Common/DataModel/PolyData"; +import vtkIncrementalOctreeNode from "../IncrementalOctreeNode"; +import vtkPoints from "../../Core/Points"; +import { IAbstractPointLocatorInitialValues } from "../AbstractPointLocator"; + +/** + * + */ +export interface IIncrementalOctreePointLocatorInitialValues + extends IAbstractPointLocatorInitialValues { + fudgeFactor: number; + octreeMaxDimSize: number; + buildCubicOctree: boolean; + maxPointsPerLeaf: number; + insertTolerance2: number; + locatorPoints: vtkPoints; + octreeRootNode: vtkIncrementalOctreeNode; + numberOfNodes: number; +} + +type vtkIncrementalOctreePointLocatorBase = vtkObject; + +export interface vtkIncrementalOctreePointLocator + extends vtkIncrementalOctreePointLocatorBase {} + +// ---------------------------------------------------------------------------- +// Static API +// ---------------------------------------------------------------------------- + +/** + * Method use to decorate a given object (publicAPI+model) with vtkIncrementalOctreePointLocator characteristics. + * + * @param publicAPI object on which methods will be bounds (public) + * @param model object on which data structure will be bounds (protected) + * @param {object} [initialValues] (default: {}) + */ +export function extend( + publicAPI: object, + model: object, + initialValues?: IIncrementalOctreePointLocatorInitialValues +): void; + +// ---------------------------------------------------------------------------- + +/** + * Method use to create a new instance of vtkIncrementalOctreePointLocator + * @param {IIncrementalOctreePointLocatorInitialValues} [initialValues] for pre-setting some of its content + */ +export function newInstance( + initialValues?: IIncrementalOctreePointLocatorInitialValues +): vtkIncrementalOctreePointLocator; + +/** + * vtkIncrementalOctreePointLocator + */ +export declare const vtkIncrementalOctreePointLocator: { + newInstance: typeof newInstance; + extend: typeof extend; +}; + +export default vtkIncrementalOctreePointLocator; diff --git a/Sources/Common/DataModel/IncrementalOctreePointLocator/index.js b/Sources/Common/DataModel/IncrementalOctreePointLocator/index.js new file mode 100644 index 00000000000..8c56d27b461 --- /dev/null +++ b/Sources/Common/DataModel/IncrementalOctreePointLocator/index.js @@ -0,0 +1,391 @@ +import macro from 'vtk.js/Sources/macros'; +import vtkMath from 'vtk.js/Sources/Common/Core/Math'; +import vtkBoundingBox from 'vtk.js/Sources/Common/DataModel/BoundingBox'; +import vtkIncrementalOctreeNode from 'vtk.js/Sources/Common/DataModel/IncrementalOctreeNode'; +import vtkAbstractPointLocator from 'vtk.js/Sources/Common/DataModel/AbstractPointLocator'; +import { VtkDataTypes } from 'vtk.js/Sources/Common/Core/DataArray/Constants'; + +const { vtkErrorMacro } = macro; + +function vtkIncrementalOctreePointLocator(publicAPI, model) { + // Set our className + model.classHierarchy.push('vtkIncrementalOctreePointLocator'); + + function getLeafContainer(node, pnt) { + return node.isLeaf() + ? node + : getLeafContainer(node.getChild(node.getChildIndex(pnt)), pnt); + } + + //------------------------------------------------------------------------------ + publicAPI.freeSearchStructure = () => { + model.octreeRootNode = null; + model.numberOfNodes = 0; + model.locatorPoints = null; + }; + + //------------------------------------------------------------------------------ + publicAPI.findClosestPointInLeafNode = (leafNode, point) => { + // NOTE: dist2 MUST be initiated with a very huge value below, but instead of + // model.octreeMaxDimSize * model.octreeMaxDimSize * 4.0, because the point + // under check may be outside the octree and hence the squared distance can + // be greater than the latter or other similar octree-based specific values. + + // TODO: Check + // let dist2 = VTK_DOUBLE_MAX; + let dist2 = Number.MAX_VALUE; + + if (leafNode.getPointIdSet() == null) { + return [-1, dist2]; + } + + let numPts = 0; + let tmpDst = 0.0; + const tmpPnt = []; + let tmpIdx = -1; + let pntIdx = -1; + let idList = leafNode.getPointIdSet(); + numPts = idList.length; + + for (let i = 0; i < numPts; i++) { + tmpIdx = idList[i]; + model.locatorPoints.getPoint(tmpIdx, tmpPnt); + tmpDst = vtkMath.distance2BetweenPoints(tmpPnt, point); + if (tmpDst < dist2) { + dist2 = tmpDst; + pntIdx = tmpIdx; + } + + if (dist2 === 0.0) { + break; + } + } + + idList = null; + + return [pntIdx, dist2]; + }; + + //------------------------------------------------------------------------------ + publicAPI.initPointInsertion = (points, bounds, estNumPts = 0) => { + let i = 0; + let bbIndex = 0; + + if (points == null) { + vtkErrorMacro('a valid vtkPoints object required for point insertion'); + return false; + } + + // destroy the existing octree, if any + publicAPI.freeSearchStructure(); + model.locatorPoints = points; + + // obtain the threshold squared distance + model.insertTolerance2 = model.tolerance * model.tolerance; + + // Fix bounds + // (1) push out a little bit if the original volume is too flat --- a slab + // (2) pull back the x, y, and z's lower bounds a little bit such that + // points are clearly "inside" the spatial region. Point p is taken as + // "inside" range r = [r1, r2] if and only if r1 < p <= r2. + model.octreeMaxDimSize = 0.0; + const tmpBbox = [...bounds]; + const dimDiff = vtkBoundingBox.getLengths(bounds); + model.octreeMaxDimSize = Math.max(...dimDiff); + + if (model.buildCubicOctree) { + // make the bounding box a cube and hence descendant octants cubes too + for (i = 0; i < 3; i++) { + if (dimDiff[i] !== model.octreeMaxDimSize) { + const delta = model.octreeMaxDimSize - dimDiff[i]; + tmpBbox[2 * i] -= 0.5 * delta; + tmpBbox[2 * i + 1] += 0.5 * delta; + dimDiff[i] = model.octreeMaxDimSize; + } + } + } + + model.fudgeFactor = model.octreeMaxDimSize * 10e-6; + const minSideSize = model.octreeMaxDimSize * 10e-2; + + for (i = 0; i < 3; i++) { + if (dimDiff[i] < minSideSize) { + // case (1) above + bbIndex = 2 * i; + const tempVal = tmpBbox[bbIndex]; + tmpBbox[bbIndex] = tmpBbox[bbIndex + 1] - minSideSize; + tmpBbox[bbIndex + 1] = tempVal + minSideSize; + } else { + // case (2) above + tmpBbox[2 * i] -= model.fudgeFactor; + } + } + + // init the octree with an empty leaf node + model.octreeRootNode = vtkIncrementalOctreeNode.newInstance(); + ++model.numberOfNodes; + + // this call internally inits the middle (center) and data range, too + model.octreeRootNode.setBounds(...tmpBbox); + + return true; + }; + + //------------------------------------------------------------------------------ + publicAPI.findDuplicateFloatTypePointInVisitedLeafNode = ( + leafNode, + point + ) => { + // TODO: Check + let tmpPnt; + let tmpIdx = -1; + let pntIdx = -1; + + // float thePnt[3]; // TODO + // thePnt[0] = static_cast(point[0]); + // thePnt[1] = static_cast(point[1]); + // thePnt[2] = static_cast(point[2]); + + const idList = leafNode.getPointIdSet(); + // float* pFloat = (static_cast(model.locatorPoints.getData())).getPointer(0); + const values = model.locatorPoints.getData(); + + for (let i = 0; i < idList.length; i++) { + tmpIdx = idList[i]; + // eslint-disable-next-line no-bitwise + tmpPnt = (tmpIdx << 1) + tmpIdx; + + if ( + point[0] === values[tmpPnt] && + point[1] === values[tmpPnt + 1] && + point[2] === values[tmpPnt + 2] + ) { + pntIdx = tmpIdx; + break; + } + } + + return pntIdx; + }; + + //------------------------------------------------------------------------------ + publicAPI.findDuplicateDoubleTypePointInVisitedLeafNode = ( + leafNode, + point + ) => { + // TODO: Check + let tmpPnt; + let tmpIdx = -1; + let pntIdx = -1; + const idList = leafNode.getPointIdSet(); + + // TODO: Check + // double* pArray = (static_cast(model.locatorPoints.getData())).getPointer(0); + const values = model.locatorPoints.getData(); + + for (let i = 0; i < idList.length; i++) { + tmpIdx = idList[i]; + // eslint-disable-next-line no-bitwise + tmpPnt = (tmpIdx << 1) + tmpIdx; + + if ( + point[0] === values[tmpPnt] && + point[1] === values[tmpPnt + 1] && + point[2] === values[tmpPnt + 2] + ) { + pntIdx = tmpIdx; + break; + } + } + + return pntIdx; + }; + + //------------------------------------------------------------------------------ + publicAPI.findDuplicatePointInLeafNode = (leafNode, point) => { + if (leafNode.getPointIdSet() == null) { + return -1; + } + + // TODO: Check + return model.locatorPoints.getDataType() === VtkDataTypes.FLOAT + ? publicAPI.findDuplicateFloatTypePointInVisitedLeafNode(leafNode, point) + : publicAPI.findDuplicateDoubleTypePointInVisitedLeafNode( + leafNode, + point + ); + }; + + //------------------------------------------------------------------------------ + publicAPI.insertPoint = (ptId, x) => { + const leafcontainer = getLeafContainer(model.octreeRootNode, x); + ({ numberOfNodes: model.numberOfNodes } = leafcontainer.insertPoint( + model.locatorPoints, + x, + model.maxPointsPerLeaf, + ptId, + 1, + model.numberOfNodes + )); + }; + + //------------------------------------------------------------------------------ + publicAPI.insertUniquePoint = (point) => { + // TODO: We have a mix of let and const here. + // eslint-disable-next-line prefer-const + let { pointIdx, leafContainer } = publicAPI.isInsertedPoint(point); + if (pointIdx > -1) { + return { success: false, idx: pointIdx }; + } + // TODO: pointIdx + let numberOfNodes; + // eslint-disable-next-line prefer-const + ({ numberOfNodes, pointIdx } = leafContainer.insertPoint( + model.locatorPoints, + point, + model.maxPointsPerLeaf, + pointIdx, + 2, + model.numberOfNodes + )); + model.numberOfNodes = numberOfNodes; + return { success: true, idx: pointIdx }; + }; + + //------------------------------------------------------------------------------ + publicAPI.insertNextPoint = (x) => { + const leafContainer = getLeafContainer(model.octreeRootNode, x); + const { numberOfNodes, pointIdx } = leafContainer.insertPoint( + model.locatorPoints, + x, + model.maxPointsPerLeaf, + -1, + 2, + model.numberOfNodes + ); + model.numberOfNodes = numberOfNodes; + return pointIdx; + }; + + //------------------------------------------------------------------------------ + publicAPI.isInsertedPointForZeroTolerance = (x) => { + // the target leaf node always exists there since the root node of the + // octree has been initialized to cover all possible points to be inserted + // and therefore we do not need to check it here + const leafContainer = getLeafContainer(model.octreeRootNode, x); + const pointIdx = publicAPI.findDuplicatePointInLeafNode(leafContainer, x); + return { pointIdx, leafContainer }; + }; + + //------------------------------------------------------------------------------ + publicAPI.isInsertedPointForNonZeroTolerance = (x) => { + // minDist2 // min distance to ALL existing points + // elseDst2 // min distance to other nodes (inner boundaries) + let dist2Ext; // min distance to an EXTended set of nodes + let pntIdExt; + + // the target leaf node always exists there since the root node of the + // octree has been initialized to cover all possible points to be inserted + // and therefore we do not need to check it here + const leafContainer = getLeafContainer(model.octreeRootNode, x); + let [pointIdx, minDist2] = publicAPI.findClosestPointInLeafNode( + leafContainer, + x + ); + + if (minDist2 === 0.0) { + return { pointIdx, leafContainer }; + } + + // As no any 'duplicate' point exists in this leaf node, we need to expand + // the search scope to capture possible closer points in other nodes. + const elseDst2 = leafContainer.getDistance2ToInnerBoundary( + x, + model.octreeRootNode + ); + + if (elseDst2 < model.insertTolerance2) { + // one or multiple closer points might exist in the neighboring nodes + // TODO: dist2Ext + pntIdExt = publicAPI.findClosestPointInSphereWithTolerance( + x, + model.insertTolerance2, + leafContainer, + dist2Ext + ); + + if (dist2Ext < minDist2) { + minDist2 = dist2Ext; + pointIdx = pntIdExt; + } + } + pointIdx = minDist2 <= model.insertTolerance2 ? pointIdx : -1; + return { pointIdx, leafContainer }; + }; + + //------------------------------------------------------------------------------ + publicAPI.isInsertedPoint = (x, leafContainer) => + model.insertTolerance2 === 0.0 + ? publicAPI.isInsertedPointForZeroTolerance(x, leafContainer) + : publicAPI.isInsertedPointForNonZeroTolerance(x, leafContainer); +} + +// ---------------------------------------------------------------------------- +// Object factory +// ---------------------------------------------------------------------------- + +// ---------------------------------------------------------------------------- +// Object factory +// ---------------------------------------------------------------------------- + +function defaultValues(initialValues) { + return { + fudgeFactor: 0, + octreeMaxDimSize: 0, + buildCubicOctree: false, + maxPointsPerLeaf: 128, + insertTolerance2: 0.000001, + locatorPoints: null, + octreeRootNode: null, + numberOfNodes: 0, + ...initialValues, + }; +} + +// ---------------------------------------------------------------------------- + +export function extend(publicAPI, model, initialValues = {}) { + vtkAbstractPointLocator.extend( + publicAPI, + model, + defaultValues(initialValues) + ); + + // Make this a VTK object + macro.obj(publicAPI, model); + + macro.setGet(publicAPI, model, [ + 'fudgeFactor', + 'octreeMaxDimSize', + 'buildCubicOctree', + 'maxPointsPerLeaf', + 'insertTolerance2', + 'locatorPoints', + 'octreeRootNode', + 'numberOfNodes', + ]); + + // Object specific methods + vtkIncrementalOctreePointLocator(publicAPI, model); +} + +// ---------------------------------------------------------------------------- + +export const newInstance = macro.newInstance( + extend, + 'vtkIncrementalOctreePointLocator' +); + +// ---------------------------------------------------------------------------- + +export default { newInstance, extend }; diff --git a/Sources/Common/DataModel/IncrementalOctreePointLocator/test/testIncrementalOctreePointLocator.js b/Sources/Common/DataModel/IncrementalOctreePointLocator/test/testIncrementalOctreePointLocator.js new file mode 100644 index 00000000000..753ec40081c --- /dev/null +++ b/Sources/Common/DataModel/IncrementalOctreePointLocator/test/testIncrementalOctreePointLocator.js @@ -0,0 +1,12 @@ +import test from 'tape-catch'; +import vtkIncrementalOctreePointLocator from 'vtk.js/Sources/Common/DataModel/IncrementalOctreePointLocator'; + +test('Test vtkIncrementalOctreePointLocator instance', (t) => { + t.ok( + vtkIncrementalOctreePointLocator, + 'Make sure the class definition exists' + ); + const instance = vtkIncrementalOctreePointLocator.newInstance(); + t.ok(instance); + t.end(); +}); diff --git a/Sources/Common/DataModel/Locator/index.d.ts b/Sources/Common/DataModel/Locator/index.d.ts new file mode 100644 index 00000000000..b362da504ce --- /dev/null +++ b/Sources/Common/DataModel/Locator/index.d.ts @@ -0,0 +1,43 @@ +import { vtkObject } from "../../../interfaces"; + +/** + * + */ +export interface ILocatorInitialValues { + dataSet?: number[]; + maxLevel?: number; + level?: number; + automatic?: boolean; + tolerance?: number; + useExistingSearchStructure?: boolean; +} + +export interface vtkLocator extends vtkObject {} + +// ---------------------------------------------------------------------------- +// Static API +// ---------------------------------------------------------------------------- + +/** + * Method use to decorate a given object (publicAPI+model) with vtkLocator characteristics. + * + * @param publicAPI object on which methods will be bounds (public) + * @param model object on which data structure will be bounds (protected) + * @param {ILocatorInitialValues} [initialValues] (default: {}) + */ +export function extend( + publicAPI: object, + model: object, + initialValues?: ILocatorInitialValues +): void; + +// ---------------------------------------------------------------------------- + +/** + * vtkLocator + */ +export declare const vtkLocator: { + extend: typeof extend; +}; + +export default vtkLocator; diff --git a/Sources/Common/DataModel/Locator/index.js b/Sources/Common/DataModel/Locator/index.js new file mode 100644 index 00000000000..2ea28fcae07 --- /dev/null +++ b/Sources/Common/DataModel/Locator/index.js @@ -0,0 +1,45 @@ +import macro from 'vtk.js/Sources/macros'; + +function vtkLocator(publicAPI, model) { + // Set our className + model.classHierarchy.push('vtkLocator'); +} + +// ---------------------------------------------------------------------------- +// Object factory +// ---------------------------------------------------------------------------- + +const DEFAULT_VALUES = { + dataSet: null, + maxLevel: 8, // TODO: clamp 0, Number.MAX_VALUE + level: 8, + automatic: false, + tolerance: 0.0, // TODO: clamp 0.0, Number.MAX_VALUE + useExistingSearchStructure: false, +}; + +// ---------------------------------------------------------------------------- + +export function extend(publicAPI, model, initialValues = {}) { + Object.assign(model, DEFAULT_VALUES, initialValues); + + // Make this a VTK object + macro.obj(publicAPI, model); + + macro.get(publicAPI, model, ['level']); + + macro.setGet(publicAPI, model, [ + 'dataSet', + 'maxLevel', + 'automatic', + 'tolerance', + 'useExistingSearchStructure', + ]); + + // Object specific methods + vtkLocator(publicAPI, model); +} + +// ---------------------------------------------------------------------------- + +export default { extend }; diff --git a/Sources/Common/DataModel/Locator/test/testLocator.js b/Sources/Common/DataModel/Locator/test/testLocator.js new file mode 100644 index 00000000000..4c09a594988 --- /dev/null +++ b/Sources/Common/DataModel/Locator/test/testLocator.js @@ -0,0 +1,8 @@ +import test from 'tape-catch'; +import vtkLocator from 'vtk.js/Sources/Common/DataModel/Locator'; + +test('Test vtkLocator instance', (t) => { + t.ok(vtkLocator, 'Make sure the class definition exists'); + t.ok(vtkLocator.newInstance === undefined, 'Make sure class is abstract'); + t.end(); +}); diff --git a/Sources/Common/DataModel/Polygon/index.d.ts b/Sources/Common/DataModel/Polygon/index.d.ts index d47cd0b794d..0f5a37cb6cb 100644 --- a/Sources/Common/DataModel/Polygon/index.d.ts +++ b/Sources/Common/DataModel/Polygon/index.d.ts @@ -1,6 +1,5 @@ -import { vtkObject } from "../../../interfaces" ; -import { Vector3 } from "../../../types"; - +import { vtkObject } from "../../../interfaces"; +import { Bounds, Vector3 } from "../../../types"; export interface IPolygonInitialValues { firstPoint?: Vector3, @@ -8,6 +7,17 @@ export interface IPolygonInitialValues { tris?: Vector3[], } +/** + * Different states which pointInPolygon could return. + */ +export enum POINT_IN_POLYGON { + FAILURE, + OUTSIDE, + INSIDE, + INTERSECTION, + ON_LINE, +} + export interface vtkPolygon extends vtkObject { /** @@ -28,6 +38,25 @@ export interface vtkPolygon extends vtkObject { * defines one triangle. */ triangulate(): void; + + /** + * 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 + * inspired by Dan Sunday's algorithm found in the book Practical Geometry + * Algorithms. + * @param {Vector3} point Point to check + * @param {Vector3[]} vertices Vertices of the polygon + * @param {Bounds} bounds Bounds of the vertices + * @param {Vector3} normal normal vector of the polygon + */ + pointInPolygon( + point: Vector3, + vertices: Vector3[], + bounds: Bounds, + normal: Vector3 + ): number; } /** diff --git a/Sources/Common/DataModel/Polygon/index.js b/Sources/Common/DataModel/Polygon/index.js index 002dcf41bc7..bd3814f33b9 100644 --- a/Sources/Common/DataModel/Polygon/index.js +++ b/Sources/Common/DataModel/Polygon/index.js @@ -3,6 +3,7 @@ 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'; // ---------------------------------------------------------------------------- @@ -10,6 +11,227 @@ import { IntersectionState } from 'vtk.js/Sources/Common/DataModel/Line/Constant // ---------------------------------------------------------------------------- const EPSILON = 1e-6; +const VTK_POLYGON_TOL = 1e-8; + +const POINT_IN_POLYGON = { + FAILURE: -1, + OUTSIDE: 0, + INSIDE: 1, + INTERSECTION: 2, + ON_LINE: 3, +}; +const FLT_EPSILON = 1.1920929e-7; + +// ---------------------------------------------------------------------------- +// Global methods +// ---------------------------------------------------------------------------- + +// Given the line (p0,p1), determine if the given point is located to the left +// of, on, or to the right of a line (with the function returning >0, ==0, or +// <0 respectively). The points are assumed 3D points, but projected into +// one of x-y-z planes; hence the indices axis0 and axis1 specify which plane +// the computation is to be performed on. +function pointLocation(axis0, axis1, p0, p1, point) { + return ( + (p1[axis0] - p0[axis0]) * (point[axis1] - p0[axis1]) - + (point[axis0] - p0[axis0]) * (p1[axis1] - p0[axis1]) + ); +} + +//------------------------------------------------------------------------------ + +function pointInPolygon(point, vertices, bounds, normal) { + // Do a quick bounds check to throw out trivial cases. + // winding plane. + if ( + point[0] < bounds[0] || + point[0] > bounds[1] || + point[1] < bounds[2] || + point[1] > bounds[3] || + point[2] < bounds[4] || + point[2] > bounds[5] + ) { + return POINT_IN_POLYGON.OUTSIDE; + } + + // Check that the normal is non-zero. + if (vtkMath.normalize(normal) <= FLT_EPSILON) { + return POINT_IN_POLYGON.FAILURE; + } + + // Assess whether the point lies on the boundary of the polygon. Points on + // the boundary are considered inside the polygon. Need to define a small + // tolerance relative to the bounding box diagonal length of the polygon. + let tol2 = + VTK_POLYGON_TOL * + ((bounds[1] - bounds[0]) * (bounds[1] - bounds[0]) + + (bounds[3] - bounds[2]) * (bounds[3] - bounds[2]) + + (bounds[5] - bounds[4]) * (bounds[5] - bounds[4])); + tol2 *= tol2; + tol2 = tol2 === 0.0 ? FLT_EPSILON : tol2; + const p0 = []; + const p1 = []; + + for (let i = 0; i < vertices.length; ) { + // Check coincidence to polygon vertices + p0[0] = vertices[i++]; + p0[1] = vertices[i++]; + p0[2] = vertices[i++]; + if (vtkMath.distance2BetweenPoints(point, p0) <= tol2) { + return POINT_IN_POLYGON.INSIDE; + } + + // Check coincidence to polygon edges + const { distance, t } = vtkLine.distanceToLine(point, p0, p1); + if (distance <= tol2 && t > 0.0 && t < 1.0) { + return POINT_IN_POLYGON.INSIDE; + } + } + + // 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. + let axis0; + let axis1; + if (Math.abs(normal[0]) > Math.abs(normal[1])) { + if (Math.abs(normal[0]) > Math.abs(normal[2])) { + axis0 = 1; + axis1 = 2; + } else { + axis0 = 0; + axis1 = 1; + } + } else if (Math.abs(normal[1]) > Math.abs(normal[2])) { + axis0 = 0; + axis1 = 2; + } else { + axis0 = 0; + axis1 = 1; + } + + // Compute the winding number wn. If after processing all polygon edges + // wn == 0, then the point is outside. Otherwise, the point is inside the + // polygon. Process all polygon edges determining if there are ascending or + // descending crossings of the line axis1=constant. + let wn = 0; + for (let i = 0; i < vertices.length; ) { + p0[0] = vertices[i++]; + p0[1] = vertices[i++]; + p0[2] = vertices[i++]; + 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]; + } + + if (p0[axis1] <= point[axis1]) { + if (p1[axis1] > point[axis1]) { + // if an upward crossing + if (pointLocation(axis0, axis1, p0, p1, point) > 0) { + // if x left of edge + ++wn; // a valid up intersect, increment the winding number + } + } + } else if (p1[axis1] <= point[axis1]) { + // if a downward crossing + if (pointLocation(axis0, axis1, p0, p1, point) < 0) { + // if x right of edge + --wn; // a valid down intersect, decrement the winding number + } + } + } // Over all polygon edges + + // A winding number == 0 is outside the polygon + return wn === 0 ? POINT_IN_POLYGON.OUTSIDE : POINT_IN_POLYGON.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 {Number[]} poly + * @param {vtkPoints} points + * @param {Bound} bounds + */ +export function getBounds(poly, points, bounds) { + const n = poly.length; + const p = []; + + points.getPoint(poly[0], p); + bounds[0] = p[0]; + bounds[1] = p[0]; + bounds[2] = p[1]; + bounds[3] = p[1]; + bounds[4] = p[2]; + bounds[5] = p[2]; + + for (let j = 1; j < n; j++) { + points.getPoint(poly[j], p); + vtkBoundingBox.addPoint(bounds, p); + } + const length = vtkBoundingBox.getLength(bounds); + return vtkMath.dot(length, length); +} + +// --------------------------------------------------- +/** + * Compute the normal of a polygon. + * + * 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 {Vector3[]} 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; + + const p0 = []; + let p1 = []; + let p2 = []; + const v1 = []; + const v2 = []; + + points.getPoint(poly[0], p0); + points.getPoint(poly[1], p1); + + for (let j = 2; j < poly.length; j++) { + points.getPoint(poly[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); + + [p1, p2] = [p2, p1]; + } + + return vtkMath.normalize(normal); +} + +// ---------------------------------------------------------------------------- +// Static API +// ---------------------------------------------------------------------------- + +const STATIC = { POINT_IN_POLYGON, pointInPolygon, getBounds, getNormal }; + +// ---------------------------------------------------------------------------- +// vtkPolygon methods +// ---------------------------------------------------------------------------- function vtkPolygon(publicAPI, model) { // Set our classname @@ -239,4 +461,4 @@ export const newInstance = macro.newInstance(extend, 'vtkPolygon'); // ---------------------------------------------------------------------------- -export default { newInstance, extend }; +export default { newInstance, extend, ...STATIC }; diff --git a/Sources/Filters/General/ClipClosedSurface/Constants.js b/Sources/Filters/General/ClipClosedSurface/Constants.js new file mode 100644 index 00000000000..1ee52322173 --- /dev/null +++ b/Sources/Filters/General/ClipClosedSurface/Constants.js @@ -0,0 +1,7 @@ +export const ScalarMode = { + NONE: 0, + COLORS: 1, + LABELS: 2, +}; + +export default { ScalarMode }; diff --git a/Sources/Filters/General/ClipClosedSurface/ccsEdgeLocator.js b/Sources/Filters/General/ClipClosedSurface/ccsEdgeLocator.js new file mode 100644 index 00000000000..3802d81ab85 --- /dev/null +++ b/Sources/Filters/General/ClipClosedSurface/ccsEdgeLocator.js @@ -0,0 +1,55 @@ +export default class CCSEdgeLocator { + constructor() { + this._edgeMap = new Map(); + } + + initialize() { + this._edgeMap.clear(); + } + + insertUniqueEdge(i0, i1, i) { + // Ensure consistent ordering of edge + if (i1 < i0) { + // eslint-disable-next-line no-param-reassign + [i0, i1] = [i1, i0]; + } + + // Generate a integer key, try to make it unique + // eslint-disable-next-line no-bitwise + const key = (i1 << 16) ^ i0; + let node = this._edgeMap.get(key); + + if (!node) { + // Didn't find key, so add a new edge entry + // node = new vtkCCSEdgeLocatorNode(); + // node.ptId0 = i0; + // node.ptId1 = i1; + node = { ptId0: i0, ptId1: i1, edgeId: -1, next: null }; + this._edgeMap.set(key, node); + return [node, i]; + } + + // Search through the list for i0 and i1 + if (node.ptId0 === i0 && node.ptId1 === i1) { + return [null, node.edgeId]; + } + + while (node.next != null) { + node = node.next; + + if (node.ptId0 === i0 && node.ptId1 === i1) { + return [null, node.edgeId]; + } + } + + // No entry for i1, so make one and return + const nnode = { + ptId0: i0, + ptId1: i1, + edgeId: this._edgeMap.size - 1, + next: null, + }; + node.next = nnode; + return [nnode, i]; + } +} diff --git a/Sources/Filters/General/ClipClosedSurface/example/index.js b/Sources/Filters/General/ClipClosedSurface/example/index.js new file mode 100644 index 00000000000..4d49a85ab9e --- /dev/null +++ b/Sources/Filters/General/ClipClosedSurface/example/index.js @@ -0,0 +1,138 @@ +import 'vtk.js/Sources/favicon'; + +// Load the rendering pieces we want to use (for both WebGL and WebGPU) +import 'vtk.js/Sources/Rendering/Profiles/Geometry'; + +import vtkFullScreenRenderWindow from 'vtk.js/Sources/Rendering/Misc/FullScreenRenderWindow'; +import vtkActor from 'vtk.js/Sources/Rendering/Core/Actor'; +import vtkCamera from 'vtk.js/Sources/Rendering/Core/Camera'; +import vtkMapper from 'vtk.js/Sources/Rendering/Core/Mapper'; +import vtkPlane from 'vtk.js/Sources/Common/DataModel/Plane'; +import vtkClipClosedSurface from 'vtk.js/Sources/Filters/General/ClipClosedSurface'; +import vtkSphereSource from 'vtk.js/Sources/Filters/Sources/SphereSource'; + +// import controlPanel from './controller.html'; + +// ---------------------------------------------------------------------------- +// Standard rendering code setup +// ---------------------------------------------------------------------------- + +const fullScreenRenderer = vtkFullScreenRenderWindow.newInstance({ + background: [0, 0, 0], +}); +const renderer = fullScreenRenderer.getRenderer(); +const renderWindow = fullScreenRenderer.getRenderWindow(); + +// ---------------------------------------------------------------------------- +// Example code +// ---------------------------------------------------------------------------- + +const NAMED_COLORS = { + BANANA: [227 / 255, 207 / 255, 87 / 255], + TOMATO: [255 / 255, 99 / 255, 71 / 255], + SANDY_BROWN: [244 / 255, 164 / 255, 96 / 255], +}; + +const actor = vtkActor.newInstance(); +renderer.addActor(actor); + +const mapper = vtkMapper.newInstance({ interpolateScalarBeforeMapping: true }); +actor.setMapper(mapper); + +const cam = vtkCamera.newInstance(); +renderer.setActiveCamera(cam); +cam.setFocalPoint(0, 0, 0); +cam.setPosition(0, 0, 10); +cam.setClippingRange(0.1, 50.0); + +// Build pipeline +const source = vtkSphereSource.newInstance({ + thetaResolution: 20, + phiResolution: 11, +}); + +const bounds = source.getOutputData().getBounds(); +const center = [ + (bounds[1] + bounds[0]) / 2, + (bounds[3] + bounds[2]) / 2, + (bounds[5] + bounds[4]) / 2, +]; +const planes = []; +const plane1 = vtkPlane.newInstance({ + origin: center, + normal: [1.0, 0.0, 0.0], +}); +planes.push(plane1); +const plane2 = vtkPlane.newInstance({ + origin: center, + normal: [0.0, 1.0, 0.0], +}); +planes.push(plane2); +const plane3 = vtkPlane.newInstance({ + origin: center, + normal: [0.0, 0.0, 1.0], +}); +planes.push(plane3); + +const filter = vtkClipClosedSurface.newInstance({ + clippingPlanes: planes, + activePlaneId: 2, + clipColor: NAMED_COLORS.BANANA, + baseColor: NAMED_COLORS.TOMATO, + activePlaneColor: NAMED_COLORS.SANDY_BROWN, + passPointData: false, +}); +filter.setInputConnection(source.getOutputPort()); +filter.setScalarModeToColors(); +filter.update(); +const filterData = filter.getOutputData(); + +mapper.setInputData(filterData); + +// ---------------------------------------------------------------------------- +// UI control handling +// ---------------------------------------------------------------------------- + +// fullScreenRenderer.addController(controlPanel); + +// Warp setup +// ['scaleFactor'].forEach((propertyName) => { +// document.querySelector(`.${propertyName}`).addEventListener('input', (e) => { +// const value = Number(e.target.value); +// filter.set({ [propertyName]: value }); +// renderWindow.render(); +// }); +// }); + +// // Checkbox +// document.querySelector('.useNormal').addEventListener('change', (e) => { +// const useNormal = !!e.target.checked; +// filter.set({ useNormal }); +// renderWindow.render(); +// }); + +// // Sphere setup +// ['radius', 'thetaResolution', 'phiResolution'].forEach((propertyName) => { +// document.querySelector(`.${propertyName}`).addEventListener('input', (e) => { +// const value = Number(e.target.value); +// source.set({ [propertyName]: value }); +// renderWindow.render(); +// }); +// }); + +// ----------------------------------------------------------- + +renderer.resetCamera(); +renderWindow.render(); + +// ----------------------------------------------------------- +// Make some variables global so that you can inspect and +// modify objects in your browser's developer console: +// ----------------------------------------------------------- + +global.source = source; +global.sourceData = source.getOutputData(); +global.filter = filter; +global.filterData = filterData; +global.mapper = mapper; +global.actor = actor; diff --git a/Sources/Filters/General/ClipClosedSurface/index.d.ts b/Sources/Filters/General/ClipClosedSurface/index.d.ts new file mode 100644 index 00000000000..31f67358821 --- /dev/null +++ b/Sources/Filters/General/ClipClosedSurface/index.d.ts @@ -0,0 +1,86 @@ +import { vtkAlgorithm, vtkObject } from "../../../interfaces"; +import { Vector3 } from "../../../types"; +import vtkCellArray from "vtk.js/Sources/Common/DataModel/CellArray"; +import vtkPolyData from "vtk.js/Sources/Common/DataModel/PolyData"; +import { ScalarMode } from "./Constants"; +import vtkPlane from "../../../Common/DataModel/Plane"; + +/** + * + */ +export interface IClipClosedSurfaceInitialValues { + clippingPlanes?: vtkPlane[]; + tolerance?: number; + passPointData?: boolean; + scalarMode?: Number; + generateOutline?: boolean; + generateFaces?: boolean; + activePlaneId?: number; + baseColor?: Vector3; + clipColor?: Vector3; + activePlaneColor?: Vector3; + triangulationErrorDisplay?: boolean; +} + +type vtkClipClosedSurfaceBase = vtkObject & vtkAlgorithm; + +export interface vtkClipClosedSurface extends vtkClipClosedSurfaceBase { + /** + * + * @param {any} inData + * @param {any} outData + */ + requestData(inData: any, outData: any): void; + + /** + * Set scalarMode to NONE. + */ + setScalarModeToNone(): void; + + /** + * Set scalarMode to COLOR. + */ + setScalarModeToColor(): void; + + /** + * Set scalarMode to LABEL. + */ + setScalarModeToLabel(): void; +} + +// ---------------------------------------------------------------------------- +// Static API +// ---------------------------------------------------------------------------- + +/** + * Method use to decorate a given object (publicAPI+model) with vtkClipClosedSurface characteristics. + * + * @param publicAPI object on which methods will be bounds (public) + * @param model object on which data structure will be bounds (protected) + * @param {object} [initialValues] (default: {}) + */ +export function extend( + publicAPI: object, + model: object, + initialValues?: IClipClosedSurfaceInitialValues +): void; + +// ---------------------------------------------------------------------------- + +/** + * Method use to create a new instance of vtkClipClosedSurface + * @param {IClipClosedSurfaceInitialValues} [initialValues] for pre-setting some of its content + */ +export function newInstance( + initialValues?: IClipClosedSurfaceInitialValues +): vtkClipClosedSurface; + +/** + * vtkClipClosedSurface + */ +export declare const vtkClipClosedSurface: { + newInstance: typeof newInstance; + extend: typeof extend; +}; + +export default vtkClipClosedSurface; diff --git a/Sources/Filters/General/ClipClosedSurface/index.js b/Sources/Filters/General/ClipClosedSurface/index.js new file mode 100644 index 00000000000..f0c9a826d0a --- /dev/null +++ b/Sources/Filters/General/ClipClosedSurface/index.js @@ -0,0 +1,1161 @@ +import macro from 'vtk.js/Sources/macros'; +import * as vtkMath from 'vtk.js/Sources/Common/Core/Math'; +import vtkCellArray from 'vtk.js/Sources/Common/Core/CellArray'; +import vtkDataArray from 'vtk.js/Sources/Common/Core/DataArray'; +import { VtkDataTypes } from 'vtk.js/Sources/Common/Core/DataArray/Constants'; +import vtkPoints from 'vtk.js/Sources/Common/Core/Points'; +import vtkDataSetAttributes from 'vtk.js/Sources/Common/DataModel/DataSetAttributes'; +import vtkPolyData from 'vtk.js/Sources/Common/DataModel/PolyData'; +import vtkContourTriangulator from 'vtk.js/Sources/Filters/General/ContourTriangulator'; + +import Constants from './Constants'; +import CCSEdgeLocator from './ccsEdgeLocator'; + +const { vtkErrorMacro, capitalize } = macro; +const { ScalarMode } = Constants; + +function vtkClipClosedSurface(publicAPI, model) { + // Set our className + model.classHierarchy.push('vtkClipClosedSurface'); + + publicAPI.getMTime = () => + model.clippingPlanes.reduce( + (a, b) => (b.getMTime() > a ? b.getMTime() : a), + model.mtime + ); + + /** + * Take three colors as doubles, and convert to unsigned char. + * + * @param {Number} color1 + * @param {Number} color2 + * @param {Number} color3 + * @param {Number[3][3]} colors + */ + function createColorValues(color1, color2, color3, colors) { + const dcolors = [color1, color2, color3]; + const clamp = (n, min, max) => Math.min(Math.max(n, min), max); + + for (let i = 0; i < 3; i++) { + for (let j = 0; j < 3; j++) { + colors[i][j] = Math.round(clamp(dcolors[i][j], 0, 1) * 255); + } + } + } + + /** + * Point interpolation for clipping and contouring, given the scalar + * values (v0, v1) for the two endpoints (p0, p1). The use of this + * function guarantees perfect consistency in the results. + * + * @param {vtkPoints} points + * @param {vtkDataArray} pointData + * @param {CCSEdgeLocator} locator + * @param {Number} tol + * @param {Number} i0 + * @param {Number} i1 + * @param {Number} v0 + * @param {Number} v1 + * @param {Number} i + * @returns {Number} + */ + function interpolateEdge(points, pointData, locator, tol, i0, i1, v0, v1, i) { + // This swap guarantees that exactly the same point is computed + // for both line directions, as long as the endpoints are the same. + if (v1 > 0) { + // eslint-disable-next-line no-param-reassign + [i0, i1] = [i1, i0]; + // eslint-disable-next-line no-param-reassign + [v0, v1] = [v1, v0]; + } + + // After the above swap, i0 will be kept, and i1 will be clipped + + // Check to see if this point has already been computed + // TODO: Check + let node; + // eslint-disable-next-line no-param-reassign + [node, i] = locator.insertUniqueEdge(i0, i1, i); + if (!node) { + return i; + } + + // Get the edge and interpolate the new point + const p0 = []; + const p1 = []; + points.getPoint(i0, p0); + points.getPoint(i1, p1); + + const f = v0 / (v0 - v1); + const s = 1.0 - f; + const t = 1.0 - s; + + const p = [ + s * p0[0] + t * p1[0], + s * p0[1] + t * p1[1], + s * p0[2] + t * p1[2], + ]; + + const tol2 = tol * tol; + + // Make sure that new point is far enough from kept point + if (vtkMath.distance2BetweenPoints(p, p0) < tol2) { + node.edgeId = i0; + return i0; + } + + if (vtkMath.distance2BetweenPoints(p, p1) < tol2) { + node.edgeId = i1; + return i1; + } + + // TODO + // eslint-disable-next-line no-param-reassign + i = points.insertNextTuple(p); + // i = points.getData().length; + // points.getData().push(...p); + // TODO: Check + pointData.interpolateData(pointData, i0, i1, i, t); + + // vtkFieldData + // Store the new index in the locator + // TODO: Check + node.edgeId = i; + + return i; + } + + /** + * Method for clipping lines and copying the scalar data. + * + * @param {vtkPoints} points + * @param {vtkDataArray} pointScalars + * @param {vtvtkDataSetAttributesk} pointData + * @param {CCSEdgeLocator} edgeLocator + * @param {vtkCellArray} inputLines + * @param {vtkCellArray} outputLines + * @param {vtkDataSetAttributes} inLineData + * @param {vtkDataSetAttributes} outLineData + */ + function clipLines( + points, + pointScalars, + pointData, + edgeLocator, + inputLines, + outputLines, + inLineData, + outLineData + ) { + let numPts; + let i0; + let i1; + let v0; + let v1; + let c0; + let c1; + const linePts = []; + + const values = inputLines.getData(); + let cellId = 0; + for (let i = 0; i < values.length; i += numPts + 1, cellId++) { + numPts = values[i]; + i1 = values[i + 1]; + v1 = pointScalars.getData()[i1]; + c1 = v1 > 0; + + for (let j = 2; j <= numPts; j++) { + i0 = i1; + v0 = v1; + c0 = c1; + + i1 = values[i + j]; + v1 = pointScalars.getData()[i1]; + c1 = v1 > 0; + + // If at least one point wasn't clipped + if (c0 || c1) { + // If only one end was clipped, interpolate new point + if (c0 ? !c1 : c1) { + linePts[c0 ? 1 : 0] = interpolateEdge( + points, + pointData, + edgeLocator, + model.tolerance, + i0, + i1, + v0, + v1, + linePts[c0 ? 1 : 0] + ); + } + + // If endpoints are different, insert the line segment + if (i0 !== i1) { + linePts[0] = i0; + linePts[1] = i1; + const newCellId = outputLines.insertNextCell(linePts); + // outLineData.copyData(inLineData, cellId, newCellId); + outLineData.passData(inLineData, cellId, newCellId); + } + } + } + } + } + + /** + * Break polylines into individual lines, copying scalar values from + * inputScalars starting at firstLineScalar. If inputScalars is zero, + * then scalars will be set to color. If scalars is zero, then no + * scalars will be generated. + * + * @param {vtkCellArray} inputLines + * @param {vtkCellArray} outputLines + * @param {vtkDataArray} inputScalars + * @param {Number} firstLineScalar + * @param {vtkDataArray} scalars + * @param {Vector3} color + */ + function breakPolylines( + inputLines, + outputLines, + inputScalars, + firstLineScalar, + scalars, + color + ) { + const cellColor = [...color]; + + let cellId = 0; + const values = inputLines.getData(); + let numPts; + for (let i = 0; i < values.length; i += numPts + 1) { + numPts = values[i]; + + if (inputScalars) { + // inputScalars.getTypedTuple(firstLineScalar + cellId++, cellColor); + inputScalars.getTuple(firstLineScalar + cellId++, cellColor); + } + + for (let j = 1; j < numPts; j++) { + // TODO: Check + // outputLines.insertNextCell(2); + // outputLines.insertCellPoint(values[i + j]); + // outputLines.insertCellPoint(values[i + j + 1]); + outputLines.insertNextCell([values[i + j], values[i + j + 1]]); + if (scalars) { + // scalars.insertNextTypedTuple(cellColor); + scalars.insertNextTuple(cellColor); + } + } + } + } + + /** + * Copy polygons and their associated scalars to a new array. + * If inputScalars is set to zero, set polyScalars to color instead. + * If polyScalars is set to zero, don't generate scalars. + * + * @param {vtkCellArray} inputPolys + * @param {vtkCellArray} outputPolys + * @param {vtkDataArray} inputScalars + * @param {Number} firstPolyScalar + * @param {vtkDataArray} polyScalars + * @param {Vector3} color + */ + function copyPolygons( + inputPolys, + outputPolys, + inputScalars, + firstPolyScalar, + polyScalars, + color + ) { + if (!inputPolys) { + return; + } + + outputPolys.deepCopy(inputPolys); + + if (polyScalars) { + const scalarValue = [...color]; + const n = outputPolys.getNumberOfCells(true); + polyScalars.insertTuple(n - 1, scalarValue); + + if (inputScalars) { + for (let i = 0; i < n; i++) { + inputScalars.getTuple(i + firstPolyScalar, scalarValue); + polyScalars.setTuple(i, scalarValue); + } + } else { + for (let i = 0; i < n; i++) { + polyScalars.setTuple(i, scalarValue); + } + } + } + } + + function breakTriangleStrips( + inputStrips, + polys, + inputScalars, + firstStripScalar, + polyScalars, + color + ) { + if (inputStrips.getNumberOfCells() === 0) { + return; + } + + const values = inputStrips.getData(); + let cellId = firstStripScalar; + let numPts; + for (let i = 0; i < values.length; i += numPts + 1, cellId++) { + numPts = values[i]; + // vtkTriangleStrip.decomposeStrip(numPts, values, polys); + let p1 = values[i + 1]; + let p2 = values[i + 2]; + + for (let j = 0; j < numPts - 2; j++) { + const p3 = values[i + j + 3]; + if (j % 2) { + polys.insertNextCell([p2, p1, p3]); + } else { + polys.insertNextCell([p1, p2, p3]); + } + p1 = p2; + p2 = p3; + } + + if (polyScalars) { + const scalarValue = [...color]; + + if (inputScalars) { + // If there are input scalars, use them instead of "color" + inputScalars.getTuple(cellId, scalarValue); + } + + const n = numPts - 3; + const m = polyScalars.getNumberOfTuples(); + if (n >= 0) { + // First insert is just to allocate space + polyScalars.insertTuple(m + n, scalarValue); + + for (let k = 0; k < n; k++) { + polyScalars.setTuple(m + k, scalarValue); + } + } + } + } + } + + /** + * Given some closed contour lines, create a triangle mesh that + * fills those lines. The input lines must be single-segment lines, + * not polylines. The input lines do not have to be in order. + * Only lines from firstLine to will be used. Specify the normal + * of the clip plane, which will be opposite the normals + * of the polys that will be produced. If outCD has scalars, then color + * scalars will be added for each poly that is created. + * + * @param {vtkPolyData} polyData + * @param {Number} firstLine + * @param {Number} numLines + * @param {vtkCellArray} outputPolys + * @param {Vector3} normal + */ + function triangulateContours( + polyData, + firstLine, + numLines, + outputPolys, + normal + ) { + // If no cut lines were generated, there's nothing to do + if (numLines <= 0) { + return; + } + + const triangulationError = !vtkContourTriangulator.triangulateContours( + polyData, + firstLine, + numLines, + outputPolys, + [-normal[0], -normal[1], -normal[2]] + ); + + if (triangulationError && model.triangulationErrorDisplay) { + vtkErrorMacro('Triangulation failed, polyData may not be watertight.'); + } + } + + /** + * Break polylines into individual lines, copying scalar values from + * inputScalars starting at firstLineScalar. If inputScalars is zero, + * then scalars will be set to color. If scalars is zero, then no + * scalars will be generated. + * + * @param {Number[]} polygon + * @param {vtkPoints} points + * @param {vtkCellArray} triangles + * @returns {Boolean} + */ + // function triangulatePolygon(polygon, points, triangles) { + // return vtkContourTriangulator.triangulatePolygon( + // polygon, + // points, + // triangles + // ); + // } + + /** + * Clip and contour polys in one step, in order to guarantee + * that the contour lines exactly match the new free edges of + * the clipped polygons. This exact correspondence is necessary + * in order to guarantee that the surface remains closed. + * + * @param {vtkPoints} points + * @param {vtkDataArray} pointScalars + * @param {vtkDataSetAttributes} pointData + * @param {CCSEdgeLocator} edgeLocator + * @param {Number} triangulate + * @param {vtkCellArray} inputPolys + * @param {vtkCellArray} outputPolys + * @param {vtkCellArray} outputLines + * @param {vtkDataSetAttributes} inCellData + * @param {vtkDataSetAttributes} outPolyData + * @param {vtkDataSetAttributes} outLineData + */ + function clipAndContourPolys( + points, + pointScalars, + pointData, + edgeLocator, + triangulate, + inputPolys, + outputPolys, + outputLines, + inCellData, + outPolyData, + outLineData + ) { + const idList = model._idList; + // How many sides for output polygons? + let polyMax = Number.MAX_VALUE; + if (triangulate) { + if (triangulate < 4) { + // triangles only + polyMax = 3; + } else if (triangulate === 4) { + // allow triangles and quads + polyMax = 4; + } + } + + // eslint-disable-next-line prefer-const + let triangulationFailure = false; + + // Go through all cells and clip them + // TODO: Check + const values = inputPolys.getData(); + const linePts = []; + let cellId = 0; + let numPts; + for (let i = 0; i < values.length; i += numPts + 1, cellId++) { + numPts = values[i]; + + let i1 = values[i + numPts]; + let v1 = pointScalars.getData()[i1]; + let c1 = v1 > 0; + + // The ids for the current edge: init j0 to -1 if i1 will be clipped + let j0 = c1 ? i1 : -1; + let j1 = 0; + + // To store the ids of the contour line + linePts[0] = 0; + linePts[1] = 0; + + let idListIdx = 0; + for (let j = 1; j <= numPts; j++) { + // Save previous point info + const i0 = i1; + const v0 = v1; + const c0 = c1; + + // Generate new point info + i1 = values[i + j]; + v1 = pointScalars.getData()[i1]; + c1 = v1 > 0; + + // If at least one edge end point wasn't clipped + if (c0 || c1) { + // If only one end was clipped, interpolate new point + if (c0 ? !c1 : c1) { + j1 = interpolateEdge( + points, + pointData, + edgeLocator, + model.tolerance, + i0, + i1, + v0, + v1, + j1 + ); + + if (j1 !== j0) { + idList[idListIdx++] = j1; + j0 = j1; + } + + // Save as one end of the contour line + linePts[c0 ? 1 : 0] = j1; + } + + if (c1) { + j1 = i1; + + if (j1 !== j0) { + idList[idListIdx++] = j1; + j0 = j1; + } + } + } + } + + // Insert the clipped poly + const numPoints = idListIdx; + idList.length = numPoints; + + if (numPoints > polyMax) { + // TODO: Support triangulatePolygon + // let newCellId = outputPolys.getNumberOfCells(); + // // Triangulate the poly and insert triangles into output. + // const success = triangulatePolygon(idList, points, outputPolys); + // if (!success) { + // triangulationFailure = true; + // } + + // // Copy the attribute data to the triangle cells + // const ncells = outputPolys.getNumberOfCells(); + // for (; newCellId < ncells; newCellId++) { + // outPolyData.passData(inCellData, cellId, newCellId); + // } + + // NOTE: This code needs to be removed if the above code is supported. + // Insert the polygon without triangulating it + const newCellId = outputPolys.insertNextCell(idList); + outPolyData.passData(inCellData, cellId, newCellId); + } else if (numPoints > 2) { + // Insert the polygon without triangulating it + const newCellId = outputPolys.insertNextCell(idList); + outPolyData.passData(inCellData, cellId, newCellId); + } + + // Insert the contour line if one was created + if (linePts[0] !== linePts[1]) { + const newCellId = outputLines.insertNextCell(linePts); + outLineData.passData(inCellData, cellId, newCellId); + } + } + + if (triangulationFailure && model.triangulationErrorDisplay) { + vtkErrorMacro('Triangulation failed, output may not be watertight'); + } + } + + /** + * Squeeze the points and store them in the output. Only the points that + * are used by the cells will be saved, and the pointIds of the cells will + * be modified. + * + * @param {vtkPolyData} output + * @param {vtkPoints} points + * @param {vtkDataSetAttributes} pointData + * @param {String} outputPointDataType + */ + function squeezeOutputPoints(output, points, pointData, outputPointDataType) { + // Create a list of points used by cells + const n = points.getNumberOfPoints(); + let numNewPoints = 0; + + const outPointData = output.getPointData(); + const pointMap = []; + pointMap.length = n; + + const cellArrays = [ + output.getVerts(), + output.getLines(), + output.getPolys(), + output.getStrips(), + ]; + + // Find all the newPoints that are used by cells + cellArrays.forEach((cellArray) => { + if (!cellArray) { + return; + } + const values = cellArray.getData(); + let numPts; + let pointId; + for (let i = 0; i < values.length; i += numPts + 1) { + numPts = values[i]; + for (let j = 1; j <= numPts; j++) { + pointId = values[i + j]; + if (pointMap[pointId] === undefined) { + pointMap[pointId] = numNewPoints++; + } + } + } + }); + + // Create exactly the number of points that are required + const newPoints = vtkPoints.newInstance({ + size: numNewPoints * 3, + dataType: outputPointDataType, + }); + // outPointData.copyAllocate(pointData, numNewPoints, 0); + + const p = []; + let newPointId; + for (let pointId = 0; pointId < n; pointId++) { + newPointId = pointMap[pointId]; + if (newPointId !== undefined) { + points.getPoint(pointId, p); + newPoints.setTuple(newPointId, p); + outPointData.passData(pointData, pointId, newPointId); + // outPointData.copyData(pointData, pointId, newPointId); + } + } + + // Change the cell pointIds to reflect the new point array + // TODO: Check + cellArrays.forEach((cellArray) => { + if (!cellArray) { + return; + } + const values = cellArray.getData(); + let numPts; + let pointId; + for (let i = 0; i < values.length; i += numPts + 1) { + numPts = values[i]; + for (let j = 1; j <= numPts; j++) { + pointId = values[i + j]; + values[i + j] = pointMap[pointId]; + } + } + }); + + output.setPoints(newPoints); + } + + publicAPI.requestData = (inData, outData) => { + // implement requestData + const input = inData[0]; + const output = vtkPolyData.newInstance(); + outData[0] = output; + + if (!input) { + vtkErrorMacro('Invalid or missing input'); + return; + } + + if (model._idList == null) { + model._idList = []; + } else { + model._idList.length = 0; + } + + // Get the input points + const inputPoints = input.getPoints(); + let numPts = 0; + let inputPointsType = VtkDataTypes.FLOAT; + if (inputPoints) { + numPts = inputPoints.getNumberOfPoints(); + inputPointsType = inputPoints.getDataType(); + } + + // Force points to double precision, copy the point attributes + const points = vtkPoints.newInstance({ + size: numPts * 3, + dataType: VtkDataTypes.DOUBLE, + }); + + const pointData = vtkDataSetAttributes.newInstance(); + let inPointData = null; + + if (model.passPointData) { + inPointData = input.getPointData(); + // pointData.interpolateAllocate(inPointData, numPts, 0); + } + + const point = []; + for (let ptId = 0; ptId < numPts; ptId++) { + inputPoints.getPoint(ptId, point); + points.setTuple(ptId, point); + if (inPointData) { + // pointData.copyData(inPointData, ptId, ptId); + pointData.passData(inPointData, ptId, ptId); + } + } + + // An edge locator to avoid point duplication while clipping + const edgeLocator = new CCSEdgeLocator(); + + // A temporary polydata for the contour lines that are triangulated + const tmpContourData = vtkPolyData.newInstance(); + + // The cell scalars + let lineScalars; + let polyScalars; + let inputScalars; + + // For input scalars: the offsets to the various cell types + let firstLineScalar = 0; + let firstPolyScalar = 0; + let firstStripScalar = 0; + + // Make the colors to be used on the data + let numberOfScalarComponents = 1; + const colors = [ + [0, 0, 0], + [0, 0, 0], + [0, 0, 0], + ]; + + if (model.scalarMode === ScalarMode.COLORS) { + numberOfScalarComponents = 3; + createColorValues( + model.baseColor, + model.clipColor, + model.activePlaneColor, + colors + ); + } else if (model.scalarMode === ScalarMode.LABELS) { + colors[0][0] = 0; + colors[1][0] = 1; + colors[2][0] = 2; + } + + // This is set if we have to work with scalars. The input scalars + // will be copied if they are unsigned char with 3 components, otherwise + // new scalars will be generated. + const numVerts = input.getVerts()?.getNumberOfCells(true) || 0; + const inputLines = input.getLines(); + const numLines = inputLines?.getNumberOfCells(true) || 0; + const inputPolys = input.getPolys(); + const numPolys = inputPolys?.getNumberOfCells(true) || 0; + const numStrips = input.getStrips()?.getNumberOfCells(true) || 0; + + if (model.scalarMode !== ScalarMode.NONE) { + lineScalars = vtkDataArray.newInstance({ + dataType: VtkDataTypes.UNSIGNED_CHAR, + empty: true, + // size: 0, + // values: new Uint8Array(numLines * 3), + numberOfComponents: numberOfScalarComponents, + }); + + const tryInputScalars = input.getCellData().getScalars(); + // Get input scalars if they are RGB color scalars + if ( + tryInputScalars && + tryInputScalars.getDataType() === VtkDataTypes.UNSIGNED_CHAR && + numberOfScalarComponents === 3 && + tryInputScalars.getNumberOfComponents() === 3 + ) { + inputScalars = input.getCellData().getScalars(); + firstLineScalar = numVerts; + firstPolyScalar = numVerts + numLines; + firstStripScalar = numVerts + numLines + numPolys; + } + } + + // Break the input lines into segments, generate scalars for lines + let lines; + if (numLines > 0) { + lines = vtkCellArray.newInstance({ + dataType: inputLines.getDataType(), + values: new Uint8Array(numLines * 3), // we will have at least that amount of lines + size: 0, + }); + breakPolylines( + inputLines, + lines, + inputScalars, + firstLineScalar, + lineScalars, + colors[0] + ); + } else { + lines = vtkCellArray.newInstance({ + empty: true, + }); + } + + let polys = null; + let polyMax = 3; + if (numPolys > 0 || numStrips > 0) { + // If there are line scalars, then poly scalars are needed too + if (lineScalars) { + polyScalars = vtkDataArray.newInstance({ + dataType: VtkDataTypes.UNSIGNED_CHAR, + empty: true, + // size: 0, + // values: new Uint8Array(inputPolys.getNumberOfCells(false) * 3), + numberOfComponents: numberOfScalarComponents, + }); + } + + polys = vtkCellArray.newInstance(); + copyPolygons( + inputPolys, + polys, + inputScalars, + firstPolyScalar, + polyScalars, + colors[0] + ); + // TODO: Support triangle strips + breakTriangleStrips( + input.getStrips(), + polys, + inputScalars, + firstStripScalar, + polyScalars, + colors[0] + ); + + // Check if the input has polys and quads or just triangles + // TODO: Improve performance + polyMax = Math.max(...inputPolys.getCellSizes()); + } + + // Arrays for storing the clipped lines and polys + let newLines = vtkCellArray.newInstance({ + dataType: lines.getDataType(), + empty: true, + }); + let newPolys = null; + if (polys) { + newPolys = vtkCellArray.newInstance({ + dataType: polys.getDataType(), + empty: true, + }); + } + + // The line scalars, for coloring the outline + let inLineData = vtkDataSetAttributes.newInstance(); + inLineData.copyScalarsOn(); + inLineData.setScalars(lineScalars); + + // The poly scalars, for coloring the faces + let inPolyData = vtkDataSetAttributes.newInstance(); + inPolyData.copyScalarsOn(); + inPolyData.setScalars(polyScalars); + + // Also create output attribute data + let outLineData = vtkDataSetAttributes.newInstance(); + outLineData.copyScalarsOn(); + + let outPolyData = vtkDataSetAttributes.newInstance(); + outPolyData.copyScalarsOn(); + + const planes = model.clippingPlanes; + + // Go through the clipping planes and clip the input with each plane + for (let planeId = 0; planeId < planes.length; planeId++) { + const plane = planes[planeId]; + + let triangulate = 5; + if (planeId === planes.length - 1) { + triangulate = polyMax; + } + + const active = planeId === model.activePlaneId; + + // Convert the plane into an easy-to-evaluate function + const pc = plane.getNormal(); + // OK to modify pc because vtkPlane.getNormal() returns a copy + pc[3] = -vtkMath.dot(pc, plane.getOrigin()); + + // Create the clip scalars by evaluating the plane at each point + const numPoints = points.getNumberOfPoints(); + + // The point scalars, needed for clipping (not for the output!) + const pointScalars = vtkDataArray.newInstance({ + dataType: VtkDataTypes.DOUBLE, + size: numPoints, + }); + const pointScalarsData = pointScalars.getData(); + const pointsData = points.getData(); + let i = 0; + for (let pointId = 0; pointId < numPoints; pointId) { + pointScalarsData[pointId++] = + pointsData[i++] * pc[0] + + pointsData[i++] * pc[1] + + pointsData[i++] * pc[2] + + pc[3]; + } + + // Prepare the output scalars + // outLineData.copyAllocate(inLineData, 0, 0); + // outPolyData.copyAllocate(inPolyData, 0, 0); + + // Reset the locator + edgeLocator.initialize(); + + // Clip the lines + clipLines( + points, + pointScalars, + pointData, + edgeLocator, + lines, + newLines, + inLineData, + outLineData + ); + + // Clip the polys + if (polys) { + // Get the number of lines remaining after the clipping + const numClipLines = newLines.getNumberOfCells(true); + + // Cut the polys to generate more lines + clipAndContourPolys( + points, + pointScalars, + pointData, + edgeLocator, + triangulate, + polys, + newPolys, + newLines, + inPolyData, + outPolyData, + outLineData + ); + + // Add scalars for the newly-created contour lines + let scalars = outLineData.getScalars(); + + if (scalars) { + // Set the color to the active color if plane is active + const color = colors[1 + (active ? 1 : 0)]; + const activeColor = colors[2]; + const numNewLines = newLines.getNumberOfCells(true); + + const oldColor = []; + for (let lineId = numClipLines; lineId < numNewLines; lineId++) { + scalars.getTuple(lineId, oldColor); + if ( + numberOfScalarComponents !== 3 || + oldColor[0] !== activeColor[0] || + oldColor[1] !== activeColor[1] || + oldColor[2] !== activeColor[2] + ) { + scalars.setTuple(lineId, color); + } + } + } + + // Generate new polys from the cut lines + let cellId = newPolys.getNumberOfCells(true); + const numClipAndContourLines = newLines.getNumberOfCells(true); + + // Create a polydata for the lines + tmpContourData.setPoints(points); + tmpContourData.setLines(newLines); + tmpContourData.buildCells(); + + triangulateContours( + tmpContourData, + numClipLines, + numClipAndContourLines - numClipLines, + newPolys, + pc + ); + + // Add scalars for the newly-created polys + scalars = outPolyData.getScalars(); + + if (scalars) { + const color = colors[1 + (active ? 1 : 0)]; + + const numCells = newPolys.getNumberOfCells(true); + if (numCells > cellId) { + // The insert allocates space up to numCells - 1 + scalars.insertTuple(numCells - 1, color); + for (; cellId < numCells; cellId++) { + scalars.setTuple(cellId, color); + } + } + } + + // Add scalars to any diagnostic lines that added by + // triangulateContours(). In usual operation, no lines are added. + scalars = outLineData.getScalars(); + + if (scalars) { + const color = [0, 255, 255]; + const numCells = newLines.getNumberOfCells(true); + if (numCells > numClipAndContourLines) { + // The insert allocates space up to numCells - 1 + scalars.insertTuple(numCells - 1, color); + for ( + let lineCellId = numClipAndContourLines; + lineCellId < numCells; + lineCellId++ + ) { + scalars.setTuple(lineCellId, color); + } + } + } + } + + // Swap the lines, points, etcetera: old output becomes new input + [lines, newLines] = [newLines, lines]; + newLines.initialize(); + + if (polys) { + [polys, newPolys] = [newPolys, polys]; + newPolys.initialize(); + } + + [inLineData, outLineData] = [outLineData, inLineData]; + outLineData.initialize(); + + [inPolyData, outPolyData] = [outPolyData, inPolyData]; + outPolyData.initialize(); + } + + // Get the line scalars + const scalars = inLineData.getScalars(); + + if (model.generateOutline) { + output.setLines(lines); + } else if (scalars) { + scalars.initialize(); + } + + if (model.generateFaces) { + output.setPolys(polys); + + if (polys && scalars) { + const pScalars = inPolyData.getScalars(); + const m = scalars.getNumberOfTuples(); + const n = pScalars.getNumberOfTuples(); + + if (n > 0) { + const color = [0, 0, 0]; + + // This is just to expand the array + scalars.insertTuple(n + m - 1, color); + + // Fill in the poly scalars + for (let i = 0; i < n; i++) { + pScalars.getTuple(i, color); + scalars.setTuple(i + m, color); + } + } + } + } + + if (scalars && model.scalarMode === ScalarMode.COLORS) { + scalars.setName('Colors'); + output.getCellData().setScalars(scalars); + } else if (model.scalarMode === ScalarMode.LABELS) { + // Don't use VTK_UNSIGNED_CHAR or they will look like color scalars + // const categories = vtkSignedCharArray.newInstance(); + // categories.deepCopy(scalars); + // categories.setName("Labels"); + // output.getCellData().setScalars(categories); + // categories.delete(); + // TODO: Check + const categories = scalars.newClone(); + categories.setData(scalars.getData().slice()); + categories.setName('Labels'); + output.getCellData().setScalars(categories); + } else { + output.getCellData().setScalars(null); + } + + // Finally, store the points in the output + squeezeOutputPoints(output, points, pointData, inputPointsType); + // TODO: Check + // output.squeeze(); + outData[0] = output; + }; + + Object.keys(ScalarMode).forEach((key) => { + const name = capitalize(key.toLowerCase()); + publicAPI[`setScalarModeTo${name}`] = () => { + model.scalarMode = ScalarMode[key]; + }; + }); +} + +// ---------------------------------------------------------------------------- +// Object factory +// ---------------------------------------------------------------------------- + +const DEFAULT_VALUES = { + clippingPlanes: null, + tolerance: 1e-6, + passPointData: false, + + scalarMode: ScalarMode.NONE, + generateOutline: false, + generateFaces: true, + activePlaneId: -1, + + baseColor: null, + clipColor: null, + activePlaneColor: null, + + triangulationErrorDisplay: false, + // _idList: null, +}; + +// ---------------------------------------------------------------------------- + +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, [ + 'clippingPlanes', + 'tolerance', + 'passPointData', + 'scalarMode', + 'generateOutline', + 'generateFaces', + 'activePlaneId', + 'triangulationErrorDisplay', + ]); + + macro.setGetArray( + publicAPI, + model, + ['baseColor', 'clipColor', 'activePlaneColor'], + 3 + ); + + // Object specific methods + vtkClipClosedSurface(publicAPI, model); +} + +// ---------------------------------------------------------------------------- + +export const newInstance = macro.newInstance(extend, 'vtkClipClosedSurface'); + +// ---------------------------------------------------------------------------- + +export default { newInstance, extend, ...Constants }; diff --git a/Sources/Filters/General/ClipClosedSurface/test/testClipClosedSurface.js b/Sources/Filters/General/ClipClosedSurface/test/testClipClosedSurface.js new file mode 100644 index 00000000000..38331a28c95 --- /dev/null +++ b/Sources/Filters/General/ClipClosedSurface/test/testClipClosedSurface.js @@ -0,0 +1,66 @@ +import test from 'tape-catch'; +import vtkClipClosedSurface from 'vtk.js/Sources/Filters/General/ClipClosedSurface'; +import vtkLineSource from 'vtk.js/Sources/Filters/Sources/LineSource'; +import vtkPlane from 'vtk.js/Sources/Common/DataModel/Plane'; +import vtkMath from 'vtk.js/Sources/Common/Core/Math'; + +test('Test vtkClipClosedSurface instance', (t) => { + t.ok(vtkClipClosedSurface, 'Make sure the class definition exists'); + const instance = vtkClipClosedSurface.newInstance(); + t.ok(instance); + t.end(); +}); + +test('Test clip a vtkLineSource', (t) => { + const resolution = 10; + const point1 = [-1, 0, 0]; + const point2 = [1, 0, 0]; + const origin = [ + (point1[0] + point2[0]) / 2, + (point1[1] + point2[1]) / 2, + (point1[2] + point2[2]) / 2, + ]; + const normal = []; + vtkMath.subtract(point2, point1, normal); + vtkMath.normalize(normal); + + const line = vtkLineSource.newInstance({ + point1, + point2, + resolution, + }); + const halfLine = vtkLineSource.newInstance({ + point1: origin, + point2, + resolution: resolution / 2, + }); + const planes = []; + const plane = vtkPlane.newInstance({ + origin, + normal, + }); + planes.push(plane); + + const clipper = vtkClipClosedSurface.newInstance({ + generateOutline: true, + clippingPlanes: planes, + }); + clipper.setInputConnection(line.getOutputPort()); + clipper.update(); + const outputData = clipper.getOutputData(); + + t.equal( + outputData.getNumberOfLines(), + resolution / 2, + 'Number of lines is half the resolution' + ); + + t.ok( + vtkMath.areEquals( + outputData.getPoints().getData(), + halfLine.getOutputData().getPoints().getData() + ), + 'Compare points with halfLine' + ); + t.end(); +}); diff --git a/Sources/Filters/General/ContourTriangulator/Constants.js b/Sources/Filters/General/ContourTriangulator/Constants.js new file mode 100644 index 00000000000..db6c7bfe685 --- /dev/null +++ b/Sources/Filters/General/ContourTriangulator/Constants.js @@ -0,0 +1,3 @@ +export const CCS_POLYGON_TOLERANCE = 1e-5; + +export default { CCS_POLYGON_TOLERANCE }; diff --git a/Sources/Filters/General/ContourTriangulator/helper.js b/Sources/Filters/General/ContourTriangulator/helper.js new file mode 100644 index 00000000000..b644b3d8796 --- /dev/null +++ b/Sources/Filters/General/ContourTriangulator/helper.js @@ -0,0 +1,2130 @@ +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'; + +const { vtkErrorMacro } = macro; + +/** + * Reverse the elements between the indices firstIdx and lastIdx of the given array arr. + * + * @param {any[]} 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 {Number[]} polys + * @param {Number[]} poly + * @param {Vector3} trids + * @param {Number[]} polyEdges + * @param {Number[]} 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 {Number[]} poly + * @param {vtkPoints} points + * @param {Number[]} polyEdges + * @param {Number[]} originalEdges + * @param {Number[]} polys + * @param {Vector3} normal + */ +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 {Number[]} newPolys + * @param {Number[]} 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]) { + // TODO: Check + // polyData.getCellPoints(lineId, npts, pts); + 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 = []; + // TODO: Check + // polyData.getPointCells(endPts[endIdx], ncells, cells); + 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] + ) { + // TODO: Check + // polyData.getCellPoints(lineId, npts, pts); + 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]; + // TODO: Check + // polyData.getCellPoints(lineId, npts, pts); + 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 === 0 && + ((endIdx === 0 && poly[npoly - 2] === pts[1]) || + (endIdx === 1 && poly[1] === pts[npts - 2]))) || + (r !== 0 && + ((endIdx === 0 && poly[npoly - 2] === pts[npts - 2]) || + (endIdx === 1 && poly[1] === pts[1]))) + ) { + // matches.erase(matches.begin() + k); + 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]; + // TODO: Check + // polyData.getCellPoints(lineId, npts, pts); + 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 = pts + 1; + ptsEnd = pts + 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 {Number[][]} polys + * @param {Number[]} 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 + */ +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 {Number[][]} polys + * @param {vtkPoints} points + * @param {Number[]} polyGroups + * @param {Number[]} 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 {Number[][]} polys + * @param {vtkPoints} points + * @param {Number[]} polyEdges + * @param {Number[]} 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 {Number[]} poly + * @param {Number[]} edges + * @param {Number[]} originalEdges + */ +export function vtkCCSReversePoly(poly, edges, originalEdges) { + // TODO: Check + reverseElements(poly, 1, poly.length - 1); + edges.reverse(); + for (let i = 0; i < edges.length; i++) { + if (edges[i] >= 0) { + // TODO: Check + // vtkIdType* pts = &originalEdges[edges[i]]; + // vtkIdType npts = *pts++; + // std::reverse(pts, pts + npts); + 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 {Number[]} 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 {Number[]} outerPoly + * @param {Number[]} 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); + + if (vtkPolygon.pointInPolygon(p, n, pp, bounds, normal)) { + 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 {Number[]} 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 {Number[]} newPolys + * @param {vtkPoints} points + * @param {Number[]} polyGroups + * @param {Number[]} polyEdges + * @param {Number[]} 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 {Number[][]} polys + * @param {vtkPoints} points + * @param {Vector3} normal + * @param {Number[]} 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 {Number[]} outerPoly + * @param {Number[]} 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 {Number[]} 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 {Number[]} polys + * @param {Number[]} polyGroup + * @param {Number} outerPolyId + * @param {Number} innerPolyId + * @param {vtkPoints} points + * @param {Vector3} normal + * @param {Number[][]} 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 {Number[][]} polys + * @param {Number[]} polyEdges + * @param {Number} outerPolyId + * @param {Number} innerPolyId + * @param {vtkPoints} points + * @param {Number[][]} 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 {Number[][]} polys + * @param {vtkPoints} points + * @param {Number[][]} polyGroups + * @param {Number[]} polyEdges + * @param {Vector3} normal + */ +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; +} diff --git a/Sources/Filters/General/ContourTriangulator/index.d.ts b/Sources/Filters/General/ContourTriangulator/index.d.ts new file mode 100644 index 00000000000..24237b2fb50 --- /dev/null +++ b/Sources/Filters/General/ContourTriangulator/index.d.ts @@ -0,0 +1,117 @@ +import { vtkAlgorithm, vtkObject } from "../../../interfaces"; +import { Vector3 } from "../../../types"; +import vtkCellArray from "vtk.js/Sources/Common/DataModel/CellArray"; +import vtkPolyData from "vtk.js/Sources/Common/DataModel/PolyData"; + +/** + * + */ +export interface IContourTriangulatorInitialValues { + triangulationErrorDisplay?: boolean; +} + +type vtkContourTriangulatorBase = vtkObject & vtkAlgorithm; + +export interface vtkContourTriangulator extends vtkContourTriangulatorBase { + /** + * + * @param {any} inData + * @param {any} outData + */ + requestData(inData: any, outData: any): void; + + /** + * + */ + triangulateContours(): boolean; + + /** + * + */ + triangulatePolygon(): 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 {Vector3} normal + */ +export function triangulateContours( + polyData: vtkPolyData, + firstLine: number, + numLines: number, + polys: vtkCellArray, + normal: Vector3 +): 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 {} polygon + * @param {} points + * @param {} triangles + */ +export function triangulatePolygon( + polygon: any, + points: any, + triangles: any +): 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; diff --git a/Sources/Filters/General/ContourTriangulator/index.js b/Sources/Filters/General/ContourTriangulator/index.js new file mode 100644 index 00000000000..56171b30484 --- /dev/null +++ b/Sources/Filters/General/ContourTriangulator/index.js @@ -0,0 +1,276 @@ +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 { + vtkCCSCutHoleyPolys, + vtkCCSFindTrueEdges, + vtkCCSJoinLooseEnds, + vtkCCSMakeHoleyPolys, + vtkCCSMakePolysFromLines, + vtkCCSSplitAtPinchPoints, + vtkCCSTriangulate, +} from './helper'; + +const { vtkErrorMacro } = macro; + +//------------------------------------------------------------------------------ +function triangulateContours(polyData, firstLine, numLines, polys, normal) { + 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 != null; + vtkCCSMakePolysFromLines( + polyData, + firstLine, + firstLine + numLines, + oriented, + newPolys, + incompletePolys + ); + + // if no normal specified, then compute one from largest contour + const computedNormal = [0.0, 0.0, 1.0]; + if (normal == null) { + 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]; + } + } + normal[0] = computedNormal[0]; + normal[1] = computedNormal[1]; + normal[2] = computedNormal[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, normal); + + // 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 = []; + 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, + normal, + 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, normal)) { + triangulationFailure = true; + } + + // Some polys might be self-intersecting. Split the polys at each intersection point. + vtkCCSSplitAtPinchPoints( + newPolys, + points, + polyGroups, + polyEdges, + normal, + 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 ( + !vtkCCSTriangulate( + newPolys[polyId], + points, + polyEdges[polyId], + originalEdges, + polys, + normal + ) + ) { + triangulationFailure = false; + // Diagnostic code: show the polys as outlines + // const lines = polyData.getLines(); + // const poly = newPolys[polyId]; + // lines.insertNextCell(poly.length + 1); + // for (let jjj = 0; jjj < poly.length; jjj++) { + // lines.insertCellPoint(poly[jjj]); + // } + // lines.insertCellPoint(poly[0]); + } + } + + return !triangulationFailure; +} + +// --------------------------------------------------- +function triangulatePolygon(polygon, points, triangles) { + const n = polygon.length; + const polys = [[]]; + const poly = polys[0]; + poly.length = n; + + for (let i = 0; i < n; i++) { + poly[i] = polygon[i]; + } + + const originalEdges = []; + const polyEdges = []; + vtkCCSFindTrueEdges(polys, points, polyEdges, originalEdges); + const edges = polyEdges[0]; + + let success = true; + const normal = []; + if (vtkPolygon.getNormal(poly, points, normal) !== 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]; + 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 polys = []; + const polysArray = vtkCellArray.newInstance({ + dataType: 'Float64Array', + empty: true, + }); + output.setPolys(polysArray); + output.setPoints(input.getPoints()); + output.getPointData().passData(input.getPointData()); + + triangulationError = !triangulateContours( + input, + input.getNumberOfVerts(), + lines.getNumberOfCells(), + polys, + null + ); + + if (triangulationError && model.triangulationErrorDisplay) { + vtkErrorMacro('Triangulation failed, output might have holes.'); + } + + polysArray.setData(polys); + + return true; + }; +} + +// ---------------------------------------------------------------------------- +// Object factory +// ---------------------------------------------------------------------------- + +const DEFAULT_VALUES = { + triangulationErrorDisplay: 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, ['triangulationErrorDisplay']); + + // Object specific methods + vtkContourTriangulator(publicAPI, model); +} + +// ---------------------------------------------------------------------------- + +export const newInstance = macro.newInstance(extend, 'vtkContourTriangulator'); + +// ---------------------------------------------------------------------------- + +export default { newInstance, extend, ...STATIC }; diff --git a/Sources/Filters/General/ContourTriangulator/test/testContourTriangulator.js b/Sources/Filters/General/ContourTriangulator/test/testContourTriangulator.js new file mode 100644 index 00000000000..349e45b817f --- /dev/null +++ b/Sources/Filters/General/ContourTriangulator/test/testContourTriangulator.js @@ -0,0 +1,27 @@ +import test from 'tape-catch'; +import vtkContourTriangulator from 'vtk.js/Sources/Filters/General/ContourTriangulator'; +import { reverseElements } from '../helper'; + +test('Test vtkContourTriangulator instance', (t) => { + t.ok(vtkContourTriangulator, 'Make sure the class definition exists'); + const instance = vtkContourTriangulator.newInstance(); + t.ok(instance); + t.end(); +}); + +test('Test reverseElements', (t) => { + const originalArr = [0, 1, 2, 3, 4, 5, 6, 7, 8, 9]; + const arr1 = [...originalArr]; + const arr2 = [...originalArr]; + + reverseElements(arr1, 2, 7); + t.deepEqual(arr1, [0, 1, 7, 6, 5, 4, 3, 2, 8, 9], 'reverse elements 2 to 7'); + + reverseElements(arr1, 2, 7); + t.deepEqual(arr1, originalArr, 'reversing again gives original array'); + + reverseElements(arr1); + arr2.reverse(); + + t.deepEqual(arr1, arr2, 'without arguments does the same as .reverse'); +});