diff --git a/Examples/Audio/AudioExample.cpp b/Examples/Audio/AudioExample.cpp
index f804b62f8d917240e572a7ddfe7b4d4bef3c7713..d4d425a1e77a08a1543655f56069e0c5b34b341c 100644
--- a/Examples/Audio/AudioExample.cpp
+++ b/Examples/Audio/AudioExample.cpp
@@ -81,7 +81,7 @@ playMusic(const std::string& filename)
 #ifdef iMSTK_AUDIO_ENABLED
     // Load an ogg music file
     sf::Music music;
-        
+
     CHECK(music.openFromFile(filename)) << "playMusic: Could not open the input music file: " << filename;
 
     // Display music informations
diff --git a/Examples/DeformableBody/DeformableBodyExample.cpp b/Examples/DeformableBody/DeformableBodyExample.cpp
index 85a89d67cfcc6fb88d66281e45eac1b4171ee972..0308baaed5fb22d2f653d047dbfa66047f19c360 100644
--- a/Examples/DeformableBody/DeformableBodyExample.cpp
+++ b/Examples/DeformableBody/DeformableBodyExample.cpp
@@ -51,7 +51,7 @@ main()
 
     // Load a tetrahedral mesh
     auto tetMesh = MeshIO::read(iMSTK_DATA_ROOT "/asianDragon/asianDragon.veg");
-    
+
     CHECK(tetMesh) << "Could not read mesh from file.";
 
     // Extract the surface mesh
@@ -59,7 +59,7 @@ main()
     auto volTetMesh = std::dynamic_pointer_cast<TetrahedralMesh>(tetMesh);
 
     CHECK(volTetMesh) << "Dynamic pointer cast from PointSet to TetrahedralMesh failed!";
-    
+
     volTetMesh->extractSurfaceMesh(surfMesh, true);
 
     // Construct a map
