diff --git a/dom/dom.py b/dom/dom.py
index 24d36b4..42be748 100644
--- a/dom/dom.py
+++ b/dom/dom.py
@@ -179,13 +179,15 @@ def sharpness_matrix(self, Im, width=2, debug=False):
         for i in range(width, domx.shape[0]-width):
             num = np.abs(domx[i-width:i+width, :]).sum(axis=0)
             dn = Cx[i-width:i+width, :].sum(axis=0)
-            Sx[i] = [(num[k]/dn[k] if dn[k] > 1e-3 else 0) for k in range(Sx.shape[1])]
+            # Sx[i] = [(num[k]/dn[k] if dn[k] > 1e-3 else 0) for k in range(Sx.shape[1])]    # This is quite slow
+            Sx[i] = np.where(dn > 1e-3, num/dn, 0)     # Much faster
         
         # Compute Sy
         for j in range(width, domy.shape[1]-width):
             num = np.abs(domy[:, j-width: j+width]).sum(axis=1)
             dn = Cy[:, j-width:j+width].sum(axis=1)
-            Sy[:, j] = [(num[k]/dn[k] if dn[k] > 1e-3 else 0) for k in range(Sy.shape[0])]
+            # Sy[:, j] = [(num[k]/dn[k] if dn[k] > 1e-3 else 0) for k in range(Sy.shape[0])]
+            Sy[:, j] = np.where(dn > 1e-3, num/dn, 0)
             
         if debug:
             print(f"domx {domx.shape}: {[(i,round(np.quantile(domx, i/100), 2)) for i in range(0, 101, 25)]}")