diff --git a/skimage/color/delta_e.py b/skimage/color/delta_e.py
index 43b4a7b2..a7f9b775 100644
--- a/skimage/color/delta_e.py
+++ b/skimage/color/delta_e.py
@@ -108,7 +108,7 @@ def deltaE_ciede94(lab1, lab2, kH=1, kC=1, kL=1, k1=0.045, k2=0.015):
 
     dL = L1 - L2
     dC = C1 - C2
-    dH = get_dH(lab1, lab2)
+    dH2 = get_dH2(lab1, lab2)
 
     SL = 1
     SC = 1 + k1 * C1
@@ -116,7 +116,7 @@ def deltaE_ciede94(lab1, lab2, kH=1, kC=1, kL=1, k1=0.045, k2=0.015):
 
     dE2 = (dL / (kL * SL)) ** 2
     dE2 += (dC / (kC * SC)) ** 2
-    dE2 += (dH / (kH * SH)) ** 2
+    dE2 += dH2 / (kH * SH) ** 2
     return np.sqrt(dE2)
 
 
@@ -289,7 +289,7 @@ def deltaE_cmc(lab1, lab2, kL=1, kC=1):
 
     dC = C1 - C2
     dL = L1 - L2
-    dH = get_dH(lab1, lab2)
+    dH2 = get_dH2(lab1, lab2)
 
     T = np.where(np.logical_and(np.rad2deg(h1) >= 164, np.rad2deg(h1) <= 345),
                  0.56 + 0.2 * np.abs(np.cos(h1 + np.deg2rad(168))),
@@ -304,29 +304,36 @@ def deltaE_cmc(lab1, lab2, kL=1, kC=1):
 
     dE2 = (dL / (kL * SL)) ** 2
     dE2 += (dC / (kC * SC)) ** 2
-    dE2 += (dH / SH) ** 2
+    dE2 += dH2 / (SH ** 2)
     return np.sqrt(dE2)
 
 
-def get_dH(lab1, lab2):
-    """numerically well behaved calculation of dH
+def get_dH2(lab1, lab2):
+    """squared hue difference term occurring in deltaE_cmc and deltaE_ciede94
 
-    Given a1, b1, a2, b2
-    c1 = sqrt(a1**2 + b1**2)
-    c2 = sqrt(a2**2 + b2**2)
-    term = (a1-a2)**2 + (b1-b2)**2 - (c1-c2)**2
-    dH = sqrt(term)
+    Despite its name "dH" is not a simple difference of hue values.  We avoid
+    working directly with the hue value directly since differencing angles is
+    troublesome.  The hue term is usually written as:
+        c1 = sqrt(a1**2 + b1**2)
+        c2 = sqrt(a2**2 + b2**2)
+        term = (a1-a2)**2 + (b1-b2)**2 - (c1-c2)**2
+        dH = sqrt(term)
+
+    However, this has poor roundoff properties when a or b is dominant.
+    Instead, r is a vector with elements a and b.  The same dH term can be
+    re-written as:
+        |r1-r2|**2 - (|r1| - |r2|)**2
+    and then simplified to:
+        2*|r1|*|r2| - 2*dot(r1, r2)
     """
     lab1 = np.asarray(lab1)
     lab2 = np.asarray(lab2)
+    a1, b1 = np.rollaxis(lab1, -1)[1:3]
+    a2, b2 = np.rollaxis(lab2, -1)[1:3]
 
-    r1 = np.rollaxis(lab1, -1)[1:3]
-    r2 = np.rollaxis(lab2, -1)[1:3]
+    # magnitude of (a, b) is the chroma
+    C1 = np.hypot(a1, b1)
+    C2 = np.hypot(a2, b2)
 
-    C1 = np.hypot(r1[0], r1[1])
-    C2 = np.hypot(r2[0], r2[1])
-
-    term = C1 * C2
-    term -= (r1 * r2).sum(0)
-    term *= 2
-    return np.sqrt(term)
+    term = (C1 * C2) - (a1 * a2 + b1 * b2)
+    return 2*term