diff --git a/Examples/PBDCutting/CMakeLists.txt b/Examples/PBDCutting/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..a2aff794533c4c81eed61e23e58e45c4b97ec27f
--- /dev/null
+++ b/Examples/PBDCutting/CMakeLists.txt
@@ -0,0 +1,37 @@
+###########################################################################
+#
+# Copyright (c) Kitware, Inc.
+#
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+#
+# http://www.apache.org/licenses/LICENSE-2.0.txt
+#
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+#
+###########################################################################
+if(iMSTK_USE_OpenHaptics)
+
+ project(Example-PBDCutting)
+
+ #-----------------------------------------------------------------------------
+ # Create executable
+ #-----------------------------------------------------------------------------
+ imstk_add_executable(${PROJECT_NAME} PBDCuttingExample.cpp)
+
+ #-----------------------------------------------------------------------------
+ # Add the target to Examples folder
+ #-----------------------------------------------------------------------------
+ SET_TARGET_PROPERTIES (${PROJECT_NAME} PROPERTIES FOLDER Examples/Haptics)
+
+ #-----------------------------------------------------------------------------
+ # Link libraries to executable
+ #-----------------------------------------------------------------------------
+ target_link_libraries(${PROJECT_NAME} SimulationManager apiUtilities Filtering)
+
+endif()
diff --git a/Examples/PBDCutting/PBDCuttingExample.cpp b/Examples/PBDCutting/PBDCuttingExample.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..e20e4feddbcc7ee77b5124d3dd394a0eb3c6556c
--- /dev/null
+++ b/Examples/PBDCutting/PBDCuttingExample.cpp
@@ -0,0 +1,278 @@
+/*=========================================================================
+
+ Library: iMSTK
+
+ Copyright (c) Kitware, Inc. & Center for Modeling, Simulation,
+ & Imaging in Medicine, Rensselaer Polytechnic Institute.
+
+ Licensed under the Apache License, Version 2.0 (the "License");
+ you may not use this file except in compliance with the License.
+ You may obtain a copy of the License at
+
+ http://www.apache.org/licenses/LICENSE-2.0.txt
+
+ Unless required by applicable law or agreed to in writing, software
+ distributed under the License is distributed on an "AS IS" BASIS,
+ WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ See the License for the specific language governing permissions and
+ limitations under the License.
+
+=========================================================================*/
+
+#include "imstkLogger.h"
+#include "imstkNew.h"
+#include "imstkScene.h"
+#include "imstkSceneManager.h"
+#include "imstkVTKViewer.h"
+
+// Objects
+#include "imstkCamera.h"
+#include "imstkCollidingObject.h"
+#include "imstkLight.h"
+#include "imstkPbdObject.h"
+
+// Geometry
+#include "imstkMeshIO.h"
+#include "imstkCapsule.h"
+#include "imstkPlane.h"
+#include "imstkSurfaceMesh.h"
+
+// Devices and controllers
+#include "imstkHapticDeviceClient.h"
+#include "imstkHapticDeviceManager.h"
+#include "imstkKeyboardDeviceClient.h"
+#include "imstkKeyboardSceneControl.h"
+//#include "imstkLaparoscopicToolController.h"
+#include "imstkMouseSceneControl.h"
+#include "imstkSceneObjectController.h"
+
+// Collisions and Models
+#include "imstkCollisionGraph.h"
+#include "imstkCollisionPair.h"
+#include "imstkCollisionDetection.h"
+#include "imstkPbdModel.h"
+#include "imstkPbdObjectCuttingPair.h"
+#include "imstkRenderMaterial.h"
+#include "imstkVisualModel.h"
+
+#include "imstkSurfaceMeshCut.h"
+
+using namespace imstk;
+
+// Parameters to play with
+const double width = 50.0;
+const double height = 50.0;
+const int nRows = 6;
+const int nCols = 6;
+
+///
+/// \brief Creates cloth geometry
+///
+static std::shared_ptr<SurfaceMesh>
+makeClothGeometry(const double width,
+ const double height,
+ const int nRows,
+ const int nCols)
+{
+ imstkNew<SurfaceMesh> clothMesh("Cloth_SurfaceMesh");
+
+ imstkNew<VecDataArray<double, 3>> verticesPtr(nRows * nCols);
+ VecDataArray<double, 3>& vertices = *verticesPtr.get();
+ const double dy = width / static_cast<double>(nCols - 1);
+ const double dx = height / static_cast<double>(nRows - 1);
+ for (int i = 0; i < nRows; i++)
+ {
+ for (int j = 0; j < nCols; j++)
+ {
+ vertices[i * nCols + j] = Vec3d(dx * static_cast<double>(i), 1.0, dy * static_cast<double>(j));
+ }
+ }
+
+ // Add connectivity data
+ imstkNew<VecDataArray<int, 3>> indicesPtr;
+ VecDataArray<int, 3>& indices = *indicesPtr.get();
+ for (int i = 0; i < nRows - 1; i++)
+ {
+ for (int j = 0; j < nCols - 1; j++)
+ {
+ const int index1 = i * nCols + j;
+ const int index2 = index1 + nCols;
+ const int index3 = index1 + 1;
+ const int index4 = index2 + 1;
+
+ // Interleave [/][\] pattern
+ if (i % 2 ^ j % 2)
+ {
+ indices.push_back(Vec3i(index1, index2, index3));
+ indices.push_back(Vec3i(index4, index3, index2));
+ }
+ else
+ {
+ indices.push_back(Vec3i(index2, index4, index1));
+ indices.push_back(Vec3i(index4, index3, index1));
+ }
+ }
+ }
+
+ clothMesh->initialize(verticesPtr, indicesPtr);
+
+ return clothMesh;
+}
+
+///
+/// \brief Creates cloth object
+///
+static std::shared_ptr<PbdObject>
+makeClothObj(const std::string& name,
+ const double width,
+ const double height,
+ const int nRows,
+ const int nCols)
+{
+ imstkNew<PbdObject> clothObj(name);
+
+ // Setup the Geometry
+ std::shared_ptr<SurfaceMesh> clothMesh(makeClothGeometry(width, height, nRows, nCols));
+
+ // Setup the Parameters
+ imstkNew<PBDModelConfig> pbdParams;
+ pbdParams->enableConstraint(PbdConstraint::Type::Distance, 1.0e3);
+ pbdParams->enableConstraint(PbdConstraint::Type::Dihedral, 1.0e2);
+ pbdParams->m_fixedNodeIds = { 0, static_cast<size_t>(nCols) - 1 };
+ pbdParams->m_uniformMassValue = width * height / ((double)nRows * (double)nCols);
+ pbdParams->m_gravity = Vec3d(0.0, -9.8, 0.0);
+ pbdParams->m_defaultDt = 0.005;
+ pbdParams->m_iterations = 5;
+
+ // Setup the Model
+ imstkNew<PbdModel> pbdModel;
+ pbdModel->setModelGeometry(clothMesh);
+ pbdModel->configure(pbdParams);
+
+ // Setup the VisualModel
+ imstkNew<RenderMaterial> material;
+ material->setBackFaceCulling(false);
+ material->setDisplayMode(RenderMaterial::DisplayMode::WireframeSurface);
+
+ imstkNew<VisualModel> visualModel(clothMesh);
+ visualModel->setRenderMaterial(material);
+
+ // Setup the Object
+ clothObj->addVisualModel(visualModel);
+ clothObj->setPhysicsGeometry(clothMesh);
+ clothObj->setCollidingGeometry(clothMesh);
+ clothObj->setDynamicalModel(pbdModel);
+
+ return clothObj;
+}
+
+///
+/// \brief This example demonstrates the concept of PBD picking
+/// for haptic interaction. NOTE: Requires GeoMagic Touch device
+///
+int
+main()
+{
+ // Setup logger (write to file and stdout)
+ Logger::startLogger();
+
+ // Scene
+ imstkNew<Scene> scene("PBDCutting");
+ //scene->getConfig()->writeTaskGraph = true;
+
+ // Create a cutting plane object in the scene
+ imstkNew<Plane> planeGeom;
+ planeGeom->setWidth(40.0);
+ planeGeom->setTranslation(Vec3d(0.0, 0.0, 20.0));
+ planeGeom->setOrientationAxis(Vec3d(-1.0, 0.0, 0.0));
+ imstkNew<CollidingObject> planeObj("Plane");
+ planeObj->setVisualGeometry(planeGeom);
+ planeObj->setCollidingGeometry(planeGeom);
+ scene->addSceneObject(planeObj);
+
+ // Create a pbd cloth object in the scene
+ std::shared_ptr<PbdObject> clothObj = makeClothObj("Cloth", width, height, nRows, nCols);
+ scene->addSceneObject(clothObj);
+
+ // Add interaction pair for pbd cutting
+
+ // Device Server
+ imstkNew<HapticDeviceManager> server;
+ std::shared_ptr<HapticDeviceClient> client = server->makeDeviceClient();
+
+ // Create the virtual coupling object controller
+ imstkNew<SceneObjectController> controller(planeObj, client);
+ scene->addController(controller);
+
+ // Camera
+ scene->getActiveCamera()->setPosition(Vec3d(1.0, 1.0, 1.0) * 100.0);
+ scene->getActiveCamera()->setFocalPoint(Vec3d(0, -50, 0));
+
+ // Light
+ imstkNew<DirectionalLight> light("light");
+ light->setFocalPoint(Vec3d(5.0, -8.0, -5.0));
+ light->setIntensity(1.0);
+ scene->addLight(light);
+
+ //Run the simulation
+ {
+ // Setup a viewer to render in its own thread
+ imstkNew<VTKViewer> viewer("Viewer");
+ viewer->setActiveScene(scene);
+
+ // Setup a scene manager to advance the scene in its own thread
+ imstkNew<SceneManager> sceneManager("Scene Manager");
+ sceneManager->setActiveScene(scene);
+ viewer->addChildThread(sceneManager); // SceneManager will start/stop with viewer
+
+ // Add server of haptic device to viewer
+ viewer->addChildThread(server);
+
+ // Add mouse and keyboard controls to the viewer
+ {
+ imstkNew<MouseSceneControl> mouseControl(viewer->getMouseDevice());
+ mouseControl->setSceneManager(sceneManager);
+ viewer->addControl(mouseControl);
+
+ imstkNew<KeyboardSceneControl> keyControl(viewer->getKeyboardDevice());
+ keyControl->setSceneManager(sceneManager);
+ keyControl->setViewer(viewer);
+ viewer->addControl(keyControl);
+ }
+
+ // Queue keypress to be called after scene thread
+ queueConnect<KeyPressEvent>(viewer->getKeyboardDevice(), EventType::KeyEvent, sceneManager, [&](KeyPressEvent* e)
+ {
+ // When i is pressed replace the PBD cloth with a cut one
+ if (e->m_key == 'i' && e->m_keyPressType == KEY_PRESS)
+ {
+ // This has a number of issues that make it not physically realistic
+ // - Mass is not conservative when interpolated from subdivide
+ // - Constraint resting lengths are not correctly reinited
+ std::shared_ptr<SurfaceMesh> clothMesh = std::dynamic_pointer_cast<SurfaceMesh>(clothObj->getPhysicsGeometry());
+ imstkNew<SurfaceMeshCut> surfCut;
+ surfCut->setInputMesh(clothMesh);
+ surfCut->setPlane(planeGeom);
+ surfCut->update();
+ std::shared_ptr<SurfaceMesh> newClothMesh = surfCut->getOutputMesh();
+
+ // RenderDelegates cannot visually have entire geometries swapped yet, so even though we could just set the geometry
+ // on the model, you would not visually see it. Instead we replace the vertex and index buffers of the existing one
+ // Another issue here is that initial geometry is not remapped so reset will not reset to undeformed config
+ clothMesh->setInitialVertexPositions(std::make_shared<VecDataArray<double, 3>>(*newClothMesh->getVertexPositions()));
+ clothMesh->setVertexPositions(newClothMesh->getVertexPositions());
+ clothMesh->setTriangleIndices(newClothMesh->getTriangleIndices());
+ clothMesh->setVertexAttribute("Velocities", newClothMesh->getVertexAttribute("Velocities"));
+ clothMesh->modified();
+
+ // Re-setup the constraints on the object
+ clothObj->initialize();
+ }
+ });
+
+ // Start viewer running, scene as paused
+ sceneManager->requestStatus(ThreadStatus::Paused);
+ viewer->start();
+ }
+ return 0;
+}
diff --git a/Source/Filtering/imstkSurfaceMeshCut.cpp b/Source/Filtering/imstkSurfaceMeshCut.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..3c7f0f11ff7b2a6b37dec2067ecc980ee4f26981
--- /dev/null
+++ b/Source/Filtering/imstkSurfaceMeshCut.cpp
@@ -0,0 +1,469 @@
+/*=========================================================================
+
+ Library: iMSTK
+
+ Copyright (c) Kitware, Inc. & Center for Modeling, Simulation,
+ & Imaging in Medicine, Rensselaer Polytechnic Institute.
+
+ Licensed under the Apache License, Version 2.0 (the "License");
+ you may not use this file except in compliance with the License.
+ You may obtain a copy of the License at
+
+ http://www.apache.org/licenses/LICENSE-2.0.txt
+
+ Unless required by applicable law or agreed to in writing, software
+ distributed under the License is distributed on an "AS IS" BASIS,
+ WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ See the License for the specific language governing permissions and
+ limitations under the License.
+
+=========================================================================*/
+
+#include "imstkSurfaceMeshCut.h"
+
+#include "imstkGeometryUtilities.h"
+#include "imstkLogger.h"
+#include "imstkPlane.h"
+
+#include <map>
+#include <set>
+
+namespace imstk
+{
+SurfaceMeshCut::SurfaceMeshCut()
+{
+ setNumberOfInputPorts(1);
+ setNumberOfOutputPorts(1);
+ setOutput(std::make_shared<SurfaceMesh>());
+}
+
+std::shared_ptr<SurfaceMesh>
+SurfaceMeshCut::getOutputMesh()
+{
+ return std::dynamic_pointer_cast<SurfaceMesh>(getOutput(0));
+}
+
+void
+SurfaceMeshCut::setInputMesh(std::shared_ptr<SurfaceMesh> inputMesh)
+{
+ setInput(inputMesh, 0);
+}
+
+void
+SurfaceMeshCut::requestUpdate()
+{
+ // input and output SurfaceMesh
+ std::shared_ptr<SurfaceMesh> inputSurf = std::dynamic_pointer_cast<SurfaceMesh>(getInput(0));
+ if (inputSurf == nullptr)
+ {
+ LOG(WARNING) << "Missing required SurfaceMesh input";
+ return;
+ }
+ std::shared_ptr<SurfaceMesh> outputSurf = std::dynamic_pointer_cast<SurfaceMesh>(getOutput(0));
+ outputSurf->deepCopy(inputSurf);
+
+ // vertices on the cutting path and whether they will be split
+ std::map<int, bool> cutVerts;
+
+ generateCutData(m_Plane, outputSurf);
+
+ // refinement
+ refinement(outputSurf, cutVerts);
+
+ // split cutting vertices
+ splitVerts(outputSurf, cutVerts);
+
+ // vtkPolyData to output SurfaceMesh
+ setOutput(outputSurf);
+}
+
+void
+SurfaceMeshCut::refinement(std::shared_ptr<SurfaceMesh> outputSurf, std::map<int, bool>& cutVerts)
+{
+ // map between one exsiting edge to the new vert generated from the cutting
+ std::map<std::pair<int, int>, int> edgeVertMap;
+
+ auto triangles = outputSurf->getTriangleIndices();
+ auto vertices = outputSurf->getVertexPositions();
+
+ for (size_t i = 0; i < m_CutData->size(); i++)
+ {
+ CutData curCutData = m_CutData->at(i);
+ auto curCutType = curCutData.cutType;
+ int triId = curCutData.triId;
+ int ptId0 = curCutData.ptIds[0];
+ int ptId1 = curCutData.ptIds[1];
+ Vec3d cood0 = curCutData.cutCoords[0];
+ Vec3d cood1 = curCutData.cutCoords[1];
+
+ if (curCutType == CutType::EDGE || curCutType == CutType::EDGE_VERT)
+ {
+ int newPtId = -1;
+ auto edge0 = std::make_pair(ptId0, ptId1);
+ auto edge1 = std::make_pair(ptId1, ptId0);
+
+ // add new point
+ if (edgeVertMap.find(edge1) != edgeVertMap.end())
+ {
+ newPtId = edgeVertMap[edge1];
+ }
+ else
+ {
+ newPtId = vertices->size();
+ vertices->push_back(cood0);
+ edgeVertMap[edge0] = newPtId;
+ }
+
+ // update triangle indices
+ Vec3i ptIds = (*triangles)[triId];
+ int ptId2 = -1;
+ if (ptIds[0] != ptId0 && ptIds[0] != ptId1)
+ {
+ ptId2 = ptIds[0];
+ }
+ else if (ptIds[1] != ptId0 && ptIds[1] != ptId1)
+ {
+ ptId2 = ptIds[1];
+ }
+ else
+ {
+ ptId2 = ptIds[2];
+ }
+ (*triangles)[triId] = Vec3i(ptId2, ptId0, newPtId);
+ triangles->push_back(Vec3i(ptId2, newPtId, ptId1));
+
+ // add vertices to cutting path
+ if (curCutType == CutType::EDGE_VERT)
+ {
+ cutVerts[ptId2] = (cutVerts.find(ptId2) != cutVerts.end()) ? true : false;
+ cutVerts[newPtId] = (cutVerts.find(newPtId) != cutVerts.end()) ? true : false;
+ }
+ }
+ else if (curCutType == CutType::EDGE_EDGE)
+ {
+ int newPtId0 = -1;
+ int newPtId1 = -1;
+
+ Vec3i ptIds = (*triangles)[triId];
+ int ptId2 = -1;
+ if (ptIds[0] != ptId0 && ptIds[0] != ptId1)
+ {
+ ptId2 = ptIds[0];
+ }
+ else if (ptIds[1] != ptId0 && ptIds[1] != ptId1)
+ {
+ ptId2 = ptIds[1];
+ }
+ else
+ {
+ ptId2 = ptIds[2];
+ }
+
+ auto edge00 = std::make_pair(ptId2, ptId0);
+ auto edge01 = std::make_pair(ptId0, ptId2);
+ auto edge10 = std::make_pair(ptId1, ptId2);
+ auto edge11 = std::make_pair(ptId2, ptId1);
+
+ // add new points
+ if (edgeVertMap.find(edge01) != edgeVertMap.end())
+ {
+ newPtId0 = edgeVertMap[edge01];
+ }
+ else
+ {
+ newPtId0 = vertices->size();
+ vertices->push_back(cood0);
+ edgeVertMap[edge00] = newPtId0;
+ }
+ if (edgeVertMap.find(edge11) != edgeVertMap.end())
+ {
+ newPtId1 = edgeVertMap[edge11];
+ }
+ else
+ {
+ newPtId1 = vertices->size();
+ vertices->push_back(cood1);
+ edgeVertMap[edge10] = newPtId1;
+ }
+
+ // update triangle indices
+ (*triangles)[triId] = Vec3i(ptId2, newPtId0, newPtId1);
+ triangles->push_back(Vec3i(newPtId0, ptId0, ptId1));
+ triangles->push_back(Vec3i(newPtId0, ptId1, newPtId1));
+
+ // add vertices to cutting path
+ cutVerts[newPtId0] = (cutVerts.find(newPtId0) != cutVerts.end()) ? true : false;
+ cutVerts[newPtId1] = (cutVerts.find(newPtId1) != cutVerts.end()) ? true : false;
+ }
+ else if (curCutType == CutType::VERT_VERT)
+ {
+ // add vertices to cutting path
+ cutVerts[ptId0] = (cutVerts.find(ptId0) != cutVerts.end()) ? true : false;
+ cutVerts[ptId1] = (cutVerts.find(ptId1) != cutVerts.end()) ? true : false;
+ }
+ else
+ {
+ //throw error
+ }
+ }
+ outputSurf->setInitialVertexPositions(std::make_shared<VecDataArray<double, 3>>(*vertices));
+ outputSurf->setVertexPositions(std::make_shared<VecDataArray<double, 3>>(*vertices));
+ outputSurf->setTriangleIndices(std::make_shared<VecDataArray<int, 3>>(*triangles));
+ outputSurf->modified();
+}
+
+void
+SurfaceMeshCut::splitVerts(std::shared_ptr<SurfaceMesh> outputSurf, std::map<int, bool>& cutVerts)
+{
+ auto triangles = outputSurf->getTriangleIndices();
+ auto vertices = outputSurf->getVertexPositions();
+
+ // build vertexToTriangleListMap
+ std::vector<std::set<int>> vertexNeighborTriangles;
+ vertexNeighborTriangles.clear();
+ vertexNeighborTriangles.resize(vertices->size());
+
+ int triangleId = 0;
+
+ for (const auto& t : *triangles)
+ {
+ vertexNeighborTriangles.at(t[0]).insert(triangleId);
+ vertexNeighborTriangles.at(t[1]).insert(triangleId);
+ vertexNeighborTriangles.at(t[2]).insert(triangleId);
+ triangleId++;
+ }
+
+ // split cutting vertices
+ for (const auto& cutVert : cutVerts)
+ {
+ if (cutVert.second == false && !vertexOnBoundary(triangles, vertexNeighborTriangles[cutVert.first]))
+ {
+ // do not split vertex since it's the cutting end in surface
+ }
+ else
+ {
+ // split vertex
+ int newPtId = vertices->size();
+ vertices->push_back((*vertices)[cutVert.first]);
+
+ for (const auto& t : vertexNeighborTriangles[cutVert.first])
+ {
+ //if triangle on the negative side
+ Vec3d pt0 = (*vertices)[(*triangles)[t][0]];
+ Vec3d pt1 = (*vertices)[(*triangles)[t][1]];
+ Vec3d pt2 = (*vertices)[(*triangles)[t][2]];
+
+ if (pointOnPlaneSide(pt0) < 0 || pointOnPlaneSide(pt1) < 0 || pointOnPlaneSide(pt2) < 0)
+ {
+ for (int i = 0; i < 3; i++)
+ {
+ if ((*triangles)[t][i] == cutVert.first)
+ {
+ (*triangles)[t][i] = newPtId;
+ }
+ }
+ }
+ }
+ }
+ }
+
+ outputSurf->setInitialVertexPositions(std::make_shared<VecDataArray<double, 3>>(*vertices));
+ outputSurf->setVertexPositions(std::make_shared<VecDataArray<double, 3>>(*vertices));
+ outputSurf->setTriangleIndices(std::make_shared<VecDataArray<int, 3>>(*triangles));
+ outputSurf->modified();
+}
+
+int
+SurfaceMeshCut::pointOnPlaneSide(Vec3d pt)
+{
+ double normalProjection = m_Plane->getNormal().dot(pt - m_Plane->getPosition());
+ if (normalProjection > m_Epsilon)
+ {
+ return 1;
+ }
+ else if (normalProjection < -m_Epsilon)
+ {
+ return -1;
+ }
+ else
+ {
+ return 0;
+ }
+ return 0;
+}
+
+bool
+SurfaceMeshCut::vertexOnBoundary(std::shared_ptr<VecDataArray<int, 3>> triangleIndices, std::set<int>& triSet)
+{
+ std::set<int> nonRepeatNeighborVerts;
+ for (const auto& t : triSet)
+ {
+ for (int i = 0; i < 3; i++)
+ {
+ int ptId = (*triangleIndices)[t][i];
+ if (nonRepeatNeighborVerts.find(ptId) != nonRepeatNeighborVerts.end())
+ {
+ nonRepeatNeighborVerts.erase(ptId);
+ }
+ else
+ {
+ nonRepeatNeighborVerts.insert(ptId);
+ }
+ }
+ }
+ if (nonRepeatNeighborVerts.size() >= 2)
+ {
+ return true;
+ }
+ else
+ {
+ return false;
+ }
+}
+
+void
+SurfaceMeshCut::generateCutData(std::shared_ptr<Plane> plane, std::shared_ptr<SurfaceMesh> inputSurf)
+{
+ auto triangles = inputSurf->getTriangleIndices();
+ auto vertices = inputSurf->getVertexPositions();
+
+ m_CutData->clear();
+ std::set<std::pair<int, int>> repeatEdges; // make sure not boundary edge in vert-vert case
+
+ for (int i = 0; i < triangles->size(); i++)
+ {
+ auto& tri = (*triangles)[i];
+ CutData newCutData;
+
+ auto ptSide = Vec3i(pointOnPlaneSide((*vertices)[tri[0]]),
+ pointOnPlaneSide((*vertices)[tri[1]]),
+ pointOnPlaneSide((*vertices)[tri[2]]));
+
+ switch (ptSide.squaredNorm())
+ {
+ case 1:
+ for (int j = 0; j < 3; j++)
+ {
+ if (ptSide[j] != 0)
+ {
+ int idx0 = (j + 1) % 3;
+ int idx1 = (j + 2) % 3;
+ int ptId0 = tri[idx0];
+ int ptId1 = tri[idx1];
+
+ std::pair<int, int> cutEdge = std::make_pair(ptId1, ptId0);
+ // if find cut triangle from the other side of the edge, add newCutData
+ if (repeatEdges.find(cutEdge) != repeatEdges.end())
+ {
+ newCutData.cutType = CutType::VERT_VERT;
+ newCutData.triId = i;
+ newCutData.ptIds[0] = ptId0;
+ newCutData.ptIds[1] = ptId1;
+ newCutData.cutCoords[0] = (*vertices)[ptId0];
+ newCutData.cutCoords[1] = (*vertices)[ptId1];
+ m_CutData->push_back(newCutData);
+ }
+ else
+ {
+ repeatEdges.insert(std::make_pair(ptId0, ptId1));
+ }
+ }
+ }
+ break;
+ case 2:
+ if (ptSide.sum() == 0)
+ {
+ for (int j = 0; j < 3; j++)
+ {
+ if (ptSide[j] == 0)
+ {
+ int idx0 = (j + 1) % 3;
+ int idx1 = (j + 2) % 3;
+ int ptId0 = tri[idx0];
+ int ptId1 = tri[idx1];
+ Vec3d pos0 = (*vertices)[ptId0];
+ Vec3d pos1 = (*vertices)[ptId1];
+ Vec3d p = m_Plane->getPosition();
+ Vec3d n = m_Plane->getNormal();
+
+ newCutData.cutType = CutType::EDGE_VERT;
+ newCutData.triId = i;
+ newCutData.ptIds[0] = ptId0;
+ newCutData.ptIds[1] = ptId1;
+ newCutData.cutCoords[0] = (p - pos0).dot(n) / (pos1 - pos0).dot(n) * (pos1 - pos0) + pos0;
+ newCutData.cutCoords[1] = (*vertices)[tri[j]];
+ m_CutData->push_back(newCutData);
+ }
+ }
+ }
+ else
+ {
+ // newCutData.cutType = CutType::VERT;
+ }
+ break;
+ case 3:
+ if (ptSide.sum() == 1)
+ {
+ for (int j = 0; j < 3; j++)
+ {
+ if (ptSide[j] < 0)
+ {
+ int idx0 = (j + 1) % 3;
+ int idx1 = (j + 2) % 3;
+ int ptId0 = tri[idx0];
+ int ptId1 = tri[idx1];
+ int ptId2 = tri[j];
+ Vec3d pos0 = (*vertices)[ptId0];
+ Vec3d pos1 = (*vertices)[ptId1];
+ Vec3d pos2 = (*vertices)[ptId2];
+ Vec3d p = m_Plane->getPosition();
+ Vec3d n = m_Plane->getNormal();
+
+ newCutData.cutType = CutType::EDGE_EDGE;
+ newCutData.triId = i;
+ newCutData.ptIds[0] = ptId0;
+ newCutData.ptIds[1] = ptId1;
+ newCutData.cutCoords[0] = (p - pos0).dot(n) / (pos2 - pos0).dot(n) * (pos2 - pos0) + pos0;
+ newCutData.cutCoords[1] = (p - pos1).dot(n) / (pos2 - pos1).dot(n) * (pos2 - pos1) + pos1;
+ m_CutData->push_back(newCutData);
+ }
+ }
+ }
+ else if (ptSide.sum() == -1)
+ {
+ for (int j = 0; j < 3; j++)
+ {
+ if (ptSide[j] > 0)
+ {
+ int idx0 = (j + 1) % 3;
+ int idx1 = (j + 2) % 3;
+ int ptId0 = tri[idx0];
+ int ptId1 = tri[idx1];
+ int ptId2 = tri[j];
+ Vec3d pos0 = (*vertices)[ptId0];
+ Vec3d pos1 = (*vertices)[ptId1];
+ Vec3d pos2 = (*vertices)[ptId2];
+ Vec3d p = m_Plane->getPosition();
+ Vec3d n = m_Plane->getNormal();
+
+ newCutData.cutType = CutType::EDGE_EDGE;
+ newCutData.triId = i;
+ newCutData.ptIds[0] = ptId0;
+ newCutData.ptIds[1] = ptId1;
+ newCutData.cutCoords[0] = (p - pos0).dot(n) / (pos2 - pos0).dot(n) * (pos2 - pos0) + pos0;
+ newCutData.cutCoords[1] = (p - pos1).dot(n) / (pos2 - pos1).dot(n) * (pos2 - pos1) + pos1;
+ m_CutData->push_back(newCutData);
+ }
+ }
+ }
+ else
+ {
+ // no intersection
+ }
+ break;
+ default:
+ break;
+ }
+ }
+}
+}
\ No newline at end of file
diff --git a/Source/Filtering/imstkSurfaceMeshCut.h b/Source/Filtering/imstkSurfaceMeshCut.h
new file mode 100644
index 0000000000000000000000000000000000000000..ceb4240798302a90da037143f586e1f8b7503aae
--- /dev/null
+++ b/Source/Filtering/imstkSurfaceMeshCut.h
@@ -0,0 +1,129 @@
+/*=========================================================================
+
+ Library: iMSTK
+
+ Copyright (c) Kitware, Inc. & Center for Modeling, Simulation,
+ & Imaging in Medicine, Rensselaer Polytechnic Institute.
+
+ Licensed under the Apache License, Version 2.0 (the "License");
+ you may not use this file except in compliance with the License.
+ You may obtain a copy of the License at
+
+ http://www.apache.org/licenses/LICENSE-2.0.txt
+
+ Unless required by applicable law or agreed to in writing, software
+ distributed under the License is distributed on an "AS IS" BASIS,
+ WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ See the License for the specific language governing permissions and
+ limitations under the License.
+
+=========================================================================*/
+#pragma once
+
+#include "imstkGeometryAlgorithm.h"
+#include "imstkSurfaceMesh.h"
+
+namespace imstk
+{
+class Plane;
+
+// vertex on the plane (0), positive side (+1), negative side (-1)
+// pt0 and pt1 follows the triangle's indexing order when tri is presented
+enum class CutType
+{
+ NONE = 0,
+ /* triangle is not cut through
+ * pt0 (-+1)
+ * / \
+ * / \
+ * / tri \
+ * pt1 (+-1)------(?)
+ */
+ EDGE,
+ /*
+ * (-+1)
+ * / \
+ * / \
+ * / tri \
+ * (-+1)------(0) pt0
+ */
+ VERT,
+ /*
+ * (+-1) pt1
+ * / \
+ * / \
+ * / tri \
+ * (-+1)------(+-1) pt0
+ */
+ EDGE_EDGE,
+ /*
+ * pt0 (+-1)
+ * / \
+ * / \
+ * / tri \
+ * pt1 (-+1)------(0)
+ */
+ EDGE_VERT,
+ /*
+ * pt0 (0)-----(+-1)
+ * / \ /
+ * / \ /
+ * / \ /
+ * (-+1)------(0) pt1
+ */
+ VERT_VERT
+};
+
+struct CutData
+{
+ public:
+ Vec3d cutCoords[2];
+ int triId = -1;
+ int ptIds[2] = { -1, -1 };
+ CutType cutType = CutType::NONE;
+};
+
+///
+/// \class SurfaceMeshCut
+///
+/// \brief This filter cuts the triangles of a SurfaceMesh into smaller
+/// triangles using input cutData
+/// \todo: test
+///
+class SurfaceMeshCut : public GeometryAlgorithm
+{
+public:
+ SurfaceMeshCut();
+ virtual ~SurfaceMeshCut() override = default;
+
+public:
+ std::shared_ptr<SurfaceMesh> getOutputMesh();
+ void setInputMesh(std::shared_ptr<SurfaceMesh> inputSurf);
+
+ imstkGetMacro(CutData, std::shared_ptr<std::vector<CutData>>);
+ imstkSetMacro(CutData, std::shared_ptr<std::vector<CutData>>);
+ imstkGetMacro(Plane, std::shared_ptr<Plane>);
+ imstkSetMacro(Plane, std::shared_ptr<Plane>);
+ imstkGetMacro(Epsilon, double);
+ imstkSetMacro(Epsilon, double);
+
+protected:
+ void requestUpdate() override;
+ void refinement(std::shared_ptr<SurfaceMesh> outputSurf,
+ std::map<int, bool>& cutVerts);
+
+ void splitVerts(std::shared_ptr<SurfaceMesh> outputSurf,
+ std::map<int, bool>& cutVerts);
+ int pointOnPlaneSide(Vec3d pt);
+ bool vertexOnBoundary(std::shared_ptr<VecDataArray<int, 3>> triangleIndices,
+ std::set<int>& triSet);
+
+ void generateCutData(std::shared_ptr<Plane> plane, std::shared_ptr<SurfaceMesh> inputSurf);
+// void generateCutData(triangle);
+
+private:
+ std::shared_ptr<std::vector<CutData>> m_CutData = std::make_shared<std::vector<CutData>>();
+ std::shared_ptr<Plane> m_Plane = std::make_shared<Plane>();
+ double m_Epsilon = 1;
+};
+}
\ No newline at end of file