1 files changed, 92 insertions, 3 deletions
diff --git a/libquanta_vquantise_wu.c b/libquanta_vquantise_wu.c
index 6e1ce94..55129f2 100644
--- a/libquanta_vquantise_wu.c
+++ b/libquanta_vquantise_wu.c
@@ -1,4 +1,4 @@
-s/* See LICENSE file for copyright and license details. */
+/* See LICENSE file for copyright and license details. */
 #include "common.h"
 
 /* TODO limitation: 32 values per channel, limits palette size */
@@ -77,8 +77,96 @@ compute_histogram(struct libquanta_image *image, struct histogram *histogram, st
 
 
 static void
-compute_moments(struct histogram *histogram, struct channel_data *channels, size_t nchannels) /* TODO */
+moments_bigint_(size_t channel, const struct channel_data *channels, size_t nchannels,
+                struct bigint *out, const struct bigint *in, size_t index)
 {
+	size_t i, mul;
+
+	if (!nchannels) {
+		if (in != out)
+			memcpy(&out[index], &in[index], sizeof(*in));
+		bigint_add_big(&out[index], &out[index - channels[channel].subindex_multiplier]);
+		return;
+	}
+
+	mul = channels[--nchannels].subindex_multiplier;
+	for (i = 0; i++ < (1U << WORKING_BITS);)
+		moments_bigint_(channel, channels, nchannels, out, in, index + i * mul);
+}
+
+static void
+moments_bigint(size_t channel, struct channel_data *channels, size_t nchannels, struct bigint *buf, size_t bufsize)
+{
+	struct bigint *in = channels[channel].histogram;
+	size_t i;
+
+	memset(buf, 0, bufsize);
+
+	for (i = nchannels; i-- > 1U;) {
+		moments_bigint_(i, channels, nchannels, buf, in, 0);
+		in = buf;
+	}
+
+	moments_bigint_(0, channels, nchannels, channels[channel].histogram, in, 0);
+}
+
+static void
+moments_histogram_(size_t channel, const struct channel_data *channels, size_t nchannels,
+                   struct histogram *out, const struct histogram *in, size_t index)
+{
+	size_t i, mul;
+
+	if (!nchannels) {
+		size_t prev = index - channels[channel].subindex_multiplier;
+		out[index].occurrences = out[prev].occurrences + in[index].occurrences;
+		out[index].square_sum = out[prev].square_sum + in[index].square_sum;
+		return;
+	}
+
+	mul = channels[--nchannels].subindex_multiplier;
+	for (i = 0; i++ < (1U << WORKING_BITS);)
+		moments_histogram_(channel, channels, nchannels, out, in, index + i * mul);
+}
+
+static void
+moments_histogram(struct histogram *histogram, const struct channel_data *channels, size_t nchannels,
+                  struct histogram *buf, size_t bufsize)
+{
+	struct histogram *in = histogram;
+	size_t i;
+
+	memset(buf, 0, bufsize);
+
+	for (i = nchannels; i-- > 1U;) {
+		moments_histogram_(i, channels, nchannels, buf, in, 0);
+		in = buf;
+	}
+
+	moments_histogram_(0, channels, nchannels, histogram, in, 0);
+}
+
+static int
+compute_cumulative_moments(struct histogram *histogram, struct channel_data *channels, size_t nchannels, size_t bufelems)
+{
+	void *buf;
+	size_t i, bufsize;
+
+	bufsize = bufelems * sizeof(struct bigint);
+	buf = malloc(bufsize);
+	if (!buf)
+		return -1;
+	for (i = 0; i < nchannels; i++)
+		moments_bigint(i, channels, nchannels, buf, bufsize);
+	free(buf);
+
+	bufsize = bufelems * sizeof(struct histogram);
+	buf = malloc(bufsize);
+	if (!buf)
+		return -1;
+	moments_histogram(histogram, channels, nchannels, buf, bufsize);
+	free(buf);
+
+	return 0;
 }
 
 
@@ -661,7 +749,8 @@ libquanta_vquantise_wu(struct libquanta_palette *palette, va_list args)
 
 	/* analyse image */
 	compute_histogram(ret, histogram, channels, nchannels);
-	compute_moments(histogram, channels, nchannels);
+	if (compute_cumulative_moments(histogram, channels, nchannels, subindex_multiplier))
+		goto fail;
 	if (partition(cubes, histogram, channels, nchannels, ncolours, &ncolours_used))
 		goto fail;