use circ mean for hue

ingest2.py

@@ -99,10 +99,13 @@ def calc_statistics(pixels: np.array) -> Stats:
   tilt = (pixels / np.sqrt(sqnorms)[:, np.newaxis]).mean(axis=0)
   # variance = mean(||p||^2) - ||mean(p)||^2
   variance = sqnorms.mean(axis=0) - sum(mean ** 2)
-  # chroma^2 = a^2 + b^2, C-bar = mean(sqrt(a^2 + b^2))
-  chroma = np.sqrt(squared[:, 1:].sum(axis=1)).mean(axis=0)
-  # hue = atan2(b, a), h-bar = mean(atan2(b, a))
-  hue = np.arctan2(pixels[:, 2], pixels[:, 1]).mean(axis=0) * 180 / math.pi
+  # chroma^2 = a^2 + b^2
+  chroma = np.sqrt(squared[:, 1:].sum(axis=1))
+  # hue = atan2(b, a), but we need a circular mean
+  # https://en.wikipedia.org/wiki/Circular_mean#Definition
+  # cos(atan2(b, a)) = a / sqrt(a^2 + b^2) = a / chroma
+  # sin(atan2(b, a)) = b / sqrt(a^2 + b^2) = b / chroma
+  hue = math.atan2(*(pixels[:, [2, 1]] / chroma[:, np.newaxis]).mean(axis=0))
   return Stats(
     size=len(pixels),
     variance=variance,
@@ -111,8 +114,8 @@ def calc_statistics(pixels: np.array) -> Stats:
     Lbar=mean[0],
     abar=mean[1],
     bbar=mean[2],
-    Cbar=chroma,
-    hbar=hue,
+    Cbar=chroma.mean(axis=0),
+    hbar=hue * 180 / math.pi,
     Lhat=tilt[0],
     ahat=tilt[1],
     bhat=tilt[2],