diff --git a/src/moc/range/mod.rs b/src/moc/range/mod.rs
index d36d383..80cc9aa 100644
--- a/src/moc/range/mod.rs
+++ b/src/moc/range/mod.rs
@@ -12,6 +12,7 @@ use std::{
   vec::IntoIter,
 };
 
+use healpix::compass_point::{Cardinal, CardinalSet};
 /// Re-export `Ordinal` not to be out-of-sync with cdshealpix version.
 pub use healpix::compass_point::{Ordinal, OrdinalMap, OrdinalSet};
 use healpix::{
@@ -970,6 +971,77 @@ impl<T: Idx> RangeMOC<T, Hpx<T>> {
         }
       })
   }
+
+  /// Same as `border_elementary_edges` but return edges vertices in an iterator of arrays of the form:
+  /// `[vertex_1.lon, vertex_1.lat, vertex_2_.lon, vertex_2.lat]`
+  pub fn border_elementary_edges_vertices(&self) -> impl Iterator<Item = [f64; 4]> + '_ {
+    let hp = healpix::nested::get(self.depth_max);
+    self
+      .internal_border_iter()
+      .flatten_to_fixed_depth_cells()
+      .filter_map(move |idx| {
+        let mut edges = OrdinalSet::new();
+        let mut vertices_set = CardinalSet::new();
+        let h = <T as Idx>::to_u64(idx);
+        let neigs = hp.neighbours(h, false);
+        if let Some(nh) = neigs.get(MainWind::NE) {
+          if !self.contains_depth_max_val(&T::from_u64(*nh)) {
+            edges.set(Ordinal::NE, true);
+            vertices_set.set(Cardinal::N, true);
+            vertices_set.set(Cardinal::E, true);
+          }
+        }
+        if let Some(nh) = neigs.get(MainWind::SE) {
+          if !self.contains_depth_max_val(&T::from_u64(*nh)) {
+            edges.set(Ordinal::SE, true);
+            vertices_set.set(Cardinal::S, true);
+            vertices_set.set(Cardinal::E, true);
+          }
+        }
+        if let Some(nh) = neigs.get(MainWind::NW) {
+          if !self.contains_depth_max_val(&T::from_u64(*nh)) {
+            edges.set(Ordinal::NW, true);
+            vertices_set.set(Cardinal::N, true);
+            vertices_set.set(Cardinal::W, true);
+          }
+        }
+        if let Some(nh) = neigs.get(MainWind::SW) {
+          if !self.contains_depth_max_val(&T::from_u64(*nh)) {
+            edges.set(Ordinal::SW, true);
+            vertices_set.set(Cardinal::S, true);
+            vertices_set.set(Cardinal::W, true);
+          }
+        }
+        if edges.is_empty() {
+          None
+        } else {
+          let vertices = hp.vertices_map(h, vertices_set);
+          Some(edges.into_iter().map(move |ordinal| match ordinal {
+            Ordinal::NE => {
+              let n = vertices.get(Cardinal::N).unwrap();
+              let e = vertices.get(Cardinal::E).unwrap();
+              [n.0, n.1, e.0, e.1]
+            }
+            Ordinal::SE => {
+              let s = vertices.get(Cardinal::S).unwrap();
+              let e = vertices.get(Cardinal::E).unwrap();
+              [s.0, s.1, e.0, e.1]
+            }
+            Ordinal::NW => {
+              let n = vertices.get(Cardinal::N).unwrap();
+              let w = vertices.get(Cardinal::W).unwrap();
+              [n.0, n.1, w.0, w.1]
+            }
+            Ordinal::SW => {
+              let s = vertices.get(Cardinal::S).unwrap();
+              let w = vertices.get(Cardinal::W).unwrap();
+              [s.0, s.1, w.0, w.1]
+            }
+          }))
+        }
+      })
+      .flatten()
+  }
 }
 
 fn from<T: Idx + TryFrom<u64, Error = TryFromIntError>>(
diff --git a/src/storage/u64idx/mod.rs b/src/storage/u64idx/mod.rs
index 05265c3..46341f9 100644
--- a/src/storage/u64idx/mod.rs
+++ b/src/storage/u64idx/mod.rs
@@ -67,8 +67,8 @@ use self::{
     tmoc_from_fits_gen,
   },
   op1::{
-    op1_1st_axis_max, op1_1st_axis_min, op1_count_split, op1_flatten_to_depth,
-    op1_flatten_to_moc_depth, op1_moc_barycenter,
+    op1_1st_axis_max, op1_1st_axis_min, op1_border_elementary_edges_vertices, op1_count_split,
+    op1_flatten_to_depth, op1_flatten_to_moc_depth, op1_moc_barycenter,
     op1_moc_largest_distance_from_coo_to_moc_vertices, Op1, Op1MultiRes,
   },
   op2::Op2,
@@ -1793,6 +1793,14 @@ impl U64MocStore {
     op1_flatten_to_depth(index, depth)
   }
 
+  /// Returns the MOC elementary edges, i.e. the edges at the deepest depth, in any order, made of
+  /// the starting vertex and the ending vertex so that each item is of the form:
+  /// `[stat_vertex_lon, starting_vertex_lat, ending_vertex_lon, ending_vertex_lat]`.
+  pub fn border_elementary_edges_vertices(&self, index: usize) -> Result<Vec<[f64; 4]>, String> {
+    op1_border_elementary_edges_vertices(index)
+  }
+  //
+
   /// Split the given disjoint S-MOC int joint S-MOCs.
   /// Split "direct", i.e. we consider 2 neighboring cells to be the same only if the share an edge.
   /// WARNING: may create a lot of new MOCs, exec `splitCount` first!!
diff --git a/src/storage/u64idx/op1.rs b/src/storage/u64idx/op1.rs
index 86ef4ed..d230835 100644
--- a/src/storage/u64idx/op1.rs
+++ b/src/storage/u64idx/op1.rs
@@ -222,6 +222,21 @@ pub(crate) fn op1_flatten_to_depth(index: usize, depth: u8) -> Result<Vec<u64>,
   })
 }
 
+pub(crate) fn op1_border_elementary_edges_vertices(index: usize) -> Result<Vec<[f64; 4]>, String> {
+  store::exec_on_one_readonly_moc(index, move |moc| match moc {
+    InternalMoc::Space(m) => Ok(m.border_elementary_edges_vertices().collect()),
+    InternalMoc::Time(_) => Err(String::from(
+      "Border elementary edges vertices not implemented for T-MOCs.",
+    )),
+    InternalMoc::Frequency(_) => Err(String::from(
+      "Border elementary edges vertices not implemented for F-MOCs.",
+    )),
+    InternalMoc::TimeSpace(_) => Err(String::from(
+      "Border elementary edges vertices not implemented for ST-MOCs.",
+    )),
+  })
+}
+
 pub(crate) fn op1_count_split(index: usize, indirect_neigh: bool) -> Result<u32, String> {
   store::exec_on_one_readonly_moc(index, move |moc| match moc {
     InternalMoc::Space(m) => Ok(m.split_into_joint_mocs(indirect_neigh).len() as u32),