diff --git a/Examples/GeometryProcessing/CMakeLists.txt b/Examples/GeometryProcessing/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..1101501cc7251c1daa7002ca3fc2f81203a7d354
--- /dev/null
+++ b/Examples/GeometryProcessing/CMakeLists.txt
@@ -0,0 +1,41 @@
+###########################################################################
+#
+# 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.
+#
+###########################################################################
+
+project(Example-GeometryProcessing)
+
+#-----------------------------------------------------------------------------
+# Create executable
+#-----------------------------------------------------------------------------
+imstk_add_executable(${PROJECT_NAME} GeometryProcessingExample.cpp)
+
+#-----------------------------------------------------------------------------
+# Add the target to Examples folder
+#-----------------------------------------------------------------------------
+SET_TARGET_PROPERTIES (${PROJECT_NAME} PROPERTIES FOLDER Examples)
+
+#-----------------------------------------------------------------------------
+# Link libraries to executable
+#-----------------------------------------------------------------------------
+target_link_libraries(${PROJECT_NAME} SimulationManager apiUtilities)
+
+#-----------------------------------------------------------------------------
+# Set MSVC working directory to the install/bin directory
+#-----------------------------------------------------------------------------
+if(MSVC) # Configure running executable out of MSVC
+  set_property(TARGET ${PROJECT_NAME} PROPERTY VS_DEBUGGER_WORKING_DIRECTORY "${iMSTK_INSTALL_BIN_DIR}")
+endif()
diff --git a/Examples/GeometryProcessing/GeometryProcessingExample.cpp b/Examples/GeometryProcessing/GeometryProcessingExample.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..7fd2c362502be9bfe834904a545f0907ea0bb727
--- /dev/null
+++ b/Examples/GeometryProcessing/GeometryProcessingExample.cpp
@@ -0,0 +1,88 @@
+/*=========================================================================
+
+   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 "imstkSimulationManager.h"
+#include "imstkSceneObject.h"
+#include "imstkPlane.h"
+#include "imstkCylinder.h"
+#include "imstkCube.h"
+#include "imstkAPIUtilities.h"
+#include "imstkGeometryUtilities.h"
+
+using namespace imstk;
+
+///
+/// \brief This example demonstrates the geometry transforms in imstk
+///
+int
+main()
+{
+    // simManager and Scene
+    auto simManager = std::make_shared<SimulationManager>();
+    auto scene      = simManager->createNewScene("GeometryTransforms");
+
+    auto coarseTetMesh  = std::dynamic_pointer_cast<TetrahedralMesh>(MeshIO::read(iMSTK_DATA_ROOT "/asianDragon/asianDragon.veg"));
+    auto coarseSurfMesh = std::make_shared<SurfaceMesh>();
+    coarseTetMesh->extractSurfaceMesh(coarseSurfMesh, true);
+
+    std::shared_ptr<SurfaceMesh> fineSurfaceMesh(std::move(GeometryUtils::loopSubdivideSurfaceMesh(*coarseSurfMesh.get())));
+
+    fineSurfaceMesh->translate(Vec3d(0., -5., 0.), Geometry::TransformType::ConcatenateToTransform);
+    coarseSurfMesh->translate(Vec3d(0., 5., 0.), Geometry::TransformType::ConcatenateToTransform);
+
+    auto material0 = std::make_shared<RenderMaterial>();
+    material0->setDisplayMode(RenderMaterial::DisplayMode::WIREFRAME_SURFACE);
+    auto surfMeshModel0 = std::make_shared<VisualModel>(coarseSurfMesh);
+    surfMeshModel0->setRenderMaterial(material0);
+
+    auto sceneObj0 = std::make_shared<VisualObject>("coarse Mesh");
+    sceneObj0->addVisualModel(surfMeshModel0);
+
+    scene->addSceneObject(sceneObj0);
+
+    auto material = std::make_shared<RenderMaterial>();
+    material->setDebugColor(imstk::Color::Red);
+    material->setDisplayMode(RenderMaterial::DisplayMode::WIREFRAME_SURFACE);
+    auto surfMeshModel = std::make_shared<VisualModel>(fineSurfaceMesh);
+    surfMeshModel->setRenderMaterial(material);
+
+    auto sceneObj = std::make_shared<VisualObject>("fine Mesh");
+    sceneObj->addVisualModel(surfMeshModel);
+
+    scene->addSceneObject(sceneObj);
+
+    // Set Camera configuration
+    auto cam = scene->getCamera();
+    cam->setPosition(Vec3d(0, 12, 12));
+    cam->setFocalPoint(Vec3d(0, 0, 0));
+
+    // Light
+    auto light = std::make_shared<DirectionalLight>("light");
+    light->setFocalPoint(Vec3d(5, -8, -5));
+    light->setIntensity(1);
+    scene->addLight(light);
+
+    // Run
+    simManager->setActiveScene(scene);
+    simManager->start(SimulationStatus::running);
+
+    return 0;
+}
diff --git a/Source/Collision/CollisionHandling/imstkPenaltyCH.cpp b/Source/Collision/CollisionHandling/imstkPenaltyCH.cpp
index e4a097affe8ac66e2bb8a1fb184546c1f0614c38..f36864259f0ecb1ff970c2d88ee8f47239dfd84e 100644
--- a/Source/Collision/CollisionHandling/imstkPenaltyCH.cpp
+++ b/Source/Collision/CollisionHandling/imstkPenaltyCH.cpp
@@ -95,7 +95,7 @@ PenaltyCH::computeContactForcesDiscreteDeformable(const std::shared_ptr<Deformab
     }
 
     CHECK(deformableObj) << "PenaltyRigidCH::computeContactForcesDiscreteDeformable error: "
-        << m_object->getName() << " is not a deformable object.";
+                         << m_object->getName() << " is not a deformable object.";
 
     // Get current force vector
     auto&       force     = deformableObj->getContactForce();
diff --git a/Source/Collision/CollisionHandling/imstkPickingCH.cpp b/Source/Collision/CollisionHandling/imstkPickingCH.cpp
index 4fd60c6fcdc4d6e9717c5adc635879caec33059c..b8328b33d562391be71ab14d27c79969d2a394cd 100644
--- a/Source/Collision/CollisionHandling/imstkPickingCH.cpp
+++ b/Source/Collision/CollisionHandling/imstkPickingCH.cpp
@@ -38,7 +38,7 @@ PickingCH::processCollisionData()
 {
     CHECK(m_object) << "PickingCH::handleCollision error: "
                     << "no picking collision handling available the object";
-    
+
     this->addPickConstraints(m_object);
 }
 
@@ -51,10 +51,10 @@ PickingCH::addPickConstraints(std::shared_ptr<DeformableObject> deformableObj)
     {
         return;
     }
-    
+
     CHECK(deformableObj) << "PenaltyRigidCH::addPickConstraints error: "
-                              << " not a deformable object.";
-        
+                         << " not a deformable object.";
+
     const auto& Uprev = deformableObj->getDisplacements();
     const auto& Vprev = deformableObj->getVelocities();
 
diff --git a/Source/Collision/Testing/imstkTetraToTetraCDTest.cpp b/Source/Collision/Testing/imstkTetraToTetraCDTest.cpp
index 31770cb52df685a9dfc65ecfa0213d504b957bcf..e13ebc41aad22395da76b3a853f94aa72233c3a5 100644
--- a/Source/Collision/Testing/imstkTetraToTetraCDTest.cpp
+++ b/Source/Collision/Testing/imstkTetraToTetraCDTest.cpp
@@ -44,9 +44,9 @@ loadMesh(std::string externalDataSuffix)
     std::string                      file = iMSTK_DATA_ROOT + externalDataSuffix;
     std::shared_ptr<TetrahedralMesh> volMesh
         = std::static_pointer_cast<TetrahedralMesh>(imstk::MeshIO::read(file));
-    
+
     CHECK(volMesh) << "Failed to read a volumetric mesh file : " << file;
-    
+
     return volMesh;
 }
 
diff --git a/Source/Collision/imstkInteractionPair.cpp b/Source/Collision/imstkInteractionPair.cpp
index 598ab08a45301b059faefee2ea39a0b5dcc3db75..5026e8fd16ea9036a545f44cc21209029bb38c4b 100644
--- a/Source/Collision/imstkInteractionPair.cpp
+++ b/Source/Collision/imstkInteractionPair.cpp
@@ -55,7 +55,7 @@ InteractionPair::InteractionPair(std::shared_ptr<CollidingObject> A,
     if (CHAType != CollisionHandling::Type::None)
     {
         CHA = CollisionHandling::make_collision_handling(CHAType, CollisionHandling::Side::A, m_colData, A, B);
-        
+
         CHECK(CHA) << "InteractionPair error: can not instantiate collision handling for '"
                    << A->getName() << "' object.";
     }
@@ -65,9 +65,9 @@ InteractionPair::InteractionPair(std::shared_ptr<CollidingObject> A,
     if (CHBType != CollisionHandling::Type::None)
     {
         CHB = CollisionHandling::make_collision_handling(CHBType, CollisionHandling::Side::B, m_colData, B, A);
-        
+
         CHECK(CHB) << "InteractionPair error: can not instantiate collision handling for '"
-            << B->getName() << "' object.";
+                   << B->getName() << "' object.";
     }
 
     // Init interactionPair
diff --git a/Source/DynamicalModels/ObjectModels/imstkFEMDeformableBodyModel.cpp b/Source/DynamicalModels/ObjectModels/imstkFEMDeformableBodyModel.cpp
index 393f7915af3964d5351f3d11c8472d101277cbbd..4063f686670852590460024e8aebe047c23bfdac 100644
--- a/Source/DynamicalModels/ObjectModels/imstkFEMDeformableBodyModel.cpp
+++ b/Source/DynamicalModels/ObjectModels/imstkFEMDeformableBodyModel.cpp
@@ -74,7 +74,7 @@ FEMDeformableBodyModel::configure(const std::string& configFileName)
     vegaConfigFileOptions.addOptionOptional("gravity", &m_FEModelConfig->m_gravity, m_FEModelConfig->m_gravity);
 
     // Parse the configuration file
-    CHECK(vegaConfigFileOptions.parseOptions(configFileName.data()) == 0) 
+    CHECK(vegaConfigFileOptions.parseOptions(configFileName.data()) == 0)
         << "ForceModelConfig::parseConfig - Unable to load the configuration file";
 
     // get the root directory of the boundary file name
@@ -386,9 +386,9 @@ FEMDeformableBodyModel::initializeDampingMatrix()
 bool
 FEMDeformableBodyModel::initializeTangentStiffness()
 {
-    CHECK(m_internalForceModel) 
+    CHECK(m_internalForceModel)
         << "DeformableBodyModel::initializeTangentStiffness: Tangent stiffness cannot be initialized without force model";
-    
+
     vega::SparseMatrix* matrix;
     m_internalForceModel->getTangentStiffnessMatrixTopology(&matrix);
 
diff --git a/Source/DynamicalModels/ObjectModels/imstkPbdModel.cpp b/Source/DynamicalModels/ObjectModels/imstkPbdModel.cpp
index aaf4ed1b5a20d3a24f8f0ac99ecce10026643526..a981dd3f7edbdc73e73087dae7f29cfeed45690d 100644
--- a/Source/DynamicalModels/ObjectModels/imstkPbdModel.cpp
+++ b/Source/DynamicalModels/ObjectModels/imstkPbdModel.cpp
@@ -33,6 +33,8 @@
 #include "imstkPbdFEHexConstraint.h"
 #include "imstkPbdConstantDensityConstraint.h"
 #include "imstkParallelUtils.h"
+#include "imstkGeometryUtilities.h"
+#include "imstkMeshIO.h"
 
 #include <g3log/g3log.hpp>
 
@@ -178,7 +180,7 @@ bool
 PbdModel::initializeFEMConstraints(PbdFEMConstraint::MaterialType type)
 {
     // Check if constraint type matches the mesh type
-    CHECK(m_mesh->getType() == Geometry::Type::TetrahedralMesh) 
+    CHECK(m_mesh->getType() == Geometry::Type::TetrahedralMesh)
         << "FEM Tetrahedral constraint should come with tetrahedral mesh";
 
     // Create constraints
@@ -296,8 +298,8 @@ bool
 PbdModel::initializeAreaConstraints(const double stiffness)
 {
     // check if constraint type matches the mesh type
-    CHECK (m_mesh->getType() == Geometry::Type::SurfaceMesh)
-        << "Area constraint should come with a triangular mesh";    
+    CHECK(m_mesh->getType() == Geometry::Type::SurfaceMesh)
+        << "Area constraint should come with a triangular mesh";
 
     // ok, now create constraints
     const auto& triMesh  = std::static_pointer_cast<SurfaceMesh>(m_mesh);
@@ -439,13 +441,12 @@ bool
 PbdModel::initializeConstantDensityConstraint(const double stiffness)
 {
     // check if constraint type matches the mesh type
-    CHECK (m_mesh->getType() == Geometry::Type::SurfaceMesh
+    CHECK(m_mesh->getType() == Geometry::Type::SurfaceMesh
         || m_mesh->getType() == Geometry::Type::TetrahedralMesh
         || m_mesh->getType() == Geometry::Type::LineMesh
         || m_mesh->getType() == Geometry::Type::HexahedralMesh
         || m_mesh->getType() == Geometry::Type::PointSet)
         << "Constant constraint should come with a mesh!";
-    
 
     auto c = std::make_shared<PbdConstantDensityConstraint>();
     c->initConstraint(*this, stiffness);
diff --git a/Source/DynamicalModels/ObjectModels/imstkSPHKernels.h b/Source/DynamicalModels/ObjectModels/imstkSPHKernels.h
index 1573a16228685d7c8da7d651d0f83db979df39d3..3bc038819f4423d241c8c84bacdf4a043107c3f7 100644
--- a/Source/DynamicalModels/ObjectModels/imstkSPHKernels.h
+++ b/Source/DynamicalModels/ObjectModels/imstkSPHKernels.h
@@ -236,8 +236,8 @@ public:
 
         CHECK(N != 2) << "Unimplemented function";
 
-        m_k = Real(32.0) / (PI * std::pow(m_radius, 9));
-        m_c = std::pow(m_radius, 6) / Real(64.0);
+        m_k  = Real(32.0) / (PI * std::pow(m_radius, 9));
+        m_c  = std::pow(m_radius, 6) / Real(64.0);
         m_W0 = W(VecXr::Zero());
     }
 
@@ -323,7 +323,7 @@ public:
 
         CHECK(N != 2) << "Unimplemented function";
 
-        m_k = Real(0.007 / std::pow(m_radius, 3.25));        
+        m_k  = Real(0.007 / std::pow(m_radius, 3.25));
         m_W0 = W(VecXr::Zero());
     }
 
diff --git a/Source/DynamicalModels/ObjectStates/imstkSPHState.cpp b/Source/DynamicalModels/ObjectStates/imstkSPHState.cpp
index a10340b54286995d95b773461675cf26935099c7..4acf9de19853cf0d1f81934cf105258e7d7ec857 100644
--- a/Source/DynamicalModels/ObjectStates/imstkSPHState.cpp
+++ b/Source/DynamicalModels/ObjectStates/imstkSPHState.cpp
@@ -56,9 +56,9 @@ void
 SPHSimulationState::initializeData()
 {
     CHECK(m_KinematicState) << "SPH basic state has not been initialized";
-    
+
     size_t numParticles = m_KinematicState->getNumParticles();
-    
+
     m_Normals.resize(numParticles);
     m_Densities.resize(numParticles);
     m_Accels.resize(numParticles);
diff --git a/Source/Geometry/Analytic/imstkAnalyticalGeometry.cpp b/Source/Geometry/Analytic/imstkAnalyticalGeometry.cpp
index 5fd1b4c7c82c1f8ccda63880054b5e025930a01a..2c5ac36dd563c7d8ef27230eb57cd29af94171b9 100644
--- a/Source/Geometry/Analytic/imstkAnalyticalGeometry.cpp
+++ b/Source/Geometry/Analytic/imstkAnalyticalGeometry.cpp
@@ -102,7 +102,7 @@ AnalyticalGeometry::applyRotation(const Mat3d r)
 }
 
 void
-AnalyticalGeometry::updatePostTransformData()
+AnalyticalGeometry::updatePostTransformData() const
 {
     m_orientationAxisPostTransform = m_transform.rotation() * m_orientationAxis;
     m_orientationAxisPostTransform.normalize();
diff --git a/Source/Geometry/Analytic/imstkAnalyticalGeometry.h b/Source/Geometry/Analytic/imstkAnalyticalGeometry.h
index e643ca5f13a93940781c3a3fe32a348e732b48f9..e346bf6b580d43f86c36879c6264a5d86d44bec7 100644
--- a/Source/Geometry/Analytic/imstkAnalyticalGeometry.h
+++ b/Source/Geometry/Analytic/imstkAnalyticalGeometry.h
@@ -61,12 +61,12 @@ protected:
 
     void applyTranslation(const Vec3d t) override;
     void applyRotation(const Mat3d r) override;
-    virtual void updatePostTransformData() override;
+    virtual void updatePostTransformData() const override;
 
-    Vec3d m_position = WORLD_ORIGIN;                  ///> position
-    Vec3d m_positionPostTransform = WORLD_ORIGIN;     ///> position once transform applied
+    Vec3d m_position = WORLD_ORIGIN;                          ///> position
+    mutable Vec3d m_positionPostTransform = WORLD_ORIGIN;     ///> position once transform applied
 
-    Vec3d m_orientationAxis = UP_VECTOR;              ///> orientation
-    Vec3d m_orientationAxisPostTransform = UP_VECTOR; ///> orientation once transform applied
+    Vec3d m_orientationAxis = UP_VECTOR;                      ///> orientation
+    mutable Vec3d m_orientationAxisPostTransform = UP_VECTOR; ///> orientation once transform applied
 };
 } //imstk
diff --git a/Source/Geometry/Analytic/imstkCapsule.cpp b/Source/Geometry/Analytic/imstkCapsule.cpp
index a2634283db9d3d4de0c7e2feca6a860bda238ab1..bbdc898ce38dad1fb01eb244ce7d6ce5d73e447d 100644
--- a/Source/Geometry/Analytic/imstkCapsule.cpp
+++ b/Source/Geometry/Analytic/imstkCapsule.cpp
@@ -89,7 +89,7 @@ Capsule::applyScaling(const double s)
 }
 
 void
-Capsule::updatePostTransformData()
+Capsule::updatePostTransformData() const
 {
     if (m_transformApplied)
     {
diff --git a/Source/Geometry/Analytic/imstkCapsule.h b/Source/Geometry/Analytic/imstkCapsule.h
index 5c11ffbd9e87ac81df0d16fb9fb44cd74374bcba..c86f1034d640be64044b2320b2ace615bf2ae78e 100644
--- a/Source/Geometry/Analytic/imstkCapsule.h
+++ b/Source/Geometry/Analytic/imstkCapsule.h
@@ -73,11 +73,11 @@ protected:
     friend class VTKCapsuleRenderDelegate;
 
     void applyScaling(const double s) override;
-    void updatePostTransformData() override;
+    void updatePostTransformData() const override;
 
-    double m_radius = 1.0;              ///> Radius of the hemispheres at the end of the capsule
-    double m_radiusPostTransform = 1.0; ///> Radius after transform
-    double m_length = 1.0;              ///> Length between the centers of two hemispheres
-    double m_lengthPostTransform = 1.0; ///> Length after transform
+    double m_radius = 1.0;                      ///> Radius of the hemispheres at the end of the capsule
+    mutable double m_radiusPostTransform = 1.0; ///> Radius after transform
+    double m_length = 1.0;                      ///> Length between the centers of two hemispheres
+    mutable double m_lengthPostTransform = 1.0; ///> Length after transform
 };
 } // imstk
diff --git a/Source/Geometry/Analytic/imstkCube.cpp b/Source/Geometry/Analytic/imstkCube.cpp
index 63c1cd78e9d488b321745ee730c9332c85f5d738..98c8d9c5d7e790742612d81d5d0563a8a2f7a723 100644
--- a/Source/Geometry/Analytic/imstkCube.cpp
+++ b/Source/Geometry/Analytic/imstkCube.cpp
@@ -71,7 +71,7 @@ Cube::applyScaling(const double s)
 }
 
 void
-Cube::updatePostTransformData()
+Cube::updatePostTransformData() const
 {
     if (m_transformApplied)
     {
diff --git a/Source/Geometry/Analytic/imstkCube.h b/Source/Geometry/Analytic/imstkCube.h
index e36a6ca2f40d5dba662cf2cabe80235824d97625..a50de6ecefef19cca1980237f86d4258f45037fa 100644
--- a/Source/Geometry/Analytic/imstkCube.h
+++ b/Source/Geometry/Analytic/imstkCube.h
@@ -60,9 +60,9 @@ protected:
     friend class VTKCubeRenderDelegate;
 
     void applyScaling(const double s) override;
-    void updatePostTransformData() override;
+    void updatePostTransformData() const override;
 
-    double m_width = 1.0;               ///> Width of the cube
-    double m_widthPostTransform = 1.0;  ///> Width of the cube once transform applied
+    double m_width = 1.0;                      ///> Width of the cube
+    mutable double m_widthPostTransform = 1.0; ///> Width of the cube once transform applied
 };
 }
diff --git a/Source/Geometry/Analytic/imstkCylinder.cpp b/Source/Geometry/Analytic/imstkCylinder.cpp
index 38a1e253c11898450a981f7a6d13b30606d6b81c..8e47fefe1947fddd4d706f2bf2a60b674ff27bcc 100644
--- a/Source/Geometry/Analytic/imstkCylinder.cpp
+++ b/Source/Geometry/Analytic/imstkCylinder.cpp
@@ -65,7 +65,7 @@ void
 Cylinder::setRadius(const double r)
 {
     CHECK(r > 0) << "Cylinder::setRadius error: radius should be positive.";
-    
+
     if (m_radius == r)
     {
         return;
@@ -101,7 +101,7 @@ Cylinder::applyScaling(const double s)
 }
 
 void
-Cylinder::updatePostTransformData()
+Cylinder::updatePostTransformData() const
 {
     if (m_transformApplied)
     {
diff --git a/Source/Geometry/Analytic/imstkCylinder.h b/Source/Geometry/Analytic/imstkCylinder.h
index 5f7e02cc5206d128a295fba7d5fb5cc713a0753e..7c190d68e927dccb19bdf233c51edfcf3607d30e 100644
--- a/Source/Geometry/Analytic/imstkCylinder.h
+++ b/Source/Geometry/Analytic/imstkCylinder.h
@@ -73,11 +73,11 @@ protected:
     friend class VTKCylinderRenderDelegate;
 
     void applyScaling(const double s) override;
-    void updatePostTransformData() override;
+    void updatePostTransformData() const override;
 
-    double m_radius = 1.;              ///> Radius of the cylinder
-    double m_length = 1.;              ///> Length of the cylinder
-    double m_radiusPostTransform = 1.; ///> Radius of the cylinder oncee transform applied
-    double m_lengthPostTransform = 1.; ///> Length of the cylinder onc transform applied
+    double m_radius = 1.;                      ///> Radius of the cylinder
+    double m_length = 1.;                      ///> Length of the cylinder
+    mutable double m_radiusPostTransform = 1.; ///> Radius of the cylinder oncee transform applied
+    mutable double m_lengthPostTransform = 1.; ///> Length of the cylinder onc transform applied
 };
 } // imstk
diff --git a/Source/Geometry/Analytic/imstkPlane.cpp b/Source/Geometry/Analytic/imstkPlane.cpp
index 1d1822a7567a7e19344460933382090aca494316..d298d07efb80080e55cfa213bc9d2225216e4681 100644
--- a/Source/Geometry/Analytic/imstkPlane.cpp
+++ b/Source/Geometry/Analytic/imstkPlane.cpp
@@ -89,7 +89,7 @@ Plane::applyScaling(const double s)
 }
 
 void
-Plane::updatePostTransformData()
+Plane::updatePostTransformData() const
 {
     if (m_transformApplied)
     {
diff --git a/Source/Geometry/Analytic/imstkPlane.h b/Source/Geometry/Analytic/imstkPlane.h
index dd7c4a6f61ceaf96d1583637b2103c591bee5c84..c4c536ba435553823796f9082a6c20cdba167218 100644
--- a/Source/Geometry/Analytic/imstkPlane.h
+++ b/Source/Geometry/Analytic/imstkPlane.h
@@ -74,9 +74,9 @@ protected:
     friend class VTKPlaneRenderDelegate;
 
     void applyScaling(const double s) override;
-    void updatePostTransformData() override;
+    void updatePostTransformData() const override;
 
-    double m_width = 1.0;               ///> Width of the plane
-    double m_widthPostTransform = 1.0;  ///> Width of the plane once transform applied
+    double m_width = 1.0;                      ///> Width of the plane
+    mutable double m_widthPostTransform = 1.0; ///> Width of the plane once transform applied
 };
 } // imstk
diff --git a/Source/Geometry/Analytic/imstkSphere.cpp b/Source/Geometry/Analytic/imstkSphere.cpp
index 05d4d9dc7818e399a2c7a678c3705eff94f39749..65af6c65b4aa4811870f89d2197f8ec6f6c2448a 100644
--- a/Source/Geometry/Analytic/imstkSphere.cpp
+++ b/Source/Geometry/Analytic/imstkSphere.cpp
@@ -88,7 +88,7 @@ Sphere::applyScaling(const double s)
 }
 
 void
-Sphere::updatePostTransformData()
+Sphere::updatePostTransformData() const
 {
     if (m_transformApplied)
     {
diff --git a/Source/Geometry/Analytic/imstkSphere.h b/Source/Geometry/Analytic/imstkSphere.h
index 4ae41fd6406aafacd3860c6741df1a4590315d25..aab7b3da6be7ebb8a859e807d827a1fdd4ff6ea8 100644
--- a/Source/Geometry/Analytic/imstkSphere.h
+++ b/Source/Geometry/Analytic/imstkSphere.h
@@ -68,9 +68,9 @@ protected:
     friend class VTKSphereRenderDelegate;
 
     void applyScaling(const double s) override;
-    void updatePostTransformData() override;
+    void updatePostTransformData() const override;
 
-    double m_radius = 1.0;              ///> Radius of the sphere
-    double m_radiusPostTransform = 1.0; ///> Radius of the sphere once transform applied
+    double m_radius = 1.0;                      ///> Radius of the sphere
+    mutable double m_radiusPostTransform = 1.0; ///> Radius of the sphere once transform applied
 };
 } // imstk
diff --git a/Source/Geometry/Decal/imstkDecalPool.h b/Source/Geometry/Decal/imstkDecalPool.h
index 3f20ae0d01469f034cdaee2d3d38ade2e6cc0061..1a2ceed70095612fbe10a3ba4145c10c14bb3714 100644
--- a/Source/Geometry/Decal/imstkDecalPool.h
+++ b/Source/Geometry/Decal/imstkDecalPool.h
@@ -61,7 +61,7 @@ protected:
     void applyTranslation(const Vec3d t) override {}
     void applyRotation(const Mat3d r) override {}
     void applyScaling(const double s) override {}
-    virtual void updatePostTransformData() override {}
+    virtual void updatePostTransformData() const override {}
 
     unsigned int m_maxNumDecals;
     unsigned int m_numDecals = 0;
diff --git a/Source/Geometry/Map/imstkIdentityMap.cpp b/Source/Geometry/Map/imstkIdentityMap.cpp
index b53441d90072aa1f82f44ea3c72ea62e444ba612..6ff7f5bd35279285a8e8a6da8a8e2f6925852979 100644
--- a/Source/Geometry/Map/imstkIdentityMap.cpp
+++ b/Source/Geometry/Map/imstkIdentityMap.cpp
@@ -35,7 +35,7 @@ IdentityMap::apply()
 
     // Check geometries
     CHECK(m_master && m_slave) << "Identity map is being applied without valid geometries";
-    
+
     // Set the follower mesh configuration to be same as that of master
     m_slave->setTranslation(m_master->getTranslation());
     m_slave->setRotation(m_master->getRotation());
diff --git a/Source/Geometry/Mesh/imstkImageData.cpp b/Source/Geometry/Mesh/imstkImageData.cpp
index 21845bbb106c5cd3a078c78910e37cc70194dbeb..34d17a2a134c0732d5dcd54beaf053bafcdd88de 100644
--- a/Source/Geometry/Mesh/imstkImageData.cpp
+++ b/Source/Geometry/Mesh/imstkImageData.cpp
@@ -123,7 +123,7 @@ ImageData::applyRotation(const Mat3d r)
 }
 
 void
-ImageData::updatePostTransformData()
+ImageData::updatePostTransformData() const
 {
     if (m_transformApplied || !this->m_data)
     {
diff --git a/Source/Geometry/Mesh/imstkImageData.h b/Source/Geometry/Mesh/imstkImageData.h
index 4d5fa73d9e674949ef5259d7d30cd93407baa416..a90b1073ddac36e8b25b246d7abb81a472eaf088 100644
--- a/Source/Geometry/Mesh/imstkImageData.h
+++ b/Source/Geometry/Mesh/imstkImageData.h
@@ -83,6 +83,6 @@ protected:
 
     void applyScaling(const double s) override;
 
-    void updatePostTransformData() override;
+    void updatePostTransformData() const override;
 };
 } // imstk
diff --git a/Source/Geometry/Mesh/imstkPointSet.cpp b/Source/Geometry/Mesh/imstkPointSet.cpp
index 2f868ccec892ae8acbedda91630685e2c67e724f..b33d9af57cbcb7a54328e4f193d0293aeef13e44 100644
--- a/Source/Geometry/Mesh/imstkPointSet.cpp
+++ b/Source/Geometry/Mesh/imstkPointSet.cpp
@@ -113,7 +113,7 @@ PointSet::setVertexPositions(const StdVectorOfVec3d& vertices)
 }
 
 const StdVectorOfVec3d&
-PointSet::getVertexPositions(DataType type /* = DataType::PostTransform */)
+PointSet::getVertexPositions(DataType type /* = DataType::PostTransform */) const
 {
     if (type == DataType::PostTransform)
     {
@@ -136,7 +136,7 @@ PointSet::setVertexPosition(const size_t vertNum, const Vec3d& pos)
 }
 
 const Vec3d&
-PointSet::getVertexPosition(const size_t vertNum, DataType type)
+PointSet::getVertexPosition(const size_t vertNum, DataType type) const
 {
 #if defined(DEBUG) || defined(_DEBUG) || !defined(NDEBUG)
     LOG_IF(FATAL, (vertNum >= getVertexPositions().size())) << "Invalid index";
@@ -275,7 +275,7 @@ PointSet::applyScaling(const double s)
 }
 
 void
-PointSet::updatePostTransformData()
+PointSet::updatePostTransformData() const
 {
     if (m_transformApplied)
     {
diff --git a/Source/Geometry/Mesh/imstkPointSet.h b/Source/Geometry/Mesh/imstkPointSet.h
index 4b7b27b4a5e0fe6c2cb8d544c932df85bd358125..f47b909e989811956a17158c2334f4fe0beb1376 100644
--- a/Source/Geometry/Mesh/imstkPointSet.h
+++ b/Source/Geometry/Mesh/imstkPointSet.h
@@ -91,7 +91,7 @@ public:
     ///
     /// \brief Returns the vector of current positions of the mesh vertices
     ///
-    const StdVectorOfVec3d& getVertexPositions(DataType type = DataType::PostTransform);
+    const StdVectorOfVec3d& getVertexPositions(DataType type = DataType::PostTransform) const;
 
     ///
     /// \brief Set the current position of a vertex given its index to certain position (this is not a thread-safe method)
@@ -101,7 +101,7 @@ public:
     ///
     /// \brief Returns the position of a vertex given its index
     ///
-    const Vec3d& getVertexPosition(const size_t vertNum, DataType type = DataType::PostTransform);
+    const Vec3d& getVertexPosition(const size_t vertNum, DataType type = DataType::PostTransform) const;
 
     ///
     /// \brief Sets the displacements of mesh vertices from an array
@@ -177,11 +177,11 @@ protected:
     void applyTranslation(const Vec3d t) override;
     void applyRotation(const Mat3d r) override;
     void applyScaling(const double s) override;
-    void updatePostTransformData() override;
+    void updatePostTransformData() const override;
 
     StdVectorOfVec3d m_initialVertexPositions;                ///> Initial positions of vertices
     StdVectorOfVec3d m_vertexPositions;                       ///> Current positions of vertices
-    StdVectorOfVec3d m_vertexPositionsPostTransform;          ///> Positions of vertices after transform
+    mutable StdVectorOfVec3d m_vertexPositionsPostTransform;  ///> Positions of vertices after transform
 
     std::map<std::string, StdVectorOfVectorf> m_pointDataMap; ///> vector of data arrays per vertice
 
diff --git a/Source/Geometry/Mesh/imstkTetrahedralMesh.cpp b/Source/Geometry/Mesh/imstkTetrahedralMesh.cpp
index 5071485e55631d3343682dfaf835154619807d17..f0aa462fe57d3d6a31c5ef7763d97a8097f316cd 100644
--- a/Source/Geometry/Mesh/imstkTetrahedralMesh.cpp
+++ b/Source/Geometry/Mesh/imstkTetrahedralMesh.cpp
@@ -103,9 +103,9 @@ bool
 TetrahedralMesh::extractSurfaceMesh(std::shared_ptr<SurfaceMesh> surfaceMesh,
                                     const bool                   enforceWindingConsistency /* = false*/)
 {
-    CHECK(surfaceMesh) 
-            << "TetrahedralMesh::extractSurfaceMesh error: the surface mesh provided is not instantiated.";
- 
+    CHECK(surfaceMesh)
+        << "TetrahedralMesh::extractSurfaceMesh error: the surface mesh provided is not instantiated.";
+
     using triArray = SurfaceMesh::TriangleArray;
     const std::vector<triArray> facePattern = {
         triArray { { 0, 1, 2 } }, triArray { { 0, 1, 3 } }, triArray { { 0, 2, 3 } }, triArray { { 1, 2, 3 } }
diff --git a/Source/Geometry/Particles/imstkRenderParticles.h b/Source/Geometry/Particles/imstkRenderParticles.h
index 56b8c751f8fe0edca3940a14dbc0364c8f109500..923c41fb5de64de75669fe6b1e798e297e29f7b9 100644
--- a/Source/Geometry/Particles/imstkRenderParticles.h
+++ b/Source/Geometry/Particles/imstkRenderParticles.h
@@ -127,6 +127,6 @@ protected:
     void applyTranslation(const Vec3d t) override {}
     void applyRotation(const Mat3d r) override {}
     void applyScaling(const double s) override {}
-    virtual void updatePostTransformData() override {}
+    virtual void updatePostTransformData() const override {}
 };
 }
diff --git a/Source/Geometry/Reader/imstkMeshIO.cpp b/Source/Geometry/Reader/imstkMeshIO.cpp
index 7e3426f12b5a490072277e494a08462d6f434ced..0f80d3f403cb81d77dc03a8e46a12227173cb977 100644
--- a/Source/Geometry/Reader/imstkMeshIO.cpp
+++ b/Source/Geometry/Reader/imstkMeshIO.cpp
@@ -36,7 +36,7 @@ MeshIO::read(const std::string& filePath)
 {
     bool isDir;
 
-    CHECK(MeshIO::fileExists(filePath, isDir))<< "MeshIO::read error: file not found: " << filePath;
+    CHECK(MeshIO::fileExists(filePath, isDir)) << "MeshIO::read error: file not found: " << filePath;
 
     if (isDir)
     {
@@ -105,8 +105,8 @@ MeshIO::getFileType(const std::string& filePath)
     MeshFileType meshType = MeshFileType::UNKNOWN;
 
     std::string extString = filePath.substr(filePath.find_last_of(".") + 1);
-    
-    CHECK (!extString.empty()) << "MeshIO::getFileType error: invalid file name";
+
+    CHECK(!extString.empty()) << "MeshIO::getFileType error: invalid file name";
 
     if (extString == "vtk" || extString == "VTK")
     {
@@ -178,6 +178,7 @@ MeshIO::write(const std::shared_ptr<imstk::PointSet> imstkMesh, const std::strin
         return VegaMeshIO::write(imstkMesh, filePath, meshType);
         break;
     case MeshFileType::VTU:
+    case MeshFileType::VTK:
     case MeshFileType::VTP:
     case MeshFileType::STL:
     case MeshFileType::PLY:
diff --git a/Source/Geometry/Reader/imstkVTKMeshIO.cpp b/Source/Geometry/Reader/imstkVTKMeshIO.cpp
index d1a02f56057eeb3f0bef92f4d7638139e9ff760d..f8b9ec312868035d123effba2c146ef142d30522 100644
--- a/Source/Geometry/Reader/imstkVTKMeshIO.cpp
+++ b/Source/Geometry/Reader/imstkVTKMeshIO.cpp
@@ -20,6 +20,7 @@
 =========================================================================*/
 
 #include "imstkVTKMeshIO.h"
+#include "imstkGeometryUtilities.h"
 
 #include "vtkSmartPointer.h"
 #include "vtkGenericDataObjectReader.h"
@@ -34,6 +35,7 @@
 #include "vtkSTLWriter.h"
 #include "vtkFloatArray.h"
 #include "vtkTriangleFilter.h"
+#include "vtkPolyDataWriter.h"
 #include "vtkDICOMImageReader.h"
 #include "vtkNrrdReader.h"
 
@@ -94,7 +96,19 @@ VTKMeshIO::write(const std::shared_ptr<PointSet> imstkMesh, const std::string& f
         switch (meshType)
         {
         case MeshFileType::VTU:
-            return VTKMeshIO::writeVtkUnstructuredGrid(vMesh, filePath);
+            if (auto tetMesh = std::dynamic_pointer_cast<TetrahedralMesh>(vMesh))
+            {
+                VTKMeshIO::writeVtkUnstructuredGrid(*tetMesh.get(), filePath);
+            }
+            else if (auto hexMesh = std::dynamic_pointer_cast<HexahedralMesh>(vMesh))
+            {
+                VTKMeshIO::writeVtkUnstructuredGrid(*hexMesh.get(), filePath);
+            }
+            else
+            {
+                return false;
+            }
+
         default:
             LOG(WARNING) << "VTKMeshIO::write error: file type not supported for volumetric mesh.";
             return false;
@@ -105,16 +119,31 @@ VTKMeshIO::write(const std::shared_ptr<PointSet> imstkMesh, const std::string& f
         switch (meshType)
         {
         case MeshFileType::VTP:
-            return VTKMeshIO::writeVtkPolyData<vtkXMLPolyDataWriter>(sMesh, filePath);
+            return VTKMeshIO::writeVtkPolyData<vtkXMLPolyDataWriter>(*sMesh.get(), filePath);
         case MeshFileType::STL:
-            return VTKMeshIO::writeVtkPolyData<vtkSTLWriter>(sMesh, filePath);
+            return VTKMeshIO::writeVtkPolyData<vtkSTLWriter>(*sMesh.get(), filePath);
         case MeshFileType::PLY:
-            return VTKMeshIO::writeVtkPolyData<vtkPLYWriter>(sMesh, filePath);
+            return VTKMeshIO::writeVtkPolyData<vtkPLYWriter>(*sMesh.get(), filePath);
+        case MeshFileType::VTK:
+            return VTKMeshIO::writeVtkPolyData<vtkPolyDataWriter>(*sMesh.get(), filePath);
         default:
             LOG(WARNING) << "VTKMeshIO::write error: file type not supported for surface mesh.";
             return false;
         }
     }
+    else if (auto lMesh = std::dynamic_pointer_cast<LineMesh>(imstkMesh))
+    {
+        switch (meshType)
+        {
+        case MeshFileType::VTK:
+            return VTKMeshIO::writeVtkPolyData<vtkPolyDataWriter>(*lMesh.get(), filePath);
+        case MeshFileType::VTP:
+            return VTKMeshIO::writeVtkPolyData<vtkXMLPolyDataWriter>(*lMesh.get(), filePath);
+        default:
+            LOG(WARNING) << "vtkMeshIO::write error: file type not supported for line mesh.";
+            return false;
+        }
+    }
     else
     {
         LOG(WARNING) << "VTKMeshIO::write error: the provided mesh is not a surface or volumetric mesh.";
@@ -130,13 +159,13 @@ VTKMeshIO::readVtkGenericFormatData(const std::string& filePath)
     reader->SetFileName(filePath.c_str());
     reader->Update();
 
-    if (vtkPolyData* vtkMesh = reader->GetPolyDataOutput())
+    if (vtkSmartPointer<vtkPolyData> vtkMesh = reader->GetPolyDataOutput())
     {
-        return VTKMeshIO::convertVtkPolyDataToSurfaceMesh(vtkMesh);
+        return GeometryUtils::convertVtkPolyDataToSurfaceMesh(vtkMesh);
     }
     else if (vtkUnstructuredGrid* vtkMesh = reader->GetUnstructuredGridOutput())
     {
-        return VTKMeshIO::convertVtkUnstructuredGridToVolumetricMesh(vtkMesh);
+        return GeometryUtils::convertVtkUnstructuredGridToVolumetricMesh(vtkMesh);
     }
     else
     {
@@ -157,15 +186,34 @@ VTKMeshIO::readVtkPolyData(const std::string& filePath)
     triFilter->SetInputData(reader->GetOutput());
     triFilter->Update();
 
-    vtkPolyData* vtkMesh = triFilter->GetOutput();
-    return VTKMeshIO::convertVtkPolyDataToSurfaceMesh(vtkMesh);
+    vtkSmartPointer<vtkPolyData> vtkMesh = triFilter->GetOutput();
+    return GeometryUtils::convertVtkPolyDataToSurfaceMesh(vtkMesh);
 }
 
 template<typename WriterType>
 bool
-VTKMeshIO::writeVtkPolyData(const std::shared_ptr<SurfaceMesh> imstkMesh, const std::string& filePath)
+VTKMeshIO::writeVtkPolyData(const SurfaceMesh& imstkMesh, const std::string& filePath)
 {
-    vtkPolyData* vtkMesh = convertSurfaceMeshToVtkPolyData(imstkMesh);
+    vtkSmartPointer<vtkPolyData> vtkMesh = GeometryUtils::convertSurfaceMeshToVtkPolyData(imstkMesh);
+    if (!vtkMesh)  //conversion unsuccessful
+    {
+        return false;
+    }
+
+    auto writer = vtkSmartPointer<WriterType>::New();
+    writer->SetInputData(vtkMesh);
+    writer->SetFileName(filePath.c_str());
+    writer->Update();
+
+    return true;
+}
+
+template<typename WriterType>
+bool
+VTKMeshIO::writeVtkPolyData(const LineMesh& imstkMesh, const std::string& filePath)
+{
+    vtkSmartPointer<vtkPolyData> vtkMesh = GeometryUtils::convertLineMeshToVtkPolyData(imstkMesh);
+    vtkMesh->PrintSelf(std::cout, vtkIndent());
     if (!vtkMesh)  //conversion unsuccessful
     {
         return false;
@@ -175,7 +223,6 @@ VTKMeshIO::writeVtkPolyData(const std::shared_ptr<SurfaceMesh> imstkMesh, const
     writer->SetInputData(vtkMesh);
     writer->SetFileName(filePath.c_str());
     writer->Update();
-    vtkMesh->Delete();
 
     return true;
 }
@@ -188,8 +235,8 @@ VTKMeshIO::readVtkUnstructuredGrid(const std::string& filePath)
     reader->SetFileName(filePath.c_str());
     reader->Update();
 
-    vtkUnstructuredGrid* vtkMesh = reader->GetOutput();
-    return VTKMeshIO::convertVtkUnstructuredGridToVolumetricMesh(vtkMesh);
+    vtkSmartPointer<vtkUnstructuredGrid> vtkMesh = reader->GetOutput();
+    return GeometryUtils::convertVtkUnstructuredGridToVolumetricMesh(vtkMesh);
 }
 
 template<typename ReaderType>
@@ -227,24 +274,9 @@ VTKMeshIO::readVtkImageDataDICOM(const std::string& filePath)
 }
 
 bool
-VTKMeshIO::writeVtkUnstructuredGrid(const std::shared_ptr<VolumetricMesh> imstkMesh, const std::string& filePath)
+VTKMeshIO::writeVtkUnstructuredGrid(const TetrahedralMesh& tetMesh, const std::string& filePath)
 {
-    auto                 tMesh   = std::dynamic_pointer_cast<TetrahedralMesh>(imstkMesh);
-    auto                 hMesh   = std::dynamic_pointer_cast<HexahedralMesh>(imstkMesh);
-    vtkUnstructuredGrid* vtkMesh = nullptr;
-    if (tMesh)
-    {
-        vtkMesh = convertTetrahedralMeshToVtkUnstructuredGrid(tMesh);
-    }
-    else if (hMesh)
-    {
-        vtkMesh = convertHexahedralMeshToVtkUnstructuredGrid(hMesh);
-    }
-    else
-    {
-        LOG(WARNING) << "VTKMeshIO::writeVtkUnstructuredGrid error: mesh is neither tetrahedral nor hexahedral";
-        return false;
-    }
+    auto vtkMesh = GeometryUtils::convertTetrahedralMeshToVtkUnstructuredGrid(tetMesh);
 
     if (!vtkMesh)
     {
@@ -256,218 +288,26 @@ VTKMeshIO::writeVtkUnstructuredGrid(const std::shared_ptr<VolumetricMesh> imstkM
     writer->SetInputData(vtkMesh);
     writer->SetFileName(filePath.c_str());
     writer->Update();
-    vtkMesh->Delete();
 
     return true;
 }
 
-std::shared_ptr<SurfaceMesh>
-VTKMeshIO::convertVtkPolyDataToSurfaceMesh(vtkPolyData* vtkMesh)
-{
-    CHECK(vtkMesh) << "VTKMeshIO::convertVtkPolyDataToSurfaceMesh error: could not read with VTK reader.";
-
-    StdVectorOfVec3d vertices;
-    VTKMeshIO::copyVerticesFromVtk(vtkMesh->GetPoints(), vertices);
-
-    std::vector<SurfaceMesh::TriangleArray> triangles;
-    VTKMeshIO::copyCellsFromVtk<3>(vtkMesh->GetPolys(), triangles);
-
-    auto mesh = std::make_shared<SurfaceMesh>();
-    mesh->initialize(vertices, triangles, true);
-
-    // Point Data
-    std::map<std::string, StdVectorOfVectorf> dataMap;
-    VTKMeshIO::copyPointData(vtkMesh->GetPointData(), dataMap);
-    if (!dataMap.empty())
-    {
-        mesh->setPointDataMap(dataMap);
-    }
-
-    // Active Texture
-    if (auto pointData = vtkMesh->GetPointData())
-    {
-        if (auto tcoords = pointData->GetTCoords())
-        {
-            mesh->setDefaultTCoords(tcoords->GetName());
-        }
-    }
-
-    return mesh;
-}
-
-vtkPolyData*
-VTKMeshIO::convertSurfaceMeshToVtkPolyData(std::shared_ptr<SurfaceMesh> imstkMesh)
-{
-    vtkSmartPointer<vtkPoints> points = vtkSmartPointer<vtkPoints>::New();
-    VTKMeshIO::copyVerticesToVtk(imstkMesh->getVertexPositions(), points.Get());
-
-    vtkSmartPointer<vtkCellArray> polys = vtkSmartPointer<vtkCellArray>::New();
-    VTKMeshIO::copyCellsToVtk<3>(imstkMesh->getTrianglesVertices(), polys.Get());
-
-    vtkPolyData* polydata = vtkPolyData::New();
-    polydata->SetPoints(points);
-    polydata->SetPolys(polys);
-    return polydata;
-}
-
-vtkUnstructuredGrid*
-VTKMeshIO::convertTetrahedralMeshToVtkUnstructuredGrid(std::shared_ptr<TetrahedralMesh> imstkMesh)
-{
-    vtkSmartPointer<vtkPoints> points = vtkSmartPointer<vtkPoints>::New();
-    VTKMeshIO::copyVerticesToVtk(imstkMesh->getVertexPositions(), points.Get());
-
-    vtkSmartPointer<vtkCellArray> tetras = vtkSmartPointer<vtkCellArray>::New();
-    VTKMeshIO::copyCellsToVtk<4>(imstkMesh->getTetrahedraVertices(), tetras.Get());
-
-    vtkUnstructuredGrid* ug = vtkUnstructuredGrid::New();
-    ug->SetPoints(points);
-    ug->SetCells(VTK_TETRA, tetras);
-    return ug;
-}
-
-vtkUnstructuredGrid*
-VTKMeshIO::convertHexahedralMeshToVtkUnstructuredGrid(std::shared_ptr<HexahedralMesh> imstkMesh)
-{
-    vtkSmartPointer<vtkPoints> points = vtkSmartPointer<vtkPoints>::New();
-    VTKMeshIO::copyVerticesToVtk(imstkMesh->getVertexPositions(), points.Get());
-
-    vtkSmartPointer<vtkCellArray> bricks = vtkSmartPointer<vtkCellArray>::New();
-    VTKMeshIO::copyCellsToVtk<8>(imstkMesh->getHexahedraVertices(), bricks.Get());
-
-    vtkUnstructuredGrid* ug = vtkUnstructuredGrid::New();
-    ug->SetPoints(points);
-    ug->SetCells(VTK_HEXAHEDRON, bricks);
-    return ug;
-}
-
-std::shared_ptr<VolumetricMesh>
-VTKMeshIO::convertVtkUnstructuredGridToVolumetricMesh(vtkUnstructuredGrid* vtkMesh)
-{
-    CHECK(vtkMesh) << "VTKMeshIO::convertVtkUnstructuredGridToVolumetricMesh error: could not read with VTK reader.";
-
-    StdVectorOfVec3d vertices;
-    VTKMeshIO::copyVerticesFromVtk(vtkMesh->GetPoints(), vertices);
-
-    int cellType = vtkMesh->GetCellType(vtkMesh->GetNumberOfCells() - 1);
-    if (cellType == VTK_TETRA)
-    {
-        std::vector<TetrahedralMesh::TetraArray> cells;
-        VTKMeshIO::copyCellsFromVtk<4>(vtkMesh->GetCells(), cells);
-
-        auto mesh = std::make_shared<TetrahedralMesh>();
-        mesh->initialize(vertices, cells, false);
-        return mesh;
-    }
-    else if (cellType == VTK_HEXAHEDRON)
-    {
-        std::vector<HexahedralMesh::HexaArray> cells;
-        VTKMeshIO::copyCellsFromVtk<8>(vtkMesh->GetCells(), cells);
-
-        auto mesh = std::make_shared<HexahedralMesh>();
-        mesh->initialize(vertices, cells, false);
-        return mesh;
-    }
-    else
-    {
-        LOG(FATAL) << "VTKMeshIO::convertVtkUnstructuredGridToVolumetricMesh error: No support for vtkCellType="
-                   << cellType << ".";
-        return nullptr;
-    }
-}
-
-void
-VTKMeshIO::copyVerticesFromVtk(vtkPoints* points, StdVectorOfVec3d& vertices)
-{
-    CHECK(points) << "VTKMeshIO::copyVerticesFromVtk error: No points found.";
-
-    vertices.reserve(points->GetNumberOfPoints());
-    for (vtkIdType i = 0; i < points->GetNumberOfPoints(); ++i)
-    {
-        double pos[3];
-        points->GetPoint(i, pos);
-        vertices.emplace_back(pos[0], pos[1], pos[2]);
-    }
-}
-
-void
-VTKMeshIO::copyVerticesToVtk(const StdVectorOfVec3d& vertices, vtkPoints* points)
-{
-    CHECK(points) << "VTKMeshIO::copyVerticesToVtk error: No points found.";
-
-    points->SetNumberOfPoints(vertices.size());
-    for (size_t i = 0; i < vertices.size(); i++)
-    {
-        points->SetPoint(i, vertices[i][0], vertices[i][1], vertices[i][2]);
-    }
-}
-
-template<size_t dim>
-void
-VTKMeshIO::copyCellsToVtk(const std::vector<std::array<size_t, dim>>& cells, vtkCellArray* vtkCells)
-{
-    CHECK(vtkCells) << "VTKMeshIO::copyCellsToVtk error: No cells found.";
-
-    for (size_t i = 0; i < cells.size(); i++)
-    {
-        vtkCells->InsertNextCell(dim);
-        for (size_t k = 0; k < dim; k++)
-        {
-            vtkCells->InsertCellPoint(cells[i][k]);
-        }
-    }
-}
-
-template<size_t dim>
-void
-VTKMeshIO::copyCellsFromVtk(vtkCellArray* vtkCells, std::vector<std::array<size_t, dim>>& cells)
+bool
+VTKMeshIO::writeVtkUnstructuredGrid(const HexahedralMesh& hMesh, const std::string& filePath)
 {
-    CHECK(vtkCells) << "VTKMeshIO::copyCellsFromVtk error: No cells found.";
+    auto vtkMesh = GeometryUtils::convertHexahedralMeshToVtkUnstructuredGrid(hMesh);
 
-    cells.reserve(vtkCells->GetNumberOfCells());
-    vtkCells->InitTraversal();
-    auto                    vtkCell = vtkSmartPointer<vtkIdList>::New();
-    std::array<size_t, dim> cell;
-    while (vtkCells->GetNextCell(vtkCell))
+    if (!vtkMesh)
     {
-        if (vtkCell->GetNumberOfIds() != dim)
-        {
-            continue;
-        }
-        for (size_t i = 0; i < dim; ++i)
-        {
-            cell[i] = vtkCell->GetId(i);
-        }
-        cells.emplace_back(cell);
+        LOG(WARNING) << "VTKMeshIO::writeVtkUnstructuredGrid error: conversion unsuccessful";
+        return false;
     }
-}
 
-void
-VTKMeshIO::copyPointData(vtkPointData* pointData, std::map<std::string, StdVectorOfVectorf>& dataMap)
-{
-    if (!pointData)
-    {
-        return;
-    }
+    auto writer = vtkSmartPointer<vtkXMLUnstructuredGridWriter>::New();
+    writer->SetInputData(vtkMesh);
+    writer->SetFileName(filePath.c_str());
+    writer->Update();
 
-    for (int i = 0; i < pointData->GetNumberOfArrays(); ++i)
-    {
-        vtkDataArray*      array       = pointData->GetArray(i);
-        std::string        name        = array->GetName();
-        int                nbrOfComp   = array->GetNumberOfComponents();
-        vtkIdType          nbrOfTuples = array->GetNumberOfTuples();
-        StdVectorOfVectorf data;
-        for (vtkIdType j = 0; j < nbrOfTuples; ++j)
-        {
-            double* tupleData = new double[nbrOfComp];
-            array->GetTuple(j, tupleData);
-            Vectorf tuple(nbrOfComp);
-            for (int k = 0; k < nbrOfComp; k++)
-            {
-                tuple[k] = tupleData[k];
-            }
-            data.push_back(tuple);
-        }
-        dataMap[name] = data;
-    }
+    return true;
 }
 } // imstk
diff --git a/Source/Geometry/Reader/imstkVTKMeshIO.h b/Source/Geometry/Reader/imstkVTKMeshIO.h
index ac611e6b724176debfedf6fe751c7f997950e7ae..bdaab1f6fbb651d3457d23f44c036300fa87e6a5 100644
--- a/Source/Geometry/Reader/imstkVTKMeshIO.h
+++ b/Source/Geometry/Reader/imstkVTKMeshIO.h
@@ -31,6 +31,7 @@
 
 #include "imstkMeshIO.h"
 #include "imstkSurfaceMesh.h"
+#include "imstkLineMesh.h"
 #include "imstkTetrahedralMesh.h"
 #include "imstkHexahedralMesh.h"
 #include "imstkImageData.h"
@@ -81,15 +82,22 @@ protected:
     static std::shared_ptr<SurfaceMesh> readVtkPolyData(const std::string& filePath);
 
     ///
-    /// \brief Writes the given surfase mesh to given file path using the provided writer type
+    /// \brief Writes the given surface mesh to given file path using the provided writer type
     ///
     template<typename WriterType>
-    static bool writeVtkPolyData(const std::shared_ptr<SurfaceMesh> imstkMesh, const std::string& filePath);
+    static bool writeVtkPolyData(const SurfaceMesh& imstkMesh, const std::string& filePath);
+
+    ///
+    /// \brief Writes the given line mesh to given file path using the provided writer type
+    ///
+    template<typename WriterType>
+    static bool writeVtkPolyData(const LineMesh& imstkMesh, const std::string& filePath);
 
     ///
     /// \brief Writes the given volumetric mesh to given file path
     ///
-    static bool writeVtkUnstructuredGrid(const std::shared_ptr<VolumetricMesh> imstkMesh, const std::string& filePath);
+    static bool VTKMeshIO::writeVtkUnstructuredGrid(const TetrahedralMesh& tetMesh, const std::string& filePath);
+    static bool VTKMeshIO::writeVtkUnstructuredGrid(const HexahedralMesh& hMesh, const std::string& filePath);
 
     ///
     /// \brief
@@ -107,57 +115,5 @@ protected:
     /// \brief
     ///
     static std::shared_ptr<ImageData> readVtkImageDataDICOM(const std::string& filePath);
-
-    ///
-    /// \brief
-    ///
-    static std::shared_ptr<SurfaceMesh> convertVtkPolyDataToSurfaceMesh(vtkPolyData* vtkMesh);
-
-    ///
-    /// \brief Converts imstk surface mesh into a vtk polydata suitable for writing to file
-    ///
-    static vtkPolyData* convertSurfaceMeshToVtkPolyData(std::shared_ptr<SurfaceMesh> imstkMesh);
-
-    ///
-    /// \brief Converts imstk tetrahedral mesh into a vtk unstructured grid suitable for writing to file
-    ///
-    static vtkUnstructuredGrid* convertTetrahedralMeshToVtkUnstructuredGrid(std::shared_ptr<TetrahedralMesh> imstkMesh);
-
-    ///
-    /// \brief Converts imstk hexahedral mesh into a vtk unstructured grid suitable for writing to file
-    ///
-    static vtkUnstructuredGrid* convertHexahedralMeshToVtkUnstructuredGrid(std::shared_ptr<HexahedralMesh> imstkMesh);
-
-    ///
-    /// \brief
-    ///
-    static std::shared_ptr<VolumetricMesh> convertVtkUnstructuredGridToVolumetricMesh(vtkUnstructuredGrid* vtkMesh);
-
-    ///
-    /// \brief
-    ///
-    static void copyVerticesFromVtk(vtkPoints* points, StdVectorOfVec3d& vertices);
-
-    ///
-    /// \brief Copies vertices from imstk structure to VTK one
-    ///
-    static void copyVerticesToVtk(const StdVectorOfVec3d& vertices, vtkPoints* points);
-
-    ///
-    /// \brief Copies cells of the given dimension from imstk structure to VTK one
-    ///
-    template<size_t dim>
-    static void copyCellsToVtk(const std::vector<std::array<size_t, dim>>& cells, vtkCellArray* vtkCells);
-
-    ///
-    /// \brief
-    ///
-    template<size_t dim>
-    static void copyCellsFromVtk(vtkCellArray* vtkCells, std::vector<std::array<size_t, dim>>& cells);
-
-    ///
-    /// \brief
-    ///
-    static void copyPointData(vtkPointData* pointData, std::map<std::string, StdVectorOfVectorf>& dataMap);
 };
 } // imstk
diff --git a/Source/Geometry/Reader/imstkVegaMeshIO.cpp b/Source/Geometry/Reader/imstkVegaMeshIO.cpp
index ec66880b2da812f795ce817beb84f8be34393b63..f604d9977109737c03c7c77abf6e17ce752a27ed 100644
--- a/Source/Geometry/Reader/imstkVegaMeshIO.cpp
+++ b/Source/Geometry/Reader/imstkVegaMeshIO.cpp
@@ -49,7 +49,7 @@ VegaMeshIO::write(const std::shared_ptr<imstk::PointSet> imstkMesh, const std::s
 
     // extract volumetric mesh
     const auto imstkVolMesh = std::dynamic_pointer_cast<imstk::VolumetricMesh>(imstkMesh);
-    
+
     CHECK(imstkVolMesh) << "VegaMeshIO::write error: imstk::Mesh is not a volumetric mesh";
 
     switch (imstkVolMesh->getType())
@@ -57,12 +57,12 @@ VegaMeshIO::write(const std::shared_ptr<imstk::PointSet> imstkMesh, const std::s
     case Geometry::Type::TetrahedralMesh:
     case Geometry::Type::HexahedralMesh:
         auto vegaMesh = convertVolumetricMeshToVegaMesh(imstkVolMesh);
-        
+
         CHECK(vegaMesh) << "VegaMeshIO::write error: failed to convert volumetric mesh to vega mesh";
 
         const auto fileName     = const_cast<char*>(filePath.c_str());
         const int  write_status = vegaMesh->save(fileName);
-        
+
         CHECK(write_status == 0) << "VegaMeshIO::write error: failed to write .veg file";
 
         return true;
@@ -179,10 +179,10 @@ VegaMeshIO::convertVolumetricMeshToVegaMesh(const std::shared_ptr<imstk::Volumet
 
         std::shared_ptr<vega::TetMesh> vegaMesh = std::make_shared<vega::TetMesh>(int(imstkVolTetMesh->getNumVertices()), &vertices[0],
                 int(imstkVolTetMesh->getNumTetrahedra()), &elements[0], E, nu, density);
-        
+
         CHECK(vegaMesh) << "VegaMeshIO::convertVolumetricMeshToVegaMesh error: Failed to create vega mesh";
-                
-        return vegaMesh;        
+
+        return vegaMesh;
     }
     else
     {
diff --git a/Source/Geometry/imstkGeometry.h b/Source/Geometry/imstkGeometry.h
index 66608d150a4020d4bdaf8e89b811949884811d8b..1f30549765faacb4460c3ba1f5f8d83d2b934af8 100644
--- a/Source/Geometry/imstkGeometry.h
+++ b/Source/Geometry/imstkGeometry.h
@@ -205,7 +205,7 @@ protected:
     virtual void applyTranslation(const Vec3d t) = 0;
     virtual void applyRotation(const Mat3d r)    = 0;
     virtual void applyScaling(const double s)    = 0;
-    virtual void updatePostTransformData()       = 0;
+    virtual void updatePostTransformData() const const = 0;
 
     Type m_type;                 ///> Type of geometry
     std::string m_name;          ///> Unique name for each geometry
@@ -213,7 +213,7 @@ protected:
 
     bool m_dataModified      = false;
     bool m_transformModified = false;
-    bool m_transformApplied  = true;
+    mutable bool m_transformApplied = true;
 
     RigidTransform3d m_transform = RigidTransform3d::Identity(); ///> Transformation matrix
     double m_scaling = 1.0;
diff --git a/Source/Geometry/imstkGeometryUtilities.cpp b/Source/Geometry/imstkGeometryUtilities.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..af245070bf21869efd2eed26eccc031ed9bd12f9
--- /dev/null
+++ b/Source/Geometry/imstkGeometryUtilities.cpp
@@ -0,0 +1,349 @@
+/*=========================================================================
+
+   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 "imstkGeometryUtilities.h"
+#include "imstkHexahedralMesh.h"
+#include "imstkLineMesh.h"
+#include "imstkTetrahedralMesh.h"
+#include <vtkAppendPolyData.h>
+#include <vtkExtractEdges.h>
+#include <vtkLinearSubdivisionFilter.h>
+#include <vtkLoopSubdivisionFilter.h>
+#include <vtkPointData.h>
+#include <vtkSmoothPolyDataFilter.h>
+#include <vtkTriangleFilter.h>
+#include <vtkUnstructuredGrid.h>
+
+namespace imstk
+{
+std::unique_ptr<SurfaceMesh>
+GeometryUtils::convertVtkPolyDataToSurfaceMesh(vtkSmartPointer<vtkPolyData> vtkMesh)
+{
+    CHECK(vtkMesh != nullptr) << "vtkPolyData provided is not valid!";
+
+    StdVectorOfVec3d vertices;
+    copyVerticesFromVtk(vtkMesh->GetPoints(), vertices);
+
+    std::vector<SurfaceMesh::TriangleArray> triangles;
+    copyCellsFromVtk<3>(vtkMesh->GetPolys(), triangles);
+
+    auto mesh = std::make_unique<SurfaceMesh>();
+    mesh->initialize(vertices, triangles, true);
+
+    // Point Data
+    std::map<std::string, StdVectorOfVectorf> dataMap;
+    copyPointDataFromVtk(vtkMesh->GetPointData(), dataMap);
+    if (!dataMap.empty())
+    {
+        mesh->setPointDataMap(dataMap);
+    }
+
+    // Active Texture
+    if (auto pointData = vtkMesh->GetPointData())
+    {
+        if (auto tcoords = pointData->GetTCoords())
+        {
+            mesh->setDefaultTCoords(tcoords->GetName());
+        }
+    }
+
+    return mesh;
+}
+
+std::unique_ptr<LineMesh>
+GeometryUtils::convertVtkPolyDataToLineMesh(const vtkSmartPointer<vtkPolyData> vtkMesh)
+{
+    CHECK(vtkMesh != nullptr) << "vtkPolyData provided is not valid!";
+
+    StdVectorOfVec3d vertices;
+    copyVerticesFromVtk(vtkMesh->GetPoints(), vertices);
+
+    std::vector<LineMesh::LineArray> segments;
+    copyCellsFromVtk<2>(vtkMesh->GetPolys(), segments);
+
+    auto mesh = std::make_unique<LineMesh>();
+    mesh->initialize(vertices, segments);
+
+    // Point Data
+    std::map<std::string, StdVectorOfVectorf> dataMap;
+    copyPointDataFromVtk(vtkMesh->GetPointData(), dataMap);
+    if (!dataMap.empty())
+    {
+        mesh->setPointDataMap(dataMap);
+    }
+
+    return mesh;
+}
+
+vtkSmartPointer<vtkPolyData>
+GeometryUtils::convertSurfaceMeshToVtkPolyData(const SurfaceMesh& imstkMesh)
+{
+    vtkNew<vtkPoints> points;
+    copyVerticesToVtk(imstkMesh.getVertexPositions(), points.Get());
+
+    vtkNew<vtkCellArray> polys;
+    copyCellsToVtk<3>(imstkMesh.getTrianglesVertices(), polys.Get());
+
+    vtkNew<vtkPolyData> polydata;
+    polydata->SetPoints(points);
+    polydata->SetPolys(polys);
+
+    return polydata;
+}
+
+vtkSmartPointer<vtkPolyData>
+GeometryUtils::convertLineMeshToVtkPolyData(const LineMesh& imstkMesh)
+{
+    vtkNew<vtkPoints> points;
+    copyVerticesToVtk(imstkMesh.getVertexPositions(), points.Get());
+
+    vtkNew<vtkCellArray> polys;
+    copyCellsToVtk<2>(imstkMesh.getLinesVertices(), polys.Get());
+
+    vtkNew<vtkPolyData> polydata;
+    polydata->SetPoints(points);
+    polydata->SetPolys(polys);
+    return polydata;
+}
+
+vtkSmartPointer<vtkUnstructuredGrid>
+GeometryUtils::convertTetrahedralMeshToVtkUnstructuredGrid(const TetrahedralMesh& imstkMesh)
+{
+    vtkNew<vtkPoints> points;
+    copyVerticesToVtk(imstkMesh.getVertexPositions(), points.Get());
+
+    vtkNew<vtkCellArray> tetras;
+    copyCellsToVtk<4>(imstkMesh.getTetrahedraVertices(), tetras.Get());
+
+    vtkNew<vtkUnstructuredGrid> ug;
+    ug->SetPoints(points);
+    ug->SetCells(VTK_TETRA, tetras);
+    return ug;
+}
+
+vtkSmartPointer<vtkUnstructuredGrid>
+GeometryUtils::convertHexahedralMeshToVtkUnstructuredGrid(const HexahedralMesh& imstkMesh)
+{
+    vtkNew<vtkPoints> points;
+    copyVerticesToVtk(imstkMesh.getVertexPositions(), points.Get());
+
+    vtkNew<vtkCellArray> bricks;
+    copyCellsToVtk<8>(imstkMesh.getHexahedraVertices(), bricks.Get());
+
+    vtkNew<vtkUnstructuredGrid> ug;
+    ug->SetPoints(points);
+    ug->SetCells(VTK_HEXAHEDRON, bricks);
+    return ug;
+}
+
+std::unique_ptr<VolumetricMesh>
+GeometryUtils::convertVtkUnstructuredGridToVolumetricMesh(const vtkSmartPointer<vtkUnstructuredGrid> vtkMesh)
+{
+    CHECK(vtkMesh != nullptr) << "vtkUnstructuredGrid provided is not valid!";
+
+    StdVectorOfVec3d vertices;
+    copyVerticesFromVtk(vtkMesh->GetPoints(), vertices);
+
+    const int cellType = vtkMesh->GetCellType(vtkMesh->GetNumberOfCells() - 1);
+    if (cellType == VTK_TETRA)
+    {
+        std::vector<TetrahedralMesh::TetraArray> cells;
+        copyCellsFromVtk<4>(vtkMesh->GetCells(), cells);
+
+        auto mesh = std::make_unique<TetrahedralMesh>();
+        mesh->initialize(vertices, cells, false);
+        return mesh;
+    }
+    else if (cellType == VTK_HEXAHEDRON)
+    {
+        std::vector<HexahedralMesh::HexaArray> cells;
+        copyCellsFromVtk<8>(vtkMesh->GetCells(), cells);
+
+        auto mesh = std::make_unique<HexahedralMesh>();
+        mesh->initialize(vertices, cells, false);
+        return mesh;
+    }
+    else
+    {
+        LOG(FATAL) << "No support for vtkCellType = "
+                   << cellType << ".";
+        return nullptr;
+    }
+}
+
+void
+GeometryUtils::copyVerticesFromVtk(vtkPoints* const points, StdVectorOfVec3d& vertices)
+{
+    CHECK(points != nullptr) << "No valid point data provided!";
+
+    vertices.reserve(points->GetNumberOfPoints());
+    for (vtkIdType i = 0; i < points->GetNumberOfPoints(); ++i)
+    {
+        double pos[3];
+        points->GetPoint(i, pos);
+        vertices.emplace_back(pos[0], pos[1], pos[2]);
+    }
+}
+
+void
+GeometryUtils::copyVerticesToVtk(const StdVectorOfVec3d& vertices, vtkPoints* points)
+{
+    CHECK(points != nullptr) << "No valid vtkPoints pointer provided!";
+
+    points->SetNumberOfPoints(vertices.size());
+    for (size_t i = 0; i < vertices.size(); i++)
+    {
+        points->SetPoint(i, vertices[i][0], vertices[i][1], vertices[i][2]);
+    }
+}
+
+template<size_t dim>
+void
+GeometryUtils::copyCellsToVtk(const std::vector<std::array<size_t, dim>>& cells, vtkCellArray* vtkCells)
+{
+    CHECK(vtkCells != nullptr) << "No valid vtkCellArray provided!";
+
+    for (size_t i = 0; i < cells.size(); i++)
+    {
+        vtkCells->InsertNextCell(dim);
+        for (size_t k = 0; k < dim; k++)
+        {
+            vtkCells->InsertCellPoint(cells[i][k]);
+        }
+    }
+}
+
+template<size_t dim>
+void
+GeometryUtils::copyCellsFromVtk(vtkCellArray* vtkCells, std::vector<std::array<size_t, dim>>& cells)
+{
+    CHECK(vtkCells != nullptr) << "No cells found!";
+
+    vtkNew<vtkIdList>       vtkCell;
+    std::array<size_t, dim> cell;
+    cells.reserve(vtkCells->GetNumberOfCells());
+    vtkCells->InitTraversal();
+    while (vtkCells->GetNextCell(vtkCell))
+    {
+        if (vtkCell->GetNumberOfIds() != dim)
+        {
+            continue;
+        }
+        for (size_t i = 0; i < dim; ++i)
+        {
+            cell[i] = vtkCell->GetId(i);
+        }
+        cells.emplace_back(cell);
+    }
+}
+
+void
+GeometryUtils::copyPointDataFromVtk(vtkPointData* const pointData, std::map<std::string, StdVectorOfVectorf>& dataMap)
+{
+    CHECK(pointData != nullptr) << "No point data provided!";
+
+    for (int i = 0; i < pointData->GetNumberOfArrays(); ++i)
+    {
+        vtkDataArray*      array       = pointData->GetArray(i);
+        std::string        name        = array->GetName();
+        int                nbrOfComp   = array->GetNumberOfComponents();
+        vtkIdType          nbrOfTuples = array->GetNumberOfTuples();
+        StdVectorOfVectorf data;
+        for (vtkIdType j = 0; j < nbrOfTuples; ++j)
+        {
+            double* tupleData = new double[nbrOfComp];
+            array->GetTuple(j, tupleData);
+            Vectorf tuple(nbrOfComp);
+            for (int k = 0; k < nbrOfComp; k++)
+            {
+                tuple[k] = tupleData[k];
+            }
+            data.push_back(tuple);
+        }
+        dataMap[name] = data;
+    }
+}
+
+std::unique_ptr<SurfaceMesh>
+GeometryUtils::combineSurfaceMesh(const SurfaceMesh& surfaceMesh1, const SurfaceMesh& surfaceMesh2)
+{
+    vtkNew<vtkAppendPolyData> filter;
+    filter->AddInputData(convertSurfaceMeshToVtkPolyData(surfaceMesh1));
+    filter->AddInputData(convertSurfaceMeshToVtkPolyData(surfaceMesh2));
+    filter->Update();
+
+    return convertVtkPolyDataToSurfaceMesh(filter->GetOutput());
+}
+
+std::unique_ptr<LineMesh>
+GeometryUtils::surfaceMeshToLineMesh(const SurfaceMesh& surfaceMesh)
+{
+    vtkNew<vtkExtractEdges> filter;
+    filter->SetInputData(convertSurfaceMeshToVtkPolyData(surfaceMesh));
+    filter->Update();
+
+    vtkNew<vtkTriangleFilter> triangleFilter;
+    triangleFilter->SetInputData(filter->GetOutput());
+    triangleFilter->Update();
+
+    return convertVtkPolyDataToLineMesh(triangleFilter->GetOutput());
+}
+
+std::unique_ptr<SurfaceMesh>
+GeometryUtils::smoothSurfaceMesh(const SurfaceMesh&          surfaceMesh,
+                                 const smoothPolydataConfig& c)
+{
+    vtkNew<vtkSmoothPolyDataFilter> filter;
+    filter->SetInputData(convertSurfaceMeshToVtkPolyData(surfaceMesh));
+    filter->SetNumberOfIterations(c.numberOfIterations);
+    filter->SetRelaxationFactor(c.relaxationFactor);
+    filter->SetConvergence(c.convergence);
+    filter->SetFeatureAngle(c.featureAngle);
+    filter->SetEdgeAngle(c.edgeAngle);
+    filter->SetFeatureEdgeSmoothing(c.featureEdgeSmoothing);
+    filter->SetBoundarySmoothing(c.boundarySmoothing);
+    filter->Update();
+
+    return convertVtkPolyDataToSurfaceMesh(filter->GetOutput());
+}
+
+std::unique_ptr<SurfaceMesh>
+GeometryUtils::linearSubdivideSurfaceMesh(const SurfaceMesh& surfaceMesh, const int numSubdivisions)
+{
+    vtkNew<vtkLinearSubdivisionFilter> filter;
+    filter->SetInputData(convertSurfaceMeshToVtkPolyData(surfaceMesh));
+    filter->SetNumberOfSubdivisions(numSubdivisions);
+    filter->Update();
+
+    return convertVtkPolyDataToSurfaceMesh(filter->GetOutput());
+}
+
+std::unique_ptr<SurfaceMesh>
+GeometryUtils::loopSubdivideSurfaceMesh(const SurfaceMesh& surfaceMesh, const int numSubdivisions)
+{
+    vtkNew<vtkLoopSubdivisionFilter> filter;
+    filter->SetInputData(convertSurfaceMeshToVtkPolyData(surfaceMesh));
+    filter->SetNumberOfSubdivisions(numSubdivisions);
+    filter->Update();
+
+    return convertVtkPolyDataToSurfaceMesh(filter->GetOutput());
+}
+}
\ No newline at end of file
diff --git a/Source/Geometry/imstkGeometryUtilities.h b/Source/Geometry/imstkGeometryUtilities.h
new file mode 100644
index 0000000000000000000000000000000000000000..55db04a8b6cd86129df266de546442a0ecfb2832
--- /dev/null
+++ b/Source/Geometry/imstkGeometryUtilities.h
@@ -0,0 +1,157 @@
+/*=========================================================================
+
+   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 "imstkMath.h"
+#include <memory>
+#include <vtkSmartPointer.h>
+
+class vtkCellArray;
+class vtkPolyData;
+class vtkPointData;
+class vtkPoints;
+class vtkUnstructuredGrid;
+
+namespace imstk
+{
+class HexahedralMesh;
+class LineMesh;
+class SurfaceMesh;
+class TetrahedralMesh;
+class VolumetricMesh;
+
+namespace GeometryUtils
+{
+///
+/// \brief Converts vtk polydata into a imstk surface mesh
+///
+std::unique_ptr<SurfaceMesh> convertVtkPolyDataToSurfaceMesh(const vtkSmartPointer<vtkPolyData> vtkMesh);
+
+///
+/// \brief Converts vtk polydata into a imstk surface mesh
+///
+std::unique_ptr<LineMesh> convertVtkPolyDataToLineMesh(const vtkSmartPointer<vtkPolyData> vtkMesh);
+
+///
+/// \brief Get imstk volumetric mesh given vtkUnstructuredGrid as input
+///
+std::unique_ptr<VolumetricMesh> convertVtkUnstructuredGridToVolumetricMesh(const vtkSmartPointer<vtkUnstructuredGrid> vtkMesh);
+
+///
+/// \brief Converts imstk surface mesh into a vtk polydata
+///
+vtkSmartPointer<vtkPolyData> convertSurfaceMeshToVtkPolyData(const SurfaceMesh& imstkMesh);
+
+///
+/// \brief Converts imstk line mesh into a vtk polydata
+///
+vtkSmartPointer<vtkPolyData> convertLineMeshToVtkPolyData(const LineMesh& imstkMesh);
+
+///
+/// \brief Converts imstk tetrahedral mesh into a vtk unstructured grid
+///
+vtkSmartPointer<vtkUnstructuredGrid> convertTetrahedralMeshToVtkUnstructuredGrid(const TetrahedralMesh& imstkMesh);
+
+///
+/// \brief Converts imstk hexahedral mesh into a vtk unstructured grid
+///
+vtkSmartPointer<vtkUnstructuredGrid> convertHexahedralMeshToVtkUnstructuredGrid(const HexahedralMesh& imstkMesh);
+
+///
+/// \brief Copy from vtk points to a imstk vertices array (StdVectorOfVec3d)
+///
+void copyVerticesFromVtk(vtkPoints* const points, StdVectorOfVec3d& vertices);
+
+///
+/// \brief Copies vertices from imstk structure to VTK one
+///
+void copyVerticesToVtk(const StdVectorOfVec3d& vertices, vtkPoints* points);
+
+///
+/// \brief Copies cells of the given dimension from imstk structure to VTK one
+///
+template<size_t dim>
+void copyCellsToVtk(const std::vector<std::array<size_t, dim>>& cells, vtkCellArray* vtkCells);
+
+///
+/// \brief
+///
+template<size_t dim>
+void copyCellsFromVtk(vtkCellArray* vtkCells, std::vector<std::array<size_t, dim>>& cells);
+
+///
+/// \brief
+///
+void copyPointDataFromVtk(vtkPointData* const pointData, std::map<std::string, StdVectorOfVectorf>& dataMap);
+
+///
+/// \brief Combines two input surface meshes
+/// Refer <a href="https://vtk.org/doc/nightly/html/classvtkAppendPolyData.html#details">vtkAppendPolyData</a> class
+/// for more details
+///
+///
+std::unique_ptr<SurfaceMesh> combineSurfaceMesh(const SurfaceMesh& surfaceMesh1, const SurfaceMesh& surfaceMesh2);
+
+///
+/// \brief Converts an imstk SurfaceMesh to a LineMesh, removing duplicate edges. Cell indices not preserved
+/// Refer <a href="https://vtk.org/doc/nightly/html/classvtkExtractEdges.html#details">vtkExtractEdges</a> class
+/// for more details
+///
+std::unique_ptr<LineMesh> surfaceMeshToLineMesh(const SurfaceMesh& surfaceMesh);
+
+///
+/// \brief Config for smooth polydata filter
+///
+struct smoothPolydataConfig
+{
+    int numberOfIterations  = 20;
+    double relaxationFactor = 0.01;
+    double convergence      = 0.0;
+    double featureAngle     = 45.0;
+    double edgeAngle = 15.0;
+    bool featureEdgeSmoothing = false;
+    bool boundarySmoothing    = true;
+};
+
+///
+/// \brief Smooths a SurfaceMesh using laplacian smoothening
+/// Refer <a href="https://vtk.org/doc/nightly/html/classvtkSmoothPolyDataFilter.html#details">vtkSmoothPolyDataFilter</a>
+/// for more details
+///
+std::unique_ptr<SurfaceMesh> smoothSurfaceMesh(const SurfaceMesh&          surfaceMesh,
+                                               const smoothPolydataConfig& c);
+
+///
+/// \brief Sub-divdes a SurfaceMesh using linear subdivision
+/// Refer <a href="https://vtk.org/doc/nightly/html/classvtkLinearSubdivisionFilter.html#details">vtk linear subdivision</a>
+/// for more details
+///
+std::unique_ptr<SurfaceMesh> linearSubdivideSurfaceMesh(const SurfaceMesh& surfaceMesh, const int numSubdivisions = 1);
+
+///
+/// \brief Sub-divides an input imstk SurfaceMesh using loop subdivision algorithm
+/// Refer <a href="https://vtk.org/doc/nightly/html/classvtkLoopSubdivisionFilter.html#details">vtk loop subdivision</a>
+/// for more details
+///
+std::unique_ptr<SurfaceMesh> loopSubdivideSurfaceMesh(const SurfaceMesh& surfaceMesh, const int numSubdivisions = 1);
+};
+}
diff --git a/Source/Scene/SceneElements/Objects/imstkPbdObject.cpp b/Source/Scene/SceneElements/Objects/imstkPbdObject.cpp
index 10407c521b5b3ae3f68ab7b0d39efefceb31ceeb..33f7b76aea628ed195f5ee43d8004aea311250d6 100644
--- a/Source/Scene/SceneElements/Objects/imstkPbdObject.cpp
+++ b/Source/Scene/SceneElements/Objects/imstkPbdObject.cpp
@@ -33,7 +33,7 @@ PbdObject::initialize()
     m_pbdModel = std::dynamic_pointer_cast<PbdModel>(m_dynamicalModel);
 
     CHECK(m_pbdModel) << "Dynamics pointer cast failure in PbdObject::initialize()";
-    
+
     return DynamicObject::initialize();
 }
 
diff --git a/Source/Scene/imstkScene.cpp b/Source/Scene/imstkScene.cpp
index 363a225e1ef1e033afa62759bbec16b5fa42ed66..e3bd72fbeff62992d82c2a18171803987e5cb04d 100644
--- a/Source/Scene/imstkScene.cpp
+++ b/Source/Scene/imstkScene.cpp
@@ -137,8 +137,8 @@ Scene::getSceneObjectControllers() const
 std::shared_ptr<SceneObject>
 Scene::getSceneObject(const std::string& sceneObjectName) const
 {
-    CHECK(this->isObjectRegistered(sceneObjectName)) 
-        << "No scene object named '"  << sceneObjectName
+    CHECK(this->isObjectRegistered(sceneObjectName))
+        << "No scene object named '" << sceneObjectName
         << "' was registered in this scene.";
 
     return m_sceneObjectsMap.at(sceneObjectName);