diff --git a/examples/triangulation_method_example.cc b/examples/triangulation_method_example.cc
new file mode 100644
index 00000000..0ec6c9cf
--- /dev/null
+++ b/examples/triangulation_method_example.cc
@@ -0,0 +1,124 @@
+// Example demonstrating the new triangulation method option in TinyUSDZ Tydra
+//
+// This example shows how to use either:
+// - Earcut algorithm (default): Robust for complex/concave polygons
+// - Triangle Fan: Fast for simple convex polygons
+
+#include <iostream>
+#include <vector>
+#include <string>
+#include "tinyusdz.hh"
+#include "tydra/render-data.hh"
+#include "tydra/scene-access.hh"
+
+int main(int argc, char** argv) {
+    if (argc < 2) {
+        std::cout << "Usage: " << argv[0] << " <input.usd> [fan|earcut]\n";
+        std::cout << "\nTriangulation methods:\n";
+        std::cout << "  earcut (default): Robust algorithm for complex polygons\n";
+        std::cout << "  fan: Fast triangle fan for convex polygons\n";
+        return 1;
+    }
+
+    std::string filename = argv[1];
+    bool use_fan = false;
+
+    if (argc >= 3) {
+        std::string method = argv[2];
+        if (method == "fan") {
+            use_fan = true;
+            std::cout << "Using triangle fan triangulation\n";
+        } else if (method == "earcut") {
+            std::cout << "Using earcut triangulation\n";
+        } else {
+            std::cerr << "Unknown triangulation method: " << method << "\n";
+            return 1;
+        }
+    } else {
+        std::cout << "Using default earcut triangulation\n";
+    }
+
+    // Load USD file
+    tinyusdz::Stage stage;
+    std::string warn, err;
+    bool ret = tinyusdz::LoadUSDFromFile(filename, &stage, &warn, &err);
+
+    if (!warn.empty()) {
+        std::cerr << "Warning: " << warn << "\n";
+    }
+
+    if (!ret) {
+        std::cerr << "Failed to load USD file: " << err << "\n";
+        return 1;
+    }
+
+    // Set up render scene converter with triangulation method
+    tinyusdz::tydra::RenderSceneConverterEnv env(stage);
+
+    // Configure mesh conversion
+    env.mesh_config.triangulate = true;  // Enable triangulation
+
+    // Set the triangulation method
+    if (use_fan) {
+        env.mesh_config.triangulation_method =
+            tinyusdz::tydra::MeshConverterConfig::TriangulationMethod::TriangleFan;
+    } else {
+        env.mesh_config.triangulation_method =
+            tinyusdz::tydra::MeshConverterConfig::TriangulationMethod::Earcut;
+    }
+
+    // Convert to render scene
+    tinyusdz::tydra::RenderSceneConverter converter;
+    tinyusdz::tydra::RenderScene render_scene;
+
+    if (!converter.ConvertToRenderScene(env, &render_scene)) {
+        std::cerr << "Failed to convert to render scene: "
+                  << converter.GetError() << "\n";
+        return 1;
+    }
+
+    // Print statistics
+    std::cout << "\nConversion complete!\n";
+    std::cout << "Number of meshes: " << render_scene.meshes.size() << "\n";
+
+    size_t total_triangles = 0;
+    size_t total_original_faces = 0;
+
+    for (const auto& mesh : render_scene.meshes) {
+        if (!mesh.triangulatedFaceVertexIndices.empty()) {
+            size_t num_triangles = mesh.triangulatedFaceVertexIndices.size() / 3;
+            size_t num_original_faces = mesh.usdFaceVertexCounts.size();
+
+            total_triangles += num_triangles;
+            total_original_faces += num_original_faces;
+
+            std::cout << "\nMesh: " << mesh.prim_name << "\n";
+            std::cout << "  Original faces: " << num_original_faces << "\n";
+            std::cout << "  Triangulated faces: " << num_triangles << "\n";
+
+            // Count polygon types
+            std::map<uint32_t, size_t> poly_counts;
+            for (auto count : mesh.usdFaceVertexCounts) {
+                poly_counts[count]++;
+            }
+
+            std::cout << "  Original polygon distribution:\n";
+            for (const auto& kv : poly_counts) {
+                std::cout << "    " << kv.first << "-gons: " << kv.second << "\n";
+            }
+        }
+    }
+
+    std::cout << "\nTotal original faces: " << total_original_faces << "\n";
+    std::cout << "Total triangles: " << total_triangles << "\n";
+
+    if (use_fan) {
+        std::cout << "\nNote: Triangle fan was used for 5+ vertex polygons.\n";
+        std::cout << "This is faster but assumes polygons are convex.\n";
+    } else {
+        std::cout << "\nNote: Earcut algorithm was used for 5+ vertex polygons.\n";
+        std::cout << "This handles complex polygons correctly.\n";
+    }
+
+    return 0;
+}
\ No newline at end of file
diff --git a/src/tydra/render-data.cc b/src/tydra/render-data.cc
index e3cd5c21..b0859c7a 100644
--- a/src/tydra/render-data.cc
+++ b/src/tydra/render-data.cc
@@ -1649,7 +1649,9 @@ bool TriangulatePolygon(
     std::vector<uint32_t> &triangulatedFaceVertexCounts,
     std::vector<uint32_t> &triangulatedFaceVertexIndices,
     std::vector<size_t> &triangulatedToOrigFaceVertexIndexMap,
-    std::vector<uint32_t> &triangulatedFaceCounts, std::string &warn, std::string &err) {
+    std::vector<uint32_t> &triangulatedFaceCounts,
+    MeshConverterConfig::TriangulationMethod triangulation_method,
+    std::string &warn, std::string &err) {
   triangulatedFaceVertexCounts.clear();
   triangulatedFaceVertexIndices.clear();
 
@@ -1720,143 +1722,172 @@ bool TriangulatePolygon(
       triangulatedFaceCounts.push_back(2);
 #endif
     } else {
-      // Use double for accuracy. `float` precision may classify small-are polygon as degenerated.
-      // Find the normal axis of the polygon using Newell's method
-      value::double3 n = {0, 0, 0};
+      // Polygon with 5+ vertices
+      if (triangulation_method == MeshConverterConfig::TriangulationMethod::TriangleFan) {
+        // Simple triangle fan triangulation
+        // This assumes the polygon is convex
+        // Creates triangles: (0,1,2), (0,2,3), (0,3,4), ...
 
-      size_t vi0;
-      size_t vi0_2;
+        size_t ntris = npolys - 2;
 
-      //std::cout << "npoly " << npolys << "\n";
-
-      for (size_t k = 0; k < npolys; ++k) {
-        vi0 = faceVertexIndices[faceIndexOffset + k];
-
-        size_t j = (k + 1) % npolys;
-        vi0_2 = faceVertexIndices[faceIndexOffset + j];
-
-        if (vi0 >= points.size()) {
-          err = fmt::format("Invalid vertex index.\n");
-          return false;
-        }
-
-        if (vi0_2 >= points.size()) {
-          err = fmt::format("Invalid vertex index.\n");
-          return false;
-        }
-
-        T v0 = points[vi0];
-        T v1 = points[vi0_2];
-
-        const T point1 = {v0[0], v0[1], v0[2]};
-        const T point2 = {v1[0], v1[1], v1[2]};
-
-        T a = {point1[0] - point2[0], point1[1] - point2[1],
-               point1[2] - point2[2]};
-        T b = {point1[0] + point2[0], point1[1] + point2[1],
-               point1[2] + point2[2]};
-
-        n[0] += double(a[1] * b[2]);
-        n[1] += double(a[2] * b[0]);
-        n[2] += double(a[0] * b[1]);
-        DCOUT("v0 " << v0);
-        DCOUT("v1 " << v1);
-        DCOUT("n " << n);
-      }
-      //BaseTy length_n = vlength(n);
-      double length_n = vlength(n);
-
-      // Check if zero length normal
-      if (std::fabs(length_n) < std::numeric_limits<double>::epsilon()) {
-        DCOUT("length_n " << length_n);
-        err = "Degenerated polygon found.\n";
-        return false;
-      }
-
-      // Negative is to flip the normal to the correct direction
-      n = vnormalize(n);
-
-      T axis_w, axis_v, axis_u;
-      axis_w[0] = BaseTy(n[0]);
-      axis_w[1] = BaseTy(n[1]);
-      axis_w[2] = BaseTy(n[2]);
-      T a;
-      if (std::fabs(axis_w[0]) > BaseTy(0.9999999)) {  // TODO: use 1.0 - eps?
-        a = {BaseTy(0), BaseTy(1), BaseTy(0)};
-      } else {
-        a = {BaseTy(1), BaseTy(0), BaseTy(0)};
-      }
-      axis_v = vnormalize(vcross(axis_w, a));
-      axis_u = vcross(axis_w, axis_v);
-
-      using Point3D = std::array<BaseTy, 3>;
-      using Point2D = std::array<BaseTy, 2>;
-      std::vector<Point2D> polyline;
-
-      // TMW change: Find best normal and project v0x and v0y to those
-      // coordinates, instead of picking a plane aligned with an axis (which
-      // can flip polygons).
-
-      // Fill polygon data.
-      for (size_t k = 0; k < npolys; k++) {
-        size_t vidx = faceVertexIndices[faceIndexOffset + k];
-
-        value::float3 v = points[vidx];
-        // Point3 polypoint = {v0[0],v0[1],v0[2]};
-
-        // world to local
-        Point3D loc = {vdot(v, axis_u), vdot(v, axis_v), vdot(v, axis_w)};
-
-        polyline.push_back({loc[0], loc[1]});
-      }
-
-      std::vector<std::vector<Point2D>> polygon_2d;
-      polygon_2d.push_back(polyline);
-      // Single polygon only(no holes)
-
-      std::vector<uint32_t> indices = mapbox::earcut<uint32_t>(polygon_2d);
-      //  => result = 3 * faces, clockwise
-
-      if (indices.empty()) {
-        warn += "Failed to triangualte a polygon. input is not CCW, have holes or invalid topology.\n";
-
-        //DumpTriangle(points, indices);
-      }
-
-      if ((indices.size() % 3) != 0) {
-        // This should not be happen, though.
-        err = "Failed to triangulate.\n";
-        return false;
-      }
-
-      size_t ntris = indices.size() / 3;
-      //std::cout << "ntris " << ntris << "\n";
-
-
-      // Up to 2GB tris.
-      if (ntris > size_t((std::numeric_limits<int32_t>::max)())) {
-        err = "Too many triangles are generated.\n";
-        return false;
-      }
-
-      if (ntris > 0) {
         for (size_t k = 0; k < ntris; k++) {
           triangulatedFaceVertexCounts.push_back(3);
-          triangulatedFaceVertexIndices.push_back(
-              faceVertexIndices[faceIndexOffset + indices[3 * k + 0]]);
-          triangulatedFaceVertexIndices.push_back(
-              faceVertexIndices[faceIndexOffset + indices[3 * k + 1]]);
-          triangulatedFaceVertexIndices.push_back(
-              faceVertexIndices[faceIndexOffset + indices[3 * k + 2]]);
 
-          triangulatedToOrigFaceVertexIndexMap.push_back(faceIndexOffset +
-                                                         indices[3 * k + 0]);
-          triangulatedToOrigFaceVertexIndexMap.push_back(faceIndexOffset +
-                                                         indices[3 * k + 1]);
-          triangulatedToOrigFaceVertexIndexMap.push_back(faceIndexOffset +
-                                                         indices[3 * k + 2]);
+          // First vertex is always the pivot (index 0)
+          triangulatedFaceVertexIndices.push_back(
+              faceVertexIndices[faceIndexOffset + 0]);
+          triangulatedFaceVertexIndices.push_back(
+              faceVertexIndices[faceIndexOffset + k + 1]);
+          triangulatedFaceVertexIndices.push_back(
+              faceVertexIndices[faceIndexOffset + k + 2]);
+
+          triangulatedToOrigFaceVertexIndexMap.push_back(faceIndexOffset + 0);
+          triangulatedToOrigFaceVertexIndexMap.push_back(faceIndexOffset + k + 1);
+          triangulatedToOrigFaceVertexIndexMap.push_back(faceIndexOffset + k + 2);
         }
+
         triangulatedFaceCounts.push_back(uint32_t(ntris));
+
+      } else {
+        // Use earcut algorithm (default, handles complex polygons)
+        // Use double for accuracy. `float` precision may classify small-are polygon as degenerated.
+        // Find the normal axis of the polygon using Newell's method
+        value::double3 n = {0, 0, 0};
+
+        size_t vi0;
+        size_t vi0_2;
+
+        //std::cout << "npoly " << npolys << "\n";
+
+        for (size_t k = 0; k < npolys; ++k) {
+          vi0 = faceVertexIndices[faceIndexOffset + k];
+
+          size_t j = (k + 1) % npolys;
+          vi0_2 = faceVertexIndices[faceIndexOffset + j];
+
+          if (vi0 >= points.size()) {
+            err = fmt::format("Invalid vertex index.\n");
+            return false;
+          }
+
+          if (vi0_2 >= points.size()) {
+            err = fmt::format("Invalid vertex index.\n");
+            return false;
+          }
+
+          T v0 = points[vi0];
+          T v1 = points[vi0_2];
+
+          const T point1 = {v0[0], v0[1], v0[2]};
+          const T point2 = {v1[0], v1[1], v1[2]};
+
+          T a = {point1[0] - point2[0], point1[1] - point2[1],
+                 point1[2] - point2[2]};
+          T b = {point1[0] + point2[0], point1[1] + point2[1],
+                 point1[2] + point2[2]};
+
+          n[0] += double(a[1] * b[2]);
+          n[1] += double(a[2] * b[0]);
+          n[2] += double(a[0] * b[1]);
+          DCOUT("v0 " << v0);
+          DCOUT("v1 " << v1);
+          DCOUT("n " << n);
+        }
+        //BaseTy length_n = vlength(n);
+        double length_n = vlength(n);
+
+        // Check if zero length normal
+        if (std::fabs(length_n) < std::numeric_limits<double>::epsilon()) {
+          DCOUT("length_n " << length_n);
+          err = "Degenerated polygon found.\n";
+          return false;
+        }
+
+        // Negative is to flip the normal to the correct direction
+        n = vnormalize(n);
+
+        T axis_w, axis_v, axis_u;
+        axis_w[0] = BaseTy(n[0]);
+        axis_w[1] = BaseTy(n[1]);
+        axis_w[2] = BaseTy(n[2]);
+        T a;
+        if (std::fabs(axis_w[0]) > BaseTy(0.9999999)) {  // TODO: use 1.0 - eps?
+          a = {BaseTy(0), BaseTy(1), BaseTy(0)};
+        } else {
+          a = {BaseTy(1), BaseTy(0), BaseTy(0)};
+        }
+        axis_v = vnormalize(vcross(axis_w, a));
+        axis_u = vcross(axis_w, axis_v);
+
+        using Point3D = std::array<BaseTy, 3>;
+        using Point2D = std::array<BaseTy, 2>;
+        std::vector<Point2D> polyline;
+
+        // TMW change: Find best normal and project v0x and v0y to those
+        // coordinates, instead of picking a plane aligned with an axis (which
+        // can flip polygons).
+
+        // Fill polygon data.
+        for (size_t k = 0; k < npolys; k++) {
+          size_t vidx = faceVertexIndices[faceIndexOffset + k];
+
+          value::float3 v = points[vidx];
+          // Point3 polypoint = {v0[0],v0[1],v0[2]};
+
+          // world to local
+          Point3D loc = {vdot(v, axis_u), vdot(v, axis_v), vdot(v, axis_w)};
+
+          polyline.push_back({loc[0], loc[1]});
+        }
+
+        std::vector<std::vector<Point2D>> polygon_2d;
+        polygon_2d.push_back(polyline);
+        // Single polygon only(no holes)
+
+        std::vector<uint32_t> indices = mapbox::earcut<uint32_t>(polygon_2d);
+        //  => result = 3 * faces, clockwise
+
+        if (indices.empty()) {
+          warn += "Failed to triangualte a polygon. input is not CCW, have holes or invalid topology.\n";
+
+          //DumpTriangle(points, indices);
+        }
+
+        if ((indices.size() % 3) != 0) {
+          // This should not be happen, though.
+          err = "Failed to triangulate.\n";
+          return false;
+        }
+
+        size_t ntris = indices.size() / 3;
+        //std::cout << "ntris " << ntris << "\n";
+
+
+        // Up to 2GB tris.
+        if (ntris > size_t((std::numeric_limits<int32_t>::max)())) {
+          err = "Too many triangles are generated.\n";
+          return false;
+        }
+
+        if (ntris > 0) {
+          for (size_t k = 0; k < ntris; k++) {
+            triangulatedFaceVertexCounts.push_back(3);
+            triangulatedFaceVertexIndices.push_back(
+                faceVertexIndices[faceIndexOffset + indices[3 * k + 0]]);
+            triangulatedFaceVertexIndices.push_back(
+                faceVertexIndices[faceIndexOffset + indices[3 * k + 1]]);
+            triangulatedFaceVertexIndices.push_back(
+                faceVertexIndices[faceIndexOffset + indices[3 * k + 2]]);
+
+            triangulatedToOrigFaceVertexIndexMap.push_back(faceIndexOffset +
+                                                           indices[3 * k + 0]);
+            triangulatedToOrigFaceVertexIndexMap.push_back(faceIndexOffset +
+                                                           indices[3 * k + 1]);
+            triangulatedToOrigFaceVertexIndexMap.push_back(faceIndexOffset +
+                                                           indices[3 * k + 2]);
+          }
+          triangulatedFaceCounts.push_back(uint32_t(ntris));
+        }
       }
     }
 
@@ -4217,6 +4248,7 @@ bool RenderSceneConverter::ConvertMesh(
             dst.points, dst.usdFaceVertexCounts, dst.usdFaceVertexIndices,
             triangulatedFaceVertexCounts, triangulatedFaceVertexIndices,
             triangulatedToOrigFaceVertexIndexMap, triangulatedFaceCounts,
+            env.mesh_config.triangulation_method,
             _warn, err)) {
       PUSH_ERROR_AND_RETURN("Triangulation failed: " + err);
     }
diff --git a/src/tydra/render-data.hh b/src/tydra/render-data.hh
index 6035e0ae..7194abed 100644
--- a/src/tydra/render-data.hh
+++ b/src/tydra/render-data.hh
@@ -1626,6 +1626,14 @@ bool DefaultTextureImageLoaderFunction(const value::AssetPath &assetPath,
 struct MeshConverterConfig {
   bool triangulate{true};
 
+  // Triangulation method for polygons with 5+ vertices
+  enum class TriangulationMethod {
+    Earcut,     // Use earcut algorithm (robust, handles complex polygons)
+    TriangleFan // Use simple triangle fan (faster, only for convex polygons)
+  };
+
+  TriangulationMethod triangulation_method{TriangulationMethod::Earcut};
+
   bool validate_geomsubset{true};  // Validate GeomSubset.
 
   // We may want texcoord data even if the Mesh does not have bound Material.