diff --git a/momentum/character/character_utility.cpp b/momentum/character/character_utility.cpp
index 08355608c3..f83ee9f1e0 100644
--- a/momentum/character/character_utility.cpp
+++ b/momentum/character/character_utility.cpp
@@ -845,6 +845,53 @@ std::vector<Eigen::Vector3i> remapFaces(
   return remappedFaces;
 }
 
+template <typename T>
+void remapPolyFaces(
+    MeshT<T>& newMesh,
+    const MeshT<T>& mesh,
+    const std::vector<bool>& activeVertices,
+    const std::vector<size_t>& reverseVertexMapping,
+    const std::vector<size_t>& reverseTexcoordMapping) {
+  if (mesh.polyFaces.empty()) {
+    return;
+  }
+
+  newMesh.polyFaceSizes.reserve(mesh.polyFaceSizes.size());
+  newMesh.polyFaces.reserve(mesh.polyFaces.size());
+  if (!mesh.polyTexcoordFaces.empty()) {
+    newMesh.polyTexcoordFaces.reserve(mesh.polyTexcoordFaces.size());
+  }
+
+  uint32_t offset = 0;
+  for (const auto polySize : mesh.polyFaceSizes) {
+    // Check if all vertices of this polygon are active
+    bool allActive = true;
+    for (uint32_t i = 0; i < polySize; ++i) {
+      const auto vtx = mesh.polyFaces[offset + i];
+      if (vtx >= activeVertices.size() || !activeVertices[vtx]) {
+        allActive = false;
+        break;
+      }
+    }
+
+    if (allActive) {
+      newMesh.polyFaceSizes.push_back(polySize);
+      for (uint32_t i = 0; i < polySize; ++i) {
+        newMesh.polyFaces.push_back(
+            static_cast<uint32_t>(reverseVertexMapping[mesh.polyFaces[offset + i]]));
+      }
+      if (!mesh.polyTexcoordFaces.empty()) {
+        for (uint32_t i = 0; i < polySize; ++i) {
+          newMesh.polyTexcoordFaces.push_back(
+              static_cast<uint32_t>(reverseTexcoordMapping[mesh.polyTexcoordFaces[offset + i]]));
+        }
+      }
+    }
+
+    offset += polySize;
+  }
+}
+
 std::vector<std::vector<int32_t>> remapLines(
     const std::vector<std::vector<int32_t>>& lines,
     const std::vector<size_t>& mapping) {
@@ -923,6 +970,53 @@ std::vector<bool> facesToVertices(
   return activeVertices;
 }
 
+std::vector<bool> verticesToPolys(
+    const std::vector<uint32_t>& polyFaces,
+    const std::vector<uint32_t>& polyFaceSizes,
+    const std::vector<bool>& activeVertices) {
+  if (polyFaceSizes.empty()) {
+    return {};
+  }
+  std::vector<bool> activePolys(polyFaceSizes.size(), false);
+  uint32_t offset = 0;
+  for (size_t polyIdx = 0; polyIdx < polyFaceSizes.size(); ++polyIdx) {
+    const auto polySize = polyFaceSizes[polyIdx];
+    bool allActive = true;
+    for (uint32_t i = 0; i < polySize; ++i) {
+      const auto vtx = polyFaces[offset + i];
+      if (vtx >= activeVertices.size() || !activeVertices[vtx]) {
+        allActive = false;
+        break;
+      }
+    }
+    activePolys[polyIdx] = allActive;
+    offset += polySize;
+  }
+  return activePolys;
+}
+
+std::vector<bool> polysToVertices(
+    const std::vector<uint32_t>& polyFaces,
+    const std::vector<uint32_t>& polyFaceSizes,
+    const std::vector<bool>& activePolys,
+    size_t vertexCount) {
+  std::vector<bool> activeVertices(vertexCount, false);
+  uint32_t offset = 0;
+  for (size_t polyIdx = 0; polyIdx < polyFaceSizes.size(); ++polyIdx) {
+    const auto polySize = polyFaceSizes[polyIdx];
+    if (activePolys[polyIdx]) {
+      for (uint32_t i = 0; i < polySize; ++i) {
+        const auto vtx = polyFaces[offset + i];
+        if (vtx < vertexCount) {
+          activeVertices[vtx] = true;
+        }
+      }
+    }
+    offset += polySize;
+  }
+  return activeVertices;
+}
+
 // Internal mesh reduction function used by both reduceMeshComponents and
 // the public reduceMeshByFaces/reduceMeshByVertices for Mesh.
 template <typename T>
@@ -954,6 +1048,7 @@ MeshT<T> reduceMeshInternal(
     newMesh.lines = remapLines(mesh.lines, reverseVertexMapping);
   }
 
+  std::vector<size_t> reverseTextureVertexMapping;
   if (!mesh.texcoord_faces.empty() && !mesh.texcoords.empty()) {
     std::vector<bool> activeTextureTriangles = activeFaces;
     activeTextureTriangles.resize(mesh.texcoord_faces.size(), false);
@@ -961,8 +1056,8 @@ MeshT<T> reduceMeshInternal(
     const std::vector<bool> activeTextureVertices =
         facesToVertices(mesh.texcoord_faces, activeTextureTriangles, mesh.texcoords.size());
 
-    auto [forwardTextureVertexMapping, reverseTextureVertexMapping] =
-        createIndexMapping(activeTextureVertices);
+    auto [forwardTextureVertexMapping, revTexMapping] = createIndexMapping(activeTextureVertices);
+    reverseTextureVertexMapping = std::move(revTexMapping);
 
     newMesh.texcoords = selectVertices(mesh.texcoords, forwardTextureVertexMapping);
     newMesh.texcoord_faces =
@@ -970,6 +1065,9 @@ MeshT<T> reduceMeshInternal(
     newMesh.texcoord_lines = remapLines(mesh.texcoord_lines, reverseTextureVertexMapping);
   }
 
+  // Remap polygon face data
+  remapPolyFaces(newMesh, mesh, activeVertices, reverseVertexMapping, reverseTextureVertexMapping);
+
   return newMesh;
 }
 
@@ -1185,6 +1283,61 @@ std::vector<bool> facesToVertices(const MeshT<T>& mesh, const std::vector<bool>&
   return facesToVertices(mesh.faces, activeFaces, mesh.vertices.size());
 }
 
+template <typename T>
+std::vector<bool> verticesToPolys(const MeshT<T>& mesh, const std::vector<bool>& activeVertices) {
+  MT_CHECK(
+      activeVertices.size() == mesh.vertices.size(),
+      "Active vertices size ({}) does not match mesh vertex count ({})",
+      activeVertices.size(),
+      mesh.vertices.size());
+
+  return verticesToPolys(mesh.polyFaces, mesh.polyFaceSizes, activeVertices);
+}
+
+template <typename T>
+std::vector<bool> polysToVertices(const MeshT<T>& mesh, const std::vector<bool>& activePolys) {
+  MT_CHECK(
+      activePolys.size() == mesh.polyFaceSizes.size(),
+      "Active polys size ({}) does not match polygon count ({})",
+      activePolys.size(),
+      mesh.polyFaceSizes.size());
+
+  return polysToVertices(mesh.polyFaces, mesh.polyFaceSizes, activePolys, mesh.vertices.size());
+}
+
+template <typename T>
+MeshT<T> reduceMeshByPolys(const MeshT<T>& mesh, const std::vector<bool>& activePolys) {
+  MT_CHECK(
+      activePolys.size() == mesh.polyFaceSizes.size(),
+      "Active polys size ({}) does not match polygon count ({})",
+      activePolys.size(),
+      mesh.polyFaceSizes.size());
+
+  // Convert polygon selection to vertex selection, then to face selection
+  const auto activeVertices = polysToVertices(mesh, activePolys);
+  const auto activeFaces = verticesToFaces(mesh, activeVertices);
+
+  return reduceMeshInternal(mesh, activeVertices, activeFaces);
+}
+
+template <typename T>
+CharacterT<T> reduceMeshByPolys(
+    const CharacterT<T>& character,
+    const std::vector<bool>& activePolys) {
+  MT_CHECK(character.mesh, "Cannot reduce mesh: character has no mesh");
+  MT_CHECK(
+      activePolys.size() == character.mesh->polyFaceSizes.size(),
+      "Active polys size ({}) does not match polygon count ({})",
+      activePolys.size(),
+      character.mesh->polyFaceSizes.size());
+
+  // Convert polygon selection to vertex selection, then to face selection
+  const auto activeVertices = polysToVertices(*character.mesh, activePolys);
+  const auto activeFaces = verticesToFaces(*character.mesh, activeVertices);
+
+  return reduceMeshComponents(character, activeVertices, activeFaces);
+}
+
 // Explicit instantiations for commonly used types
 template CharacterT<float> reduceMeshByVertices<float>(
     const CharacterT<float>& character,
@@ -1234,4 +1387,36 @@ template std::vector<bool> facesToVertices<double>(
     const MeshT<double>& mesh,
     const std::vector<bool>& activeFaces);
 
+template std::vector<bool> verticesToPolys<float>(
+    const MeshT<float>& mesh,
+    const std::vector<bool>& activeVertices);
+
+template std::vector<bool> verticesToPolys<double>(
+    const MeshT<double>& mesh,
+    const std::vector<bool>& activeVertices);
+
+template std::vector<bool> polysToVertices<float>(
+    const MeshT<float>& mesh,
+    const std::vector<bool>& activePolys);
+
+template std::vector<bool> polysToVertices<double>(
+    const MeshT<double>& mesh,
+    const std::vector<bool>& activePolys);
+
+template MeshT<float> reduceMeshByPolys<float>(
+    const MeshT<float>& mesh,
+    const std::vector<bool>& activePolys);
+
+template MeshT<double> reduceMeshByPolys<double>(
+    const MeshT<double>& mesh,
+    const std::vector<bool>& activePolys);
+
+template CharacterT<float> reduceMeshByPolys<float>(
+    const CharacterT<float>& character,
+    const std::vector<bool>& activePolys);
+
+template CharacterT<double> reduceMeshByPolys<double>(
+    const CharacterT<double>& character,
+    const std::vector<bool>& activePolys);
+
 } // namespace momentum
diff --git a/momentum/character/character_utility.h b/momentum/character/character_utility.h
index 505cfc020b..4618e85531 100644
--- a/momentum/character/character_utility.h
+++ b/momentum/character/character_utility.h
@@ -145,7 +145,7 @@ template <typename T>
 
 /// Converts face selection to vertex selection
 ///
-/// @param[in] character Character containing the mesh
+/// @param[in] mesh The mesh containing the faces
 /// @param[in] activeFaces Boolean vector indicating which faces are active
 /// @return Boolean vector indicating which vertices are used by active faces
 template <typename T>
@@ -153,4 +153,52 @@ template <typename T>
     const MeshT<T>& mesh,
     const std::vector<bool>& activeFaces);
 
+/// Converts vertex selection to polygon selection
+///
+/// A polygon is active if all its vertices are active.
+///
+/// @param[in] mesh The mesh containing the polygon data
+/// @param[in] activeVertices Boolean vector indicating which vertices are active
+/// @return Boolean vector indicating which polygons only contain active vertices
+template <typename T>
+[[nodiscard]] std::vector<bool> verticesToPolys(
+    const MeshT<T>& mesh,
+    const std::vector<bool>& activeVertices);
+
+/// Converts polygon selection to vertex selection
+///
+/// @param[in] mesh The mesh containing the polygon data
+/// @param[in] activePolys Boolean vector indicating which polygons are active
+/// @return Boolean vector indicating which vertices are used by active polygons
+template <typename T>
+[[nodiscard]] std::vector<bool> polysToVertices(
+    const MeshT<T>& mesh,
+    const std::vector<bool>& activePolys);
+
+/// Reduces a standalone mesh to only include the specified polygons and associated vertices.
+///
+/// Converts polygon selection to vertex selection and face selection, then reduces the mesh.
+/// Both triangulated faces and polygon data are filtered and remapped.
+///
+/// @param[in] mesh The mesh to reduce
+/// @param[in] activePolys Boolean vector indicating which polygons to keep
+/// @return A new mesh containing only the specified polygons and their referenced vertices
+template <typename T>
+[[nodiscard]] MeshT<T> reduceMeshByPolys(
+    const MeshT<T>& mesh,
+    const std::vector<bool>& activePolys);
+
+/// Reduces the mesh to only include the specified polygons and associated vertices
+///
+/// Converts polygon selection to vertex selection and face selection, then reduces the mesh.
+/// Both triangulated faces and polygon data are filtered and remapped.
+///
+/// @param[in] character Character to be reduced
+/// @param[in] activePolys Boolean vector indicating which polygons to keep
+/// @return A new character with mesh reduced to the specified polygons
+template <typename T>
+[[nodiscard]] CharacterT<T> reduceMeshByPolys(
+    const CharacterT<T>& character,
+    const std::vector<bool>& activePolys);
+
 } // namespace momentum
diff --git a/momentum/io/fbx/fbx_io.cpp b/momentum/io/fbx/fbx_io.cpp
index 5d76690e9c..2516404ac0 100644
--- a/momentum/io/fbx/fbx_io.cpp
+++ b/momentum/io/fbx/fbx_io.cpp
@@ -561,6 +561,56 @@ void addMetaData(::fbxsdk::FbxNode* skeletonRootNode, const Character& character
   }
 }
 
+void writePolygonsToFbxMesh(const Mesh& mesh, ::fbxsdk::FbxMesh* lMesh) {
+  if (!mesh.polyFaces.empty() && !mesh.polyFaceSizes.empty()) {
+    // Write original polygon topology
+    uint32_t offset = 0;
+    for (const auto polySize : mesh.polyFaceSizes) {
+      lMesh->BeginPolygon();
+      for (uint32_t i = 0; i < polySize; ++i) {
+        lMesh->AddPolygon(mesh.polyFaces[offset + i]);
+      }
+      lMesh->EndPolygon();
+      offset += polySize;
+    }
+  } else {
+    // Fall back to triangulated faces
+    for (const auto& face : mesh.faces) {
+      lMesh->BeginPolygon();
+      for (int i = 0; i < 3; i++) {
+        lMesh->AddPolygon(face[i]);
+      }
+      lMesh->EndPolygon();
+    }
+  }
+}
+
+void writeTextureUVIndicesToFbxMesh(
+    const Mesh& mesh,
+    ::fbxsdk::FbxMesh* lMesh,
+    fbxsdk::FbxLayerElement::EType uvType) {
+  if (!mesh.polyFaces.empty() && !mesh.polyFaceSizes.empty()) {
+    uint32_t offset = 0;
+    int polyIdx = 0;
+    for (const auto polySize : mesh.polyFaceSizes) {
+      for (uint32_t i = 0; i < polySize; ++i) {
+        const uint32_t texIdx = mesh.polyTexcoordFaces.empty() ? mesh.polyFaces[offset + i]
+                                                               : mesh.polyTexcoordFaces[offset + i];
+        lMesh->SetTextureUVIndex(polyIdx, i, texIdx, uvType);
+      }
+      offset += polySize;
+      polyIdx++;
+    }
+  } else {
+    for (int faceIdx = 0; faceIdx < static_cast<int>(mesh.texcoord_faces.size()); ++faceIdx) {
+      const auto& texcoords = mesh.texcoord_faces[faceIdx];
+      lMesh->SetTextureUVIndex(faceIdx, 0, texcoords[0], uvType);
+      lMesh->SetTextureUVIndex(faceIdx, 1, texcoords[1], uvType);
+      lMesh->SetTextureUVIndex(faceIdx, 2, texcoords[2], uvType);
+    }
+  }
+}
+
 void saveFbxCommon(
     const filesystem::path& filename,
     const Character& character,
@@ -687,7 +737,6 @@ void saveFbxCommon(
   if (saveMesh && character.mesh != nullptr) {
     // Add the mesh
     const int numVertices = character.mesh.get()->vertices.size();
-    const int numFaces = character.mesh.get()->faces.size();
     ::fbxsdk::FbxNode* meshNode = ::fbxsdk::FbxNode::Create(scene, "body_mesh");
     ::fbxsdk::FbxMesh* lMesh = ::fbxsdk::FbxMesh::Create(scene, "mesh");
     lMesh->SetControlPointCount(numVertices);
@@ -704,14 +753,7 @@ void saveFbxCommon(
       lMesh->SetControlPointAt(point, normal, i);
     }
     // Add polygons to lMesh
-    for (int iFace = 0; iFace < numFaces; iFace++) {
-      lMesh->BeginPolygon();
-      for (int i = 0; i < 3; i++) { // We have tris for models. This could be extended for
-                                    // supporting Quads or npoly if needed.
-        lMesh->AddPolygon(character.mesh.get()->faces[iFace][i]);
-      }
-      lMesh->EndPolygon();
-    }
+    writePolygonsToFbxMesh(*character.mesh, lMesh);
     lMesh->BuildMeshEdgeArray();
     meshNode->SetNodeAttribute(lMesh);
 
@@ -733,13 +775,8 @@ void saveFbxCommon(
         lMesh->AddTextureUV(::fbxsdk::FbxVector2(texcoords[0], 1.0f - texcoords[1]), uvType);
       }
 
-      // Set UV indices for each face. We only have triangles.
-      int faceCount = 0;
-      for (const auto& texcoords : character.mesh->texcoord_faces) {
-        lMesh->SetTextureUVIndex(faceCount, 0, texcoords[0], uvType);
-        lMesh->SetTextureUVIndex(faceCount, 1, texcoords[1], uvType);
-        lMesh->SetTextureUVIndex(faceCount++, 2, texcoords[2], uvType);
-      }
+      // Set UV indices for each face.
+      writeTextureUVIndicesToFbxMesh(*character.mesh, lMesh, uvType);
     }
 
     // ---------------------------------------------
diff --git a/momentum/io/fbx/openfbx_loader.cpp b/momentum/io/fbx/openfbx_loader.cpp
index e34b55d66a..7755012dc4 100644
--- a/momentum/io/fbx/openfbx_loader.cpp
+++ b/momentum/io/fbx/openfbx_loader.cpp
@@ -501,6 +501,28 @@ parseSkeleton(const ofbx::Object* sceneRoot, const std::string& skelRoot, Permis
   return {skeleton, jointFbxNodes, locators, collision};
 }
 
+void populatePolyFaceData(
+    Mesh& mesh,
+    const PolygonData& vertices,
+    size_t vertexOffset,
+    size_t textureCoordOffset) {
+  mesh.polyFaces.reserve(mesh.polyFaces.size() + vertices.indices.size());
+  for (const auto& idx : vertices.indices) {
+    mesh.polyFaces.push_back(idx + static_cast<uint32_t>(vertexOffset));
+  }
+  mesh.polyFaceSizes.reserve(
+      mesh.polyFaceSizes.size() + (vertices.offsets.empty() ? 0 : vertices.offsets.size() - 1));
+  for (size_t i = 0; i + 1 < vertices.offsets.size(); ++i) {
+    mesh.polyFaceSizes.push_back(vertices.offsets[i + 1] - vertices.offsets[i]);
+  }
+  if (!vertices.textureIndices.empty()) {
+    mesh.polyTexcoordFaces.reserve(mesh.polyTexcoordFaces.size() + vertices.textureIndices.size());
+    for (const auto& idx : vertices.textureIndices) {
+      mesh.polyTexcoordFaces.push_back(idx + static_cast<uint32_t>(textureCoordOffset));
+    }
+  }
+}
+
 void parseSkinnedModel(
     const ofbx::Mesh* meshRoot,
     const std::vector<const ofbx::Object*>& boneFbxNodes,
@@ -659,6 +681,9 @@ void parseSkinnedModel(
     mesh.texcoord_faces.emplace_back(t + Eigen::Vector3i::Constant(textureCoordOffset));
   }
 
+  // Populate polygon face data (preserves original topology)
+  populatePolyFaceData(mesh, vertices, vertexOffset, textureCoordOffset);
+
   const auto* blendshapes = geometry->getBlendShape();
   if (loadBlendShapes == LoadBlendShapes::Yes && blendshapes) {
     if (!blendShape) {
diff --git a/momentum/math/mesh.cpp b/momentum/math/mesh.cpp
index a94b83b48f..14a33bb6c4 100644
--- a/momentum/math/mesh.cpp
+++ b/momentum/math/mesh.cpp
@@ -69,6 +69,9 @@ MeshT<T2> MeshT<T>::cast() const {
   result.texcoords = this->texcoords;
   result.texcoord_faces = this->texcoord_faces;
   result.texcoord_lines = this->texcoord_lines;
+  result.polyFaces = this->polyFaces;
+  result.polyFaceSizes = this->polyFaceSizes;
+  result.polyTexcoordFaces = this->polyTexcoordFaces;
   return result;
 }
 
@@ -83,6 +86,9 @@ void MeshT<T>::reset() {
   texcoords.clear();
   texcoord_faces.clear();
   texcoord_lines.clear();
+  polyFaces.clear();
+  polyFaceSizes.clear();
+  polyTexcoordFaces.clear();
 }
 
 template MeshT<float> MeshT<float>::cast<float>() const;
diff --git a/momentum/math/mesh.h b/momentum/math/mesh.h
index 01c2e54076..9eb7da9d3f 100644
--- a/momentum/math/mesh.h
+++ b/momentum/math/mesh.h
@@ -56,6 +56,26 @@ struct MeshT {
   /// Maps each line to its corresponding texture coordinates
   std::vector<std::vector<int32_t>> texcoord_lines;
 
+  /// Packed polygon face vertex indices (all polygons concatenated back-to-back).
+  ///
+  /// For example, a quad followed by a triangle would be: [v0, v1, v2, v3, v4, v5, v6].
+  /// This is an optional, redundant representation — the triangulated @ref faces field is always
+  /// required and used by all downstream code. Polygon data preserves the original topology (quads,
+  /// n-gons) from source files like FBX.
+  std::vector<uint32_t> polyFaces;
+
+  /// Number of vertices in each polygon face.
+  ///
+  /// For example, [4, 3] means the first polygon is a quad, the second is a triangle.
+  /// May be empty if polygon data is not available; see @ref polyFaces.
+  std::vector<uint32_t> polyFaceSizes;
+
+  /// Packed polygon face texture coordinate indices (same layout as polyFaces).
+  ///
+  /// Uses the same polyFaceSizes array — polygon sizes are identical for geometry and texcoords.
+  /// May be empty if texture coordinates are not available; see @ref polyFaces.
+  std::vector<uint32_t> polyTexcoordFaces;
+
   /// Compute vertex normals by averaging connected face normals.
   ///
   /// This method calculates normals for each vertex by averaging the normals of all
diff --git a/momentum/test/character/character_helpers_gtest.cpp b/momentum/test/character/character_helpers_gtest.cpp
index b726680f2a..7147603cb6 100644
--- a/momentum/test/character/character_helpers_gtest.cpp
+++ b/momentum/test/character/character_helpers_gtest.cpp
@@ -100,6 +100,12 @@ void compareMeshes(const Mesh_u& refMesh, const Mesh_u& mesh) {
       refMesh->confidence, testing::Pointwise(testing::DoubleNear(0.0001), mesh->confidence));
   EXPECT_THAT(
       refMesh->texcoord_faces, testing::Pointwise(IntExactPointwise(), mesh->texcoord_faces));
+  // Only compare polygon data when both meshes have it (e.g. glTF doesn't populate poly fields)
+  if (!refMesh->polyFaces.empty() && !mesh->polyFaces.empty()) {
+    EXPECT_EQ(refMesh->polyFaces, mesh->polyFaces);
+    EXPECT_EQ(refMesh->polyFaceSizes, mesh->polyFaceSizes);
+    EXPECT_EQ(refMesh->polyTexcoordFaces, mesh->polyTexcoordFaces);
+  }
 }
 
 void compareBlendShapes(const BlendShape_const_p& refShapes, const BlendShape_const_p& shapes) {
diff --git a/momentum/test/character/character_test.cpp b/momentum/test/character/character_test.cpp
index b3d5ccab23..a36a5161c3 100644
--- a/momentum/test/character/character_test.cpp
+++ b/momentum/test/character/character_test.cpp
@@ -1568,3 +1568,157 @@ TYPED_TEST(CharacterTest, MeshReductionSkinningConsistency) {
         << ", reduced=" << reducedVertex.transpose();
   }
 }
+
+// Helper to add polygon data to a mesh for testing.
+// Creates polygons from existing triangle faces by merging adjacent triangle pairs into quads.
+void addTestPolyData(Mesh& mesh) {
+  mesh.polyFaces.clear();
+  mesh.polyFaceSizes.clear();
+  mesh.polyTexcoordFaces.clear();
+
+  // Merge pairs of triangles into quads; leave any leftover triangle as-is.
+  // For each pair of adjacent triangles (2*i, 2*i+1), create a quad from the unique vertices.
+  for (size_t i = 0; i + 1 < mesh.faces.size(); i += 2) {
+    const auto& f0 = mesh.faces[i];
+    const auto& f1 = mesh.faces[i + 1];
+
+    // Quad from the 4 unique vertices of two triangles sharing an edge
+    mesh.polyFaceSizes.push_back(4);
+    mesh.polyFaces.push_back(static_cast<uint32_t>(f0[0]));
+    mesh.polyFaces.push_back(static_cast<uint32_t>(f0[2]));
+    mesh.polyFaces.push_back(static_cast<uint32_t>(f1[2]));
+    mesh.polyFaces.push_back(static_cast<uint32_t>(f1[0]));
+  }
+
+  // If odd number of triangles, add the last one as a 3-gon
+  if (mesh.faces.size() % 2 == 1) {
+    const auto& f = mesh.faces.back();
+    mesh.polyFaceSizes.push_back(3);
+    mesh.polyFaces.push_back(static_cast<uint32_t>(f[0]));
+    mesh.polyFaces.push_back(static_cast<uint32_t>(f[1]));
+    mesh.polyFaces.push_back(static_cast<uint32_t>(f[2]));
+  }
+}
+
+// Test that reduceMeshByVertices preserves and correctly remaps polygon data
+TYPED_TEST(CharacterTest, ReduceMeshByVerticesPreservesPolygonData) {
+  // Add polygon data to the mesh
+  auto meshWithPolys = std::make_unique<Mesh>(*this->character.mesh);
+  addTestPolyData(*meshWithPolys);
+  const size_t originalPolyCount = meshWithPolys->polyFaceSizes.size();
+  ASSERT_GT(originalPolyCount, 0);
+
+  // Use the standalone mesh version for simplicity
+  const size_t vertexCount = meshWithPolys->vertices.size();
+
+  // Keep the first half of vertices
+  std::vector<bool> activeVertices(vertexCount, false);
+  for (size_t i = 0; i < vertexCount / 2; ++i) {
+    activeVertices[i] = true;
+  }
+
+  MeshT<float> reduced = reduceMeshByVertices(*meshWithPolys, activeVertices);
+
+  // Verify polygon data was filtered
+  EXPECT_LE(reduced.polyFaceSizes.size(), originalPolyCount);
+
+  // Verify all polygon vertex indices are valid in the reduced mesh
+  uint32_t offset = 0;
+  for (size_t polyIdx = 0; polyIdx < reduced.polyFaceSizes.size(); ++polyIdx) {
+    const auto polySize = reduced.polyFaceSizes[polyIdx];
+    EXPECT_GE(polySize, 3u);
+    for (uint32_t i = 0; i < polySize; ++i) {
+      EXPECT_LT(reduced.polyFaces[offset + i], static_cast<uint32_t>(reduced.vertices.size()));
+    }
+    offset += polySize;
+  }
+  EXPECT_EQ(offset, reduced.polyFaces.size());
+}
+
+// Test that reduceMeshByFaces preserves and correctly remaps polygon data
+TYPED_TEST(CharacterTest, ReduceMeshByFacesPreservesPolygonData) {
+  auto meshWithPolys = std::make_unique<Mesh>(*this->character.mesh);
+  addTestPolyData(*meshWithPolys);
+  const size_t originalPolyCount = meshWithPolys->polyFaceSizes.size();
+  ASSERT_GT(originalPolyCount, 0);
+
+  // Keep the first half of faces
+  std::vector<bool> activeFaces(meshWithPolys->faces.size(), false);
+  for (size_t i = 0; i < meshWithPolys->faces.size() / 2; ++i) {
+    activeFaces[i] = true;
+  }
+
+  MeshT<float> reduced = reduceMeshByFaces(*meshWithPolys, activeFaces);
+
+  // Verify all polygon vertex indices are valid
+  uint32_t offset = 0;
+  for (size_t polyIdx = 0; polyIdx < reduced.polyFaceSizes.size(); ++polyIdx) {
+    const auto polySize = reduced.polyFaceSizes[polyIdx];
+    EXPECT_GE(polySize, 3u);
+    for (uint32_t i = 0; i < polySize; ++i) {
+      EXPECT_LT(reduced.polyFaces[offset + i], static_cast<uint32_t>(reduced.vertices.size()));
+    }
+    offset += polySize;
+  }
+  EXPECT_EQ(offset, reduced.polyFaces.size());
+}
+
+// Test reduceMeshByPolys: select some polygons and verify triangle faces and polygon data
+TYPED_TEST(CharacterTest, ReduceMeshByPolys) {
+  auto meshWithPolys = std::make_unique<Mesh>(*this->character.mesh);
+  addTestPolyData(*meshWithPolys);
+  const size_t polyCount = meshWithPolys->polyFaceSizes.size();
+  ASSERT_GT(polyCount, 2u);
+
+  // Select the first two polygons only
+  std::vector<bool> activePolys(polyCount, false);
+  activePolys[0] = true;
+  activePolys[1] = true;
+
+  MeshT<float> reduced = reduceMeshByPolys(*meshWithPolys, activePolys);
+
+  // Should have exactly 2 polygons
+  EXPECT_EQ(reduced.polyFaceSizes.size(), 2u);
+
+  // Verify polygon vertex indices are valid
+  uint32_t offset = 0;
+  for (size_t polyIdx = 0; polyIdx < reduced.polyFaceSizes.size(); ++polyIdx) {
+    const auto polySize = reduced.polyFaceSizes[polyIdx];
+    for (uint32_t i = 0; i < polySize; ++i) {
+      EXPECT_LT(reduced.polyFaces[offset + i], static_cast<uint32_t>(reduced.vertices.size()));
+    }
+    offset += polySize;
+  }
+  EXPECT_EQ(offset, reduced.polyFaces.size());
+
+  // Should also have some triangle faces (the triangulated versions of those polys)
+  EXPECT_GT(reduced.faces.size(), 0u);
+
+  // All triangle face indices should be valid
+  for (const auto& face : reduced.faces) {
+    EXPECT_LT(face[0], static_cast<int>(reduced.vertices.size()));
+    EXPECT_LT(face[1], static_cast<int>(reduced.vertices.size()));
+    EXPECT_LT(face[2], static_cast<int>(reduced.vertices.size()));
+    EXPECT_GE(face[0], 0);
+    EXPECT_GE(face[1], 0);
+    EXPECT_GE(face[2], 0);
+  }
+
+  // Verify that vertex count is compact (no unused vertices)
+  std::vector<bool> usedVertices(reduced.vertices.size(), false);
+  for (const auto& face : reduced.faces) {
+    usedVertices[face[0]] = true;
+    usedVertices[face[1]] = true;
+    usedVertices[face[2]] = true;
+  }
+  offset = 0;
+  for (size_t polyIdx = 0; polyIdx < reduced.polyFaceSizes.size(); ++polyIdx) {
+    for (uint32_t i = 0; i < reduced.polyFaceSizes[polyIdx]; ++i) {
+      usedVertices[reduced.polyFaces[offset + i]] = true;
+    }
+    offset += reduced.polyFaceSizes[polyIdx];
+  }
+  for (size_t i = 0; i < usedVertices.size(); ++i) {
+    EXPECT_TRUE(usedVertices[i]) << "Vertex " << i << " is unused in reduced mesh";
+  }
+}
diff --git a/momentum/test/math/mesh_test.cpp b/momentum/test/math/mesh_test.cpp
index 34c3824a63..d0e64fe540 100644
--- a/momentum/test/math/mesh_test.cpp
+++ b/momentum/test/math/mesh_test.cpp
@@ -38,6 +38,9 @@ TYPED_TEST(MeshTest, DefaultConstructor) {
   EXPECT_TRUE(mesh.texcoords.empty());
   EXPECT_TRUE(mesh.texcoord_faces.empty());
   EXPECT_TRUE(mesh.texcoord_lines.empty());
+  EXPECT_TRUE(mesh.polyFaces.empty());
+  EXPECT_TRUE(mesh.polyFaceSizes.empty());
+  EXPECT_TRUE(mesh.polyTexcoordFaces.empty());
 }
 
 // Test updateNormals with a tetrahedron
@@ -166,6 +169,11 @@ TYPED_TEST(MeshTest, Cast) {
 
   mesh.confidence = {0.5, 0.6, 0.7};
 
+  // Add polygon data: a quad (v0,v1,v2,v3) + a triangle (v0,v1,v2)
+  mesh.polyFaces = {0, 1, 2, 0, 0, 1, 2};
+  mesh.polyFaceSizes = {4, 3};
+  mesh.polyTexcoordFaces = {10, 11, 12, 13, 14, 15, 16};
+
   // Cast to the same type (should be a copy)
   auto sameMesh = mesh.template cast<T>();
 
@@ -179,6 +187,9 @@ TYPED_TEST(MeshTest, Cast) {
     EXPECT_TRUE(sameMesh.normals[i].isApprox(mesh.normals[i]));
     EXPECT_NEAR(sameMesh.confidence[i], mesh.confidence[i], Eps<T>(1e-5f, 1e-13));
   }
+  EXPECT_EQ(sameMesh.polyFaces, mesh.polyFaces);
+  EXPECT_EQ(sameMesh.polyFaceSizes, mesh.polyFaceSizes);
+  EXPECT_EQ(sameMesh.polyTexcoordFaces, mesh.polyTexcoordFaces);
 
   // Cast to the other type
   using OtherT = typename std::conditional<std::is_same<T, float>::value, double, float>::type;
@@ -195,6 +206,9 @@ TYPED_TEST(MeshTest, Cast) {
     // Use a more relaxed tolerance for float-to-double or double-to-float conversions
     EXPECT_NEAR(static_cast<T>(otherMesh.confidence[i]), mesh.confidence[i], Eps<T>(1e-4f, 1e-6));
   }
+  EXPECT_EQ(otherMesh.polyFaces, mesh.polyFaces);
+  EXPECT_EQ(otherMesh.polyFaceSizes, mesh.polyFaceSizes);
+  EXPECT_EQ(otherMesh.polyTexcoordFaces, mesh.polyTexcoordFaces);
 }
 
 // Test reset method
@@ -225,6 +239,10 @@ TYPED_TEST(MeshTest, Reset) {
 
   mesh.texcoord_lines = {{0, 1}, {1, 2}, {2, 0}};
 
+  mesh.polyFaces = {0, 1, 2, 0, 0, 1, 2};
+  mesh.polyFaceSizes = {4, 3};
+  mesh.polyTexcoordFaces = {10, 11, 12, 13, 14, 15, 16};
+
   // Reset the mesh
   mesh.reset();
 
@@ -238,6 +256,9 @@ TYPED_TEST(MeshTest, Reset) {
   EXPECT_TRUE(mesh.texcoords.empty());
   EXPECT_TRUE(mesh.texcoord_faces.empty());
   EXPECT_TRUE(mesh.texcoord_lines.empty());
+  EXPECT_TRUE(mesh.polyFaces.empty());
+  EXPECT_TRUE(mesh.polyFaceSizes.empty());
+  EXPECT_TRUE(mesh.polyTexcoordFaces.empty());
 }
 
 // Test with a complex mesh
@@ -351,6 +372,9 @@ TYPED_TEST(MeshTest, EmptyMesh) {
   EXPECT_TRUE(castedMesh.texcoords.empty());
   EXPECT_TRUE(castedMesh.texcoord_faces.empty());
   EXPECT_TRUE(castedMesh.texcoord_lines.empty());
+  EXPECT_TRUE(castedMesh.polyFaces.empty());
+  EXPECT_TRUE(castedMesh.polyFaceSizes.empty());
+  EXPECT_TRUE(castedMesh.polyTexcoordFaces.empty());
 
   // Reset an empty mesh
   EXPECT_NO_THROW(mesh.reset());
@@ -462,3 +486,66 @@ TYPED_TEST(MeshTest, Lines) {
   mesh.reset();
   EXPECT_TRUE(mesh.lines.empty());
 }
+
+// Test polygon face representation
+TYPED_TEST(MeshTest, PolygonFaces) {
+  using T = TypeParam;
+  using MeshType = typename TestFixture::MeshType;
+  using Vector3 = Eigen::Vector3<T>;
+
+  MeshType mesh;
+
+  // Create a mesh with 5 vertices
+  mesh.vertices = {
+      Vector3(0, 0, 0), Vector3(1, 0, 0), Vector3(1, 1, 0), Vector3(0, 1, 0), Vector3(0.5, 0.5, 1)};
+
+  // Polygon representation: a quad (4 verts) + a triangle (3 verts)
+  mesh.polyFaces = {0, 1, 2, 3, 0, 1, 4};
+  mesh.polyFaceSizes = {4, 3};
+  mesh.polyTexcoordFaces = {0, 1, 2, 3, 4, 5, 6};
+
+  // Verify sizes sum matches polyFaces length
+  uint32_t totalSize = 0;
+  for (auto s : mesh.polyFaceSizes) {
+    totalSize += s;
+  }
+  EXPECT_EQ(totalSize, mesh.polyFaces.size());
+  EXPECT_EQ(totalSize, mesh.polyTexcoordFaces.size());
+
+  // Cast and verify preservation
+  using OtherT = typename std::conditional<std::is_same<T, float>::value, double, float>::type;
+  auto otherMesh = mesh.template cast<OtherT>();
+  EXPECT_EQ(otherMesh.polyFaces, mesh.polyFaces);
+  EXPECT_EQ(otherMesh.polyFaceSizes, mesh.polyFaceSizes);
+  EXPECT_EQ(otherMesh.polyTexcoordFaces, mesh.polyTexcoordFaces);
+
+  // Reset and verify cleared
+  mesh.reset();
+  EXPECT_TRUE(mesh.polyFaces.empty());
+  EXPECT_TRUE(mesh.polyFaceSizes.empty());
+  EXPECT_TRUE(mesh.polyTexcoordFaces.empty());
+}
+
+// Test polygon faces with empty texcoord faces
+TYPED_TEST(MeshTest, PolygonFacesNoTexcoords) {
+  using T = TypeParam;
+  using MeshType = typename TestFixture::MeshType;
+  using Vector3 = Eigen::Vector3<T>;
+
+  MeshType mesh;
+
+  mesh.vertices = {Vector3(0, 0, 0), Vector3(1, 0, 0), Vector3(0, 1, 0)};
+  mesh.polyFaces = {0, 1, 2};
+  mesh.polyFaceSizes = {3};
+  // polyTexcoordFaces intentionally left empty
+
+  EXPECT_EQ(mesh.polyFaces.size(), 3u);
+  EXPECT_EQ(mesh.polyFaceSizes.size(), 1u);
+  EXPECT_TRUE(mesh.polyTexcoordFaces.empty());
+
+  // Cast preserves the empty texcoord faces
+  auto castedMesh = mesh.template cast<T>();
+  EXPECT_EQ(castedMesh.polyFaces, mesh.polyFaces);
+  EXPECT_EQ(castedMesh.polyFaceSizes, mesh.polyFaceSizes);
+  EXPECT_TRUE(castedMesh.polyTexcoordFaces.empty());
+}
diff --git a/pymomentum/geometry/mesh_pybind.cpp b/pymomentum/geometry/mesh_pybind.cpp
index 6e57114331..559b01ebda 100644
--- a/pymomentum/geometry/mesh_pybind.cpp
+++ b/pymomentum/geometry/mesh_pybind.cpp
@@ -29,6 +29,52 @@ namespace pymomentum {
 
 namespace {
 
+void validateAndSetPolyData(
+    mm::Mesh& mesh,
+    const std::vector<uint32_t>& poly_faces,
+    const std::vector<uint32_t>& poly_face_sizes,
+    const std::vector<uint32_t>& poly_texcoord_faces,
+    int nVerts,
+    int nTextureCoords) {
+  MT_THROW_IF(
+      poly_faces.empty() != poly_face_sizes.empty(),
+      "poly_faces and poly_face_sizes must both be empty or both be non-empty");
+  if (!poly_face_sizes.empty()) {
+    uint32_t totalSize = 0;
+    for (const auto s : poly_face_sizes) {
+      MT_THROW_IF(s < 3, "poly_face_sizes entries must be >= 3, got {}", s);
+      totalSize += s;
+    }
+    MT_THROW_IF(
+        totalSize != poly_faces.size(),
+        "poly_face_sizes sum ({}) must equal poly_faces length ({})",
+        totalSize,
+        poly_faces.size());
+    for (const auto idx : poly_faces) {
+      MT_THROW_IF(
+          idx >= static_cast<uint32_t>(nVerts),
+          "poly_faces index ({}) exceeded vertex count ({})",
+          idx,
+          nVerts);
+    }
+    MT_THROW_IF(
+        !poly_texcoord_faces.empty() && poly_texcoord_faces.size() != poly_faces.size(),
+        "poly_texcoord_faces must be empty or same length as poly_faces ({} vs {})",
+        poly_texcoord_faces.size(),
+        poly_faces.size());
+    for (const auto idx : poly_texcoord_faces) {
+      MT_THROW_IF(
+          idx >= static_cast<uint32_t>(nTextureCoords),
+          "poly_texcoord_faces index ({}) exceeded texcoord count ({})",
+          idx,
+          nTextureCoords);
+    }
+  }
+  mesh.polyFaces = poly_faces;
+  mesh.polyFaceSizes = poly_face_sizes;
+  mesh.polyTexcoordFaces = poly_texcoord_faces;
+}
+
 mm::Mesh createMesh(
     const py::array_t<float>& vertices,
     const py::array_t<int>& faces,
@@ -38,7 +84,10 @@ mm::Mesh createMesh(
     const std::vector<float>& confidence,
     std::optional<py::array_t<float>> texcoords,
     std::optional<py::array_t<int>> texcoord_faces,
-    const std::vector<std::vector<int32_t>>& texcoord_lines) {
+    const std::vector<std::vector<int32_t>>& texcoord_lines,
+    const std::vector<uint32_t>& poly_faces,
+    const std::vector<uint32_t>& poly_face_sizes,
+    const std::vector<uint32_t>& poly_texcoord_faces) {
   mm::Mesh mesh;
   MT_THROW_IF(vertices.ndim() != 2, "vertices must be a 2D array");
   MT_THROW_IF(vertices.shape(1) != 3, "vertices must have size n x 3");
@@ -110,6 +159,9 @@ mm::Mesh createMesh(
 
   mesh.texcoord_lines = texcoord_lines;
 
+  validateAndSetPolyData(
+      mesh, poly_faces, poly_face_sizes, poly_texcoord_faces, static_cast<int>(nVerts), nTextureCoords);
+
   return mesh;
 }
 
@@ -136,6 +188,9 @@ void registerMeshBindings(py::class_<mm::Mesh>& meshClass) {
 :param texcoords: Optional n x 2 array of texture coordinates.
 :param texcoord_faces: Optional n x 3 array of triangles in the texture map.  Each triangle corresponds to a triangle on the mesh, but indices should refer to the texcoord array.
 :param texcoord_lines: Optional list of lines, where each line is a list of texture coordinate indices.
+:param poly_faces: Optional list of packed polygon face vertex indices (all polygons concatenated). The triangulated ``faces`` parameter is always required; polygon data is optional and preserves original topology.
+:param poly_face_sizes: Optional list of vertex counts per polygon face.
+:param poly_texcoord_faces: Optional list of packed polygon face texture coordinate indices (same layout as poly_faces).
           )",
           py::arg("vertices"),
           py::arg("faces"),
@@ -146,7 +201,10 @@ void registerMeshBindings(py::class_<mm::Mesh>& meshClass) {
           py::arg("confidence") = std::vector<float>{},
           py::arg("texcoords") = std::optional<py::array_t<float>>{},
           py::arg("texcoord_faces") = std::optional<py::array_t<int>>{},
-          py::arg("texcoord_lines") = std::vector<std::vector<int32_t>>{})
+          py::arg("texcoord_lines") = std::vector<std::vector<int32_t>>{},
+          py::arg("poly_faces") = std::vector<uint32_t>{},
+          py::arg("poly_face_sizes") = std::vector<uint32_t>{},
+          py::arg("poly_texcoord_faces") = std::vector<uint32_t>{})
       .def_property_readonly(
           "n_vertices",
           [](const mm::Mesh& mesh) { return mesh.vertices.size(); },
@@ -188,6 +246,24 @@ void registerMeshBindings(py::class_<mm::Mesh>& meshClass) {
           "texcoord_lines",
           &mm::Mesh::texcoord_lines,
           "Texture coordinate indices for each line.  ")
+      .def_readonly(
+          "poly_faces",
+          &mm::Mesh::polyFaces,
+          "Packed polygon face vertex indices (all polygons concatenated back-to-back). "
+          "This is optional — the triangulated :attr:`faces` field is always required.")
+      .def_readonly(
+          "poly_face_sizes",
+          &mm::Mesh::polyFaceSizes,
+          "Number of vertices in each polygon face. May be empty if polygon data is not available.")
+      .def_readonly(
+          "poly_texcoord_faces",
+          &mm::Mesh::polyTexcoordFaces,
+          "Packed polygon face texture coordinate indices (same layout as :attr:`poly_faces`). "
+          "May be empty if texture coordinates are not available.")
+      .def_property_readonly(
+          "n_poly_faces",
+          [](const mm::Mesh& mesh) { return mesh.polyFaceSizes.size(); },
+          ":return: The number of polygon faces in the mesh.")
       .def(
           "self_intersections",
           [](const mm::Mesh& mesh) {
@@ -204,12 +280,13 @@ void registerMeshBindings(py::class_<mm::Mesh>& meshClass) {
           })
       .def("__repr__", [](const mm::Mesh& m) {
         return fmt::format(
-            "Mesh(vertices={}, faces={}, has_normals={}, has_colors={}, has_texcoords={})",
+            "Mesh(vertices={}, faces={}, has_normals={}, has_colors={}, has_texcoords={}, poly_faces={})",
             m.vertices.size(),
             m.faces.size(),
             !m.normals.empty() ? "True" : "False",
             !m.colors.empty() ? "True" : "False",
-            !m.texcoords.empty() ? "True" : "False");
+            !m.texcoords.empty() ? "True" : "False",
+            m.polyFaceSizes.size());
       });
 }