diff --git a/gdist.py b/gdist.py
index 8322bcd..b88cb3d 100644
--- a/gdist.py
+++ b/gdist.py
@@ -32,6 +32,7 @@
     np.ctypeslib.ndpointer(dtype=np.int32),
     np.ctypeslib.ndpointer(dtype=np.float64),
     ctypes.c_double,
+    ctypes.c_uint,
 ]
 lib.compute_gdist.restype = None
 
@@ -42,6 +43,7 @@
     np.ctypeslib.ndpointer(dtype=np.int32),
     ctypes.POINTER(ctypes.c_uint),
     ctypes.c_double,
+    ctypes.c_uint,
 ]
 lib.local_gdist_matrix.restype = ctypes.POINTER(ctypes.c_double)
 
@@ -63,6 +65,7 @@ def compute_gdist(
         source_indices_array,
         target_indices_array,
         distance_limit,
+        is_one_indexed,
     ):
         target_indices_size = target_indices_array.size
         distance = np.empty(target_indices_size, dtype=np.float64)
@@ -77,6 +80,7 @@ def compute_gdist(
             target_indices_array,
             distance,
             distance_limit,
+            is_one_indexed,
         )
         return distance
 
@@ -87,6 +91,7 @@ def local_gdist_matrix(
         vertices,
         triangles,
         max_distance,
+        is_one_indexed,
     ):
         sparse_matrix_size = ctypes.c_uint(0)
         data = lib.local_gdist_matrix(
@@ -96,6 +101,7 @@ def local_gdist_matrix(
             triangles,
             ctypes.byref(sparse_matrix_size),
             max_distance,
+            is_one_indexed,
         )
 
         np_data = np.fromiter(
@@ -113,6 +119,7 @@ def compute_gdist(
     source_indices=None,
     target_indices=None,
     max_distance=1e100,
+    is_one_indexed=False,
 ):
     vertices = vertices.ravel()
     triangles = triangles.ravel()
@@ -130,6 +137,7 @@ def compute_gdist(
         source_indices_array=source_indices,
         target_indices_array=target_indices,
         distance_limit=max_distance,
+        is_one_indexed=int(is_one_indexed),
     )
     return np.fromiter(distance, dtype=np.float64, count=target_indices.size)
 
@@ -138,17 +146,19 @@ def local_gdist_matrix(
     vertices,
     triangles,
     max_distance=1e100,
+    is_one_indexed=False,
 ):
     vertices = vertices.ravel()
     triangles = triangles.ravel()
 
     g = Gdist()
     data = g.local_gdist_matrix(
-        vertices.size,
-        triangles.size,
-        vertices,
-        triangles,
-        max_distance,
+        number_of_vertices=vertices.size,
+        number_of_triangles=triangles.size,
+        vertices=vertices,
+        triangles=triangles,
+        max_distance=max_distance,
+        is_one_indexed=is_one_indexed,
     )
     sizes = data.size // 3
     rows = data[:sizes]
diff --git a/gdist_c_api.cpp b/gdist_c_api.cpp
index 5635871..a0f0939 100644
--- a/gdist_c_api.cpp
+++ b/gdist_c_api.cpp
@@ -11,7 +11,8 @@ void compute_gdist_impl(
     unsigned *source_indices_array,
     unsigned *target_indices_array,
     double *distance,
-    double distance_limit
+    double distance_limit,
+    unsigned is_one_indexed
 ) {
 
     std::vector<double> points (vertices, vertices + number_of_vertices);
@@ -20,7 +21,7 @@ void compute_gdist_impl(
     std::vector<unsigned> target_indices (target_indices_array, target_indices_array + number_of_target_indices);
 
     geodesic::Mesh mesh;
-    mesh.initialize_mesh_data(points, faces); // create internal mesh data structure including edges
+    mesh.initialize_mesh_data(points, faces, is_one_indexed); // create internal mesh data structure including edges
 
     geodesic::GeodesicAlgorithmExact algorithm(&mesh); // create exact algorithm for the mesh
     
@@ -47,13 +48,14 @@ double* local_gdist_matrix_impl(
     double *vertices,
     unsigned *triangles,
     unsigned *sparse_matrix_size,
-    double max_distance
+    double max_distance,
+    unsigned is_one_indexed
 ) {
     std::vector<double> points (vertices, vertices + number_of_vertices);
     std::vector<unsigned> faces (triangles, triangles + number_of_triangles);
     
     geodesic::Mesh mesh;
-    mesh.initialize_mesh_data(points, faces); // create internal mesh data structure including edges
+    mesh.initialize_mesh_data(points, faces, is_one_indexed); // create internal mesh data structure including edges
     geodesic::GeodesicAlgorithmExact algorithm(&mesh); // create exact algorithm for the mesh
     std::vector <unsigned> rows_vector, columns_vector;
     std::vector <double> data_vector;
@@ -108,7 +110,8 @@ extern "C" {
         unsigned *source_indices_array,
         unsigned *target_indices_array,
         double *distance,
-        double distance_limit
+        double distance_limit,
+        unsigned is_one_indexed
     ) {
         compute_gdist_impl(
             number_of_vertices,
@@ -120,7 +123,8 @@ extern "C" {
             source_indices_array,
             target_indices_array,
             distance,
-            distance_limit
+            distance_limit,
+            is_one_indexed
         );
     }
 
@@ -130,7 +134,8 @@ extern "C" {
         double *vertices,
         unsigned *triangles,
         unsigned *sparse_matrix_size,
-        double max_distance
+        double max_distance,
+        unsigned is_one_indexed
     ) {
         return local_gdist_matrix_impl(
             number_of_vertices,
@@ -138,7 +143,8 @@ extern "C" {
             vertices,
             triangles,
             sparse_matrix_size,
-            max_distance
+            max_distance,
+            is_one_indexed
         );
     }
 
diff --git a/gdist_c_api.h b/gdist_c_api.h
index dfd2a7f..717d24e 100644
--- a/gdist_c_api.h
+++ b/gdist_c_api.h
@@ -26,7 +26,8 @@ void compute_gdist_impl(
     unsigned *source_indices_array,
     unsigned *target_indices_array,
     double *distance,
-    double distance_limit
+    double distance_limit,
+    unsigned is_one_indexed
 );
 
 double* local_gdist_matrix_impl(
@@ -35,7 +36,8 @@ double* local_gdist_matrix_impl(
     double *vertices,
     unsigned *triangles,
     unsigned *sparse_matrix_size,
-    double max_distance
+    double max_distance,
+    unsigned is_one_indexed
 );
 
 void free_memory_impl(double *ptr);
@@ -51,7 +53,8 @@ extern "C" {
         unsigned *source_indices_array,
         unsigned *target_indices_array,
         double *distance,
-        double distance_limit
+        double distance_limit,
+        unsigned is_one_indexed
     );
 
     DLL_EXPORT_API double* local_gdist_matrix(
@@ -60,7 +63,8 @@ extern "C" {
         double *vertices,
         unsigned *triangles,
         unsigned *sparse_matrix_size,
-        double max_distance
+        double max_distance,
+        unsigned is_one_indexed
     );
 
     DLL_EXPORT_API void free_memory(double *ptr);
diff --git a/geodesic_library/geodesic_mesh.h b/geodesic_library/geodesic_mesh.h
index a1a7a60..bf895b2 100644
--- a/geodesic_library/geodesic_mesh.h
+++ b/geodesic_library/geodesic_mesh.h
@@ -34,10 +34,11 @@ class Mesh
     void initialize_mesh_data(unsigned num_vertices,
                               Points& p, 
                               unsigned num_faces,
-                              Faces& tri);    //build mesh from regular point-triangle representation
+                              Faces& tri,
+                              bool is_one_indexed=false);    //build mesh from regular point-triangle representation
 
     template<class Points, class Faces>
-    void initialize_mesh_data(Points& p, Faces& tri);    //build mesh from regular point-triangle representation
+    void initialize_mesh_data(Points& p, Faces& tri, bool is_one_indexed=false);    //build mesh from regular point-triangle representation
 
     std::vector<Vertex>& vertices(){return m_vertices;};
     std::vector<Edge>& edges(){return m_edges;};
@@ -111,21 +112,22 @@ inline unsigned Mesh::closest_vertices(SurfacePoint* p,
 }
 
 template<class Points, class Faces>
-void Mesh::initialize_mesh_data(Points& p, Faces& tri)    //build mesh from regular point-triangle representation
+void Mesh::initialize_mesh_data(Points& p, Faces& tri, bool is_one_indexed)    //build mesh from regular point-triangle representation
 {
     assert(p.size() % 3 == 0);
     unsigned const num_vertices = p.size() / 3;
     assert(tri.size() % 3 == 0);
     unsigned const num_faces = tri.size() / 3; 
 
-    initialize_mesh_data(num_vertices, p, num_faces, tri);
+    initialize_mesh_data(num_vertices, p, num_faces, tri, is_one_indexed);
 }
 
 template<class Points, class Faces>
 void Mesh::initialize_mesh_data(unsigned num_vertices,
                                 Points& p,
                                 unsigned num_faces,
-                                Faces& tri)
+                                Faces& tri,
+                                bool is_one_indexed)
 {
     unsigned const approximate_number_of_internal_pointers = (num_vertices + num_faces)*4;
     unsigned const max_number_of_pointer_blocks = 100; 
@@ -154,7 +156,7 @@ void Mesh::initialize_mesh_data(unsigned num_vertices,
         unsigned shift = 3*i;
         for(unsigned j=0; j<3; ++j)
         {
-            unsigned vertex_index = tri[shift + j];
+            unsigned vertex_index = tri[shift + j] - (unsigned)is_one_indexed;
             assert(vertex_index < num_vertices);
             f.adjacent_vertices()[j] = &m_vertices[vertex_index];
         }