diff --git a/Cargo.lock b/Cargo.lock
index 8adc702ba4de8ead2b2c0f086b667df02e12e0f0..6599b68f9a2920eae9ca5f5e10fbb06e25ae45c4 100644
--- a/Cargo.lock
+++ b/Cargo.lock
@@ -90,6 +90,7 @@ version = "0.1.0"
 dependencies = [
  "itertools",
  "nauty-Traces-sys",
+ "union-find",
 ]
 
 [[package]]
@@ -267,6 +268,12 @@ version = "1.0.14"
 source = "registry+https://github.com/rust-lang/crates.io-index"
 checksum = "adb9e6ca4f869e1180728b7950e35922a7fc6397f7b641499e8f3ef06e50dc83"
 
+[[package]]
+name = "union-find"
+version = "0.4.3"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "039142448432983c34b64739f8526f8f233a1eec7a66e61b6ab29acfa781194e"
+
 [[package]]
 name = "windows-targets"
 version = "0.52.6"
diff --git a/Cargo.toml b/Cargo.toml
index 49a2a0d985ec4b6948fbf5c4e8cfa674c38d4988..b41ed5c6ee55094cec9b93d747068c0ba14c7afb 100644
--- a/Cargo.toml
+++ b/Cargo.toml
@@ -6,3 +6,4 @@ edition = "2021"
 [dependencies]
 itertools = "0.13.0"
 nauty-Traces-sys = "0.8.0"
+union-find = "0.4.3"
diff --git a/src/main.rs b/src/main.rs
index 63812cbf162e9d750fb642b3b26487ea140b422a..3d3e8ef1738cb6a74f12e18620a82e532bd076b6 100644
--- a/src/main.rs
+++ b/src/main.rs
@@ -6,9 +6,14 @@ fn main() {
     println!("g = {:?}", g);
     println!("g.to_nauty() = {:?}", g.to_nauty());
     println!("UndirectedGraph::from_nauty(&g.to_nauty()) = {:?}", UndirectedGraph::from_nauty(&g.to_nauty()));
-    let (_sign, normal_g, orbits) = g.normal_form();
+    let (_sign, normal_g, orbits,_) = g.normal_form(false);
     println!("normal_g = {:?}", normal_g);
-    println!("normal_g.expanding_differential(orbits) = {:?}", normal_g.expanding_differential(&orbits));
+    println!("normal_g.expanding_differential(&orbits) = {:?}", normal_g.expanding_differential(&orbits));
+    println!("normal_g.contractions() = {:?}", normal_g.contractions(&vec![0_usize]));
+    let h: UndirectedGraph = UndirectedGraph::new(7,vec![(0, 1), (0, 2), (0, 6), (1, 4), (1, 6), (2, 5), (2, 6), (3, 4), (3, 5), (3, 6), (4, 5)]);
+    let edge_orbit_reps_h: Vec<usize> = vec![0, 2, 3, 4, 7, 9, 10];
+    println!("normal_h.contractions(&edge_orbits_h) = {:?}", h.contractions(&edge_orbit_reps_h));
     let triangle: UndirectedGraph = UndirectedGraph::new(3, vec![(0,1), (1,2), (0,2)]);
-    println!("triangle.normal_form() = {:?}", triangle.normal_form());
+    println!("triangle.normal_form() = {:?}", triangle.normal_form(true));
+    println!("triangle.contractions() = {:?}", triangle.contractions(&vec![0_usize]));
 }
diff --git a/src/undirected_graph.rs b/src/undirected_graph.rs
index 5194186a642c618542ba29184f7056bb9c77609b..0aef71cc458c1e9fc0a7c70606ce4412e485001c 100644
--- a/src/undirected_graph.rs
+++ b/src/undirected_graph.rs
@@ -120,7 +120,10 @@ impl UndirectedGraph {
         }
     }
 
-    pub fn normal_form(&self) -> (i32, UndirectedGraph, HashMap<usize, usize>) {
+    pub fn normal_form(
+        &self,
+        compute_edge_orbits: bool,
+    ) -> (i32, UndirectedGraph, HashMap<usize, usize>, Vec<usize>) {
         let mut nauty_graph: Vec<u64> = self.to_nauty();
 
         let mut options = nauty_Traces_sys::optionblk::default();
@@ -172,27 +175,60 @@ impl UndirectedGraph {
             }
         }
 
-        let mut is_odd: bool = false;
+        let mut sign: i32 = 1;
+        let mut normal_self_edge_orbits: Vec<usize> = vec![];
         Self::AUTOMORPHSM_GROUP_GENERATORS.with_borrow(|gens| {
+            // NOTE: Use Union-Find to compute edge orbits.
+            use union_find::{QuickUnionUf, UnionByRank, UnionFind};
+            let mut uf: QuickUnionUf<UnionByRank> =
+                QuickUnionUf::<UnionByRank>::new(self.edges.len());
+
             for automorphism in gens {
                 // NOTE: The automorphisms produced by nauty act on the original graph.
                 let sigma: Vec<i32> = Self::induced_edge_permutation(&self, &self, automorphism);
+
+                if compute_edge_orbits {
+                    // Take powers of sigma and act on edges of self.
+                    let mut tau: Vec<i32> = sigma.clone();
+                    let num_edges: i32 = self.edges.len() as i32;
+                    let identity_permutation: Vec<i32> = (0..num_edges).collect();
+                    while tau != identity_permutation {
+                        for k in 0..self.edges.len() {
+                            uf.union(k, tau[k] as usize);
+                        }
+                        for e in tau.iter_mut() {
+                            *e = sigma[*e as usize];
+                        }
+                    }
+                }
+
                 if Self::permutation_sign(sigma) == -1 {
-                    is_odd = true;
+                    sign = 0;
+                    // NOTE: Giving up computation of edge orbits; we will not use them anyway.
                     break;
                 }
             }
-        });
 
-        let sign: i32 = if is_odd {
-            0
-        } else {
-            let edge_permutation: Vec<i32> =
-                Self::induced_edge_permutation(&self, &normal_self, &lab);
-            Self::permutation_sign(edge_permutation)
-        };
+            if sign != 0 {
+                let edge_permutation: Vec<i32> =
+                    Self::induced_edge_permutation(&normal_self, &self, &lab_inv);
+                if compute_edge_orbits {
+                    normal_self_edge_orbits = (0..self.edges.len()).map(|e| uf.find(e)).collect();
+                    // So far these are the edge orbits of self. Finally we apply the isomorphism to get the edge orbits of normal_self.
+                    for e in &mut normal_self_edge_orbits {
+                        *e = edge_permutation[*e] as usize;
+                    }
+                }
+                sign = Self::permutation_sign(edge_permutation);
+            }
+        });
 
-        return (sign, normal_self, normal_self_orbits);
+        return (
+            sign,
+            normal_self,
+            normal_self_orbits,
+            normal_self_edge_orbits,
+        );
     }
 
     fn expand_vertex(&self, &position: &usize) -> Vec<(usize, UndirectedGraph)> {
@@ -260,27 +296,28 @@ impl UndirectedGraph {
     pub fn expanding_differential(
         &self,
         orbits: &HashMap<usize, usize>,
-    ) -> HashMap<UndirectedGraph, i64> {
-        let mut diff: HashMap<UndirectedGraph, i64> = HashMap::new();
+    ) -> HashMap<UndirectedGraph, (i64, Vec<usize>)> {
+        let mut diff: HashMap<UndirectedGraph, (i64, Vec<usize>)> = HashMap::new();
 
         for (&orbit_rep, &orbit_size) in orbits {
             for (c, dg) in self.expand_vertex(&orbit_rep) {
                 // TODO: Avoid computing orbits when not used.
-                let (normal_dg_sign, normal_dg, _normal_dg_orbits) = dg.normal_form();
+                let (normal_dg_sign, normal_dg, _normal_dg_orbits, normal_dg_edge_orbits) =
+                    dg.normal_form(true);
                 if normal_dg_sign == 0 {
                     continue;
                 }
                 // TODO: Fix coefficient types.
                 let normal_dg_coeff: i64 =
                     (normal_dg_sign as i64) * (orbit_size as i64) * (c as i64);
-                if let Some(diff_coeff) = diff.get_mut(&normal_dg) {
+                if let Some((diff_coeff, _)) = diff.get_mut(&normal_dg) {
                     *diff_coeff += normal_dg_coeff;
                 } else {
-                    diff.insert(normal_dg, normal_dg_coeff);
+                    diff.insert(normal_dg, (normal_dg_coeff, normal_dg_edge_orbits));
                 }
             }
         }
-        diff.retain(|_dg, dg_coeff| *dg_coeff != 0);
+        diff.retain(|_dg, (dg_coeff, _)| *dg_coeff != 0);
 
         return diff;
     }
@@ -323,7 +360,7 @@ impl UndirectedGraph {
                 .all(|new_edge| unique_edges.insert(*new_edge));
             if has_no_double_edges {
                 let cg: UndirectedGraph = UndirectedGraph::new(self.num_vertices - 1, new_edges);
-                let (_normal_cg_sign, normal_cg, normal_cg_orbits) = cg.normal_form();
+                let (_normal_cg_sign, normal_cg, normal_cg_orbits, _) = cg.normal_form(false);
                 contractions.insert(normal_cg, normal_cg_orbits);
             }
         }