First commit

Signed-off-by: Mattias Andrée <m@maandree.se>

30 files changed, 1863 insertions, 0 deletions

diff --git a/.gitignore b/.gitignore
new file mode 100644
index 0000000..a071ed4
--- /dev/null
+++ b/.gitignore
@@ -0,0 +1,14 @@
+*\#*
+*~
+*.o
+*.a
+*.lo
+*.su
+*.so
+*.so.*
+*.dll
+*.dylib
+*.gch
+*.gcov
+*.gcno
+*.gcda
diff --git a/LICENSE b/LICENSE
new file mode 100644
index 0000000..0e6be1c
--- /dev/null
+++ b/LICENSE
@@ -0,0 +1,15 @@
+ISC License
+
+© 2025 Mattias Andrée <m@maandree.se>
+
+Permission to use, copy, modify, and/or distribute this software for any
+purpose with or without fee is hereby granted, provided that the above
+copyright notice and this permission notice appear in all copies.
+
+THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES
+WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF
+MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR
+ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
+WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN
+ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF
+OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.
diff --git a/Makefile b/Makefile
new file mode 100644
index 0000000..26127d3
--- /dev/null
+++ b/Makefile
@@ -0,0 +1,99 @@
+.POSIX:
+
+CONFIGFILE = config.mk
+include $(CONFIGFILE)
+
+OS = linux
+# Linux:   linux
+# Mac OS:  macos
+# Windows: windows
+include mk/$(OS).mk
+
+
+LIB_MAJOR = 1
+LIB_MINOR = 0
+LIB_VERSION = $(LIB_MAJOR).$(LIB_MINOR)
+LIB_NAME = quanta
+
+
+OBJ_PUBLIC =\
+	libquanta_palette_size.o\
+	libquanta_malloc_palette.o\
+	libquanta_vquantise.o\
+	libquanta_quantise.o\
+	libquanta_vquantise_wu.o\
+	libquanta_quantise_wu.o
+
+OBJ_PRIVATE =\
+	libquanta_bigint_divmod_small__.o
+
+OBJ =\
+	$(OBJ_PUBLIC)\
+	$(OBJ_PRIVATE)
+
+HDR =\
+	libquanta.h\
+	common.h
+
+LOBJ = $(OBJ:.o=.lo)
+
+MAN3 =\
+	$(OBJ_PUBLIC:.o=.3)\
+	struct_libquanta_image.3 libquanta_image.3\
+	struct_libquanta_channel.3 libquanta_channel.3\
+	struct_libquanta_palette.3 libquanta_palette.3\
+
+MAN7 = libquanta.7
+
+
+all: libquanta.a libquanta.$(LIBEXT)
+$(OBJ): $(HDR)
+$(LOBJ): $(HDR)
+
+.c.o:
+	$(CC) -c -o $@ $< $(CFLAGS) $(CPPFLAGS)
+
+.c.lo:
+	$(CC) -fPIC -c -o $@ $< $(CFLAGS) $(CPPFLAGS)
+
+libquanta.a: $(OBJ)
+	@rm -f -- $@
+	$(AR) rc $@ $(OBJ)
+	$(AR) ts $@ > /dev/null
+
+libquanta.$(LIBEXT): $(LOBJ)
+	$(CC) $(LIBFLAGS) -o $@ $(LOBJ) $(LDFLAGS)
+
+install: libquanta.a libquanta.$(LIBEXT)
+	mkdir -p -- "$(DESTDIR)$(PREFIX)/lib"
+	mkdir -p -- "$(DESTDIR)$(PREFIX)/include"
+	mkdir -p -- "$(DESTDIR)$(MANPREFIX)/man3"
+	mkdir -p -- "$(DESTDIR)$(MANPREFIX)/man7"
+	cp -- libquanta.a "$(DESTDIR)$(PREFIX)/lib/"
+	cp -- libquanta.$(LIBEXT) "$(DESTDIR)$(PREFIX)/lib/libquanta.$(LIBMINOREXT)"
+	$(FIX_INSTALL_NAME) "$(DESTDIR)$(PREFIX)/lib/libquanta.$(LIBMINOREXT)"
+	ln -sf -- libquanta.$(LIBMINOREXT) "$(DESTDIR)$(PREFIX)/lib/libquanta.$(LIBMAJOREXT)"
+	ln -sf -- libquanta.$(LIBMAJOREXT) "$(DESTDIR)$(PREFIX)/lib/libquanta.$(LIBEXT)"
+	cp -- libquanta.h "$(DESTDIR)$(PREFIX)/include/"
+	cp -P -- $(MAN3) "$(DESTDIR)$(MANPREFIX)/man3"
+	cp -P -- $(MAN7) "$(DESTDIR)$(MANPREFIX)/man7"
+
+uninstall:
+	-rm -f -- "$(DESTDIR)$(PREFIX)/lib/libquanta.a"
+	-rm -f -- "$(DESTDIR)$(PREFIX)/lib/libquanta.$(LIBMAJOREXT)"
+	-rm -f -- "$(DESTDIR)$(PREFIX)/lib/libquanta.$(LIBMINOREXT)"
+	-rm -f -- "$(DESTDIR)$(PREFIX)/lib/libquanta.$(LIBEXT)"
+	-rm -f -- "$(DESTDIR)$(PREFIX)/include/libquanta.h"
+	-cd -- "$(DESTDIR)$(MANPREFIX)/man3" && rm -f -- $(MAN3)
+	-cd -- "$(DESTDIR)$(MANPREFIX)/man7" && rm -f -- $(MAN7)
+	-rmdir -- "$(DESTDIR)$(MANPREFIX)/man3"
+	-rmdir -- "$(DESTDIR)$(MANPREFIX)/man7"
+
+clean:
+	-rm -f -- *.o *.a *.lo *.su *.so *.so.* *.dll *.dylib
+	-rm -f -- *.gch *.gcov *.gcno *.gcda *.$(LIBEXT)
+
+.SUFFIXES:
+.SUFFIXES: .lo .o .c
+
+.PHONY: all install uninstall clean
diff --git a/README b/README
new file mode 100644
index 0000000..f0e4f1e
--- /dev/null
+++ b/README
@@ -0,0 +1,44 @@
+NAME
+	libquanta - Colour quantisation library
+
+SYNOPSIS
+	#include <libquanta.h>
+
+	Link with -lquanta.
+
+DESCRIPTION
+	The libquanta software library is a C library for colour
+	quantisation of images, meaning it lets you create a colour
+	palette for an image and convert the image to a palette-based
+	image.
+
+	The palette may have lower colour resolution and fewer colours
+	than the original image.
+
+	libquanta is not designed to convert an image to a
+	palette-based image based on a pre-existing palette. It is
+	only designed to create a palette suitable for a provided
+	image.
+
+	The libquanta software library implements the following
+	functions:
+
+   Colour palette functions
+	libquanta_palette_size(3)
+		Calculate the size of a colour palette.
+
+	libquanta_malloc_palette(3)
+		Allocate a colour palette.
+
+   Colour quantisation functions
+	libquanta_quantise(3),
+	libquanta_vquantise(3)
+		Colour-quantise an image using an algorithm selected
+		by the library.
+
+	libquanta_quantise_wu(3),
+	libquanta_vquantise_wu(3)
+		Colour-quantise an image using Wu's Colour Quantiser.
+
+SEE ALSO
+	None.
diff --git a/common.h b/common.h
new file mode 100644
index 0000000..b299fd1
--- /dev/null
+++ b/common.h
@@ -0,0 +1,91 @@
+/* See LICENSE file for copyright and license details. */
+#include "libquanta.h"
+#include <errno.h>
+#include <limits.h>
+#include <math.h>
+#include <stdlib.h>
+
+
+#define PALETTE_BASE_SIZE       offsetof(struct libquanta_palette, palette)
+#define PALETTE_VALUE_TYPE      uint64_t
+#define PALETTE_VALUE_SIZE      sizeof(PALETTE_VALUE_TYPE)
+#define PALETTE_VALUE_MAX_BITS  64U
+
+
+struct bigint {
+	uintmax_t high, low;
+};
+
+
+#if defined(__GNUC__)
+__attribute__((__visibility__("hidden")))
+#endif
+uintmax_t libquanta_bigint_divmod_small__(struct bigint *big, uintmax_t small);
+#define bigint_divmod_small libquanta_bigint_divmod_small__
+
+
+static inline void
+bigint_add_small(struct bigint *big, uintmax_t small)
+{
+#if defined(__GNUC__)
+	if (__builtin_add_overflow(big->low, small, &big->low))
+		big->high += 1U;
+#else
+	if (big->low > UINTMAX_MAX - small)
+		big->high += 1U;
+	big->low += small;
+#endif
+}
+
+
+static inline void
+bigint_sub_small(struct bigint *big, uintmax_t small)
+{
+#if defined(__GNUC__)
+	if (__builtin_sub_overflow(big->low, small, &big->low))
+		big->high -= 1U;
+#else
+	if (big->low < small)
+		big->high -= 1U;
+	big->low -= small;
+#endif
+}
+
+
+static inline void
+bigint_add_big(struct bigint *res, const struct bigint *restrict other)
+{
+	bigint_add_small(res, other->low);
+	res->high += other->high;
+}
+
+
+static inline void
+bigint_sub_big(struct bigint *res, const struct bigint *restrict other)
+{
+	bigint_sub_small(res, other->low);
+	res->high -= other->high;
+}
+
+
+static inline void
+bigint_rsub_big(struct bigint *res, const struct bigint *restrict other)
+{
+	res->high = other->high - res->high;
+#if defined(__GNUC__)
+	if (__builtin_sub_overflow(other->low, res->low, &res->low))
+		res->high -= 1U;
+#else
+	if (res->low > other->low)
+		res->high -= 1U;
+	res->low = other->low - res->low;
+#endif
+}
+
+
+static inline uintmax_t
+bigint_div_small(const struct bigint *big, uintmax_t small)
+{
+	struct bigint big_copy = *big;
+	return bigint_divmod_small(&big_copy, small);
+}
diff --git a/config.mk b/config.mk
new file mode 100644
index 0000000..037fe35
--- /dev/null
+++ b/config.mk
@@ -0,0 +1,8 @@
+PREFIX    = /usr
+MANPREFIX = $(PREFIX)/share/man
+
+CC = c99
+
+CPPFLAGS = -D_DEFAULT_SOURCE -D_BSD_SOURCE -D_XOPEN_SOURCE=700 -D_GNU_SOURCE
+CFLAGS   =
+LDFLAGS  = -lm
diff --git a/libquanta.7 b/libquanta.7
new file mode 100644
index 0000000..d68b73e
--- /dev/null
+++ b/libquanta.7
@@ -0,0 +1,55 @@
+.TH LIBQUANTA 7 libquanta
+.SH NAME
+libquanta \- Colour quantisation library
+
+.SH SYNPOSIS
+.nf
+#include <libquanta.h>
+.fi
+.PP
+Link with
+.I -lquanta
+.IR -lm .
+
+.SH DESCRIPTION
+The
+.B libquanta
+software library is a C library for colour quantisation
+of images, meaning it lets you create a colour palette for
+an image and convert the image to a palette-based image.
+.PP
+The palette may have lower colour resolution and fewer
+colours than the original image.
+.PP
+.B libquanta
+is not designed to convert an image to a palette-based
+image based on a pre-existing palette. It is only designed
+to create a palette suitable for a provided image.
+.PP
+The
+.B libquanta
+software library implements the following functions:
+
+.SS Colour palette functions
+.TP
+.BR libquanta_palette_size (3)
+Calculate the size of a colour palette.
+.TP
+.BR libquanta_malloc_palette (3)
+Allocate a colour palette.
+
+.SS Colour quantisation functions
+.TP
+.BR libquanta_quantise (3),
+.TQ
+.BR libquanta_vquantise (3)
+Colour-quantise an image using an algorithm selected by the library.
+.TP
+.TP
+.BR libquanta_quantise_wu (3),
+.TQ
+.BR libquanta_vquantise_wu (3)
+Colour-quantise an image using Wu's Colour Quantiser.
+
+.SH SEE ALSO
+None.
diff --git a/libquanta.h b/libquanta.h
new file mode 100644
index 0000000..89c2fa0
--- /dev/null
+++ b/libquanta.h
@@ -0,0 +1,253 @@
+/* See LICENSE file for copyright and license details. */
+#ifndef LIBQUANTA_H
+#define LIBQUANTA_H
+
+#include <stdarg.h>
+#include <stdint.h>
+#include <stddef.h>
+
+
+/**
+ * Colour-quantised image
+ */
+struct libquanta_image {
+	/**
+	 * The image's width
+	 */
+	size_t width;
+
+	/**
+	 * The image's height
+	 */
+	size_t height;
+
+	/**
+	 * `.image[y * .width + x]` for `y` in [0, `.height`), `x` in [0, `.width`),
+	 * will be the index of the colour at pixel (`x`, `y`). This index is the
+	 * colour index in the colour palette stored separatedly in the `.palette`
+	 * field of a `struct libquanta_palette`.
+	 */
+	size_t image[];
+};
+
+
+/**
+ * A colour channel from an image
+ */
+struct libquanta_channel {
+	/**
+	 * The number of bits per input image colour value,
+	 * must be between 1 and 64 (inclusive)
+	 */
+	size_t bits_in;
+
+	/**
+	 * The number of bits per palette colour value,
+	 * must be between 1 and `.bits_in` (inclusive)
+	 */
+	size_t bits_out;
+
+	/**
+	 * The image's width in cells
+	 */
+	size_t image_width;
+
+	/**
+	 * The image's height in cells
+	 */
+	size_t image_height;
+
+	/**
+	 * If image is row-major:
+	 *    The number of bytes per row in the image
+	 * 
+	 * If image is column-major:
+	 *    The number of bytes per cell in the image
+	 */
+	size_t image_row_size;
+
+	/**
+	 * If image is row-major:
+	 *    The number of bytes per cell in the image
+	 * 
+	 * If image is column-major:
+	 *    The number of bytes per column in the image
+	 */
+	size_t image_cell_size;
+
+	/**
+	 * The image
+	 * 
+	 * Must not be `NULL` unless `.image_width` or `.image_height` is 0
+	 * 
+	 * `&((const char *).image)[y * .image_row_size + x * .image_cell_size]`
+	 * shall point to the colour value for cell (`x`, `y`), for the colour
+	 * channel`, encoded in host-native endian; for
+	 * 1 ≤ `.bits_in` ≤ 8, it shall be a `uint8_t`, for
+	 * 9 ≤ `.bits_in` ≤ 16, it shall be a `uint16_t`, for
+	 * 17 ≤ `.bits_in` ≤ 32, it shall be a `uint32_t`, and for
+	 * 33 ≤ `.bits_in` ≤ 64, it shall be a `uint64_t`
+	 */
+	const void *image;
+};
+
+
+/**
+ * Colour palette for an image
+ */
+struct libquanta_palette {
+	/**
+	 * The number of colours stored in the palette, must be at least 1
+	 */
+	size_t size;
+
+	/**
+	 * The colour palette
+	 * 
+	 * If the image has N colour channels, colour channel value i (within [0, N))
+	 * for palette colour j (within [0, `.size`)), is `.palette[j * N + i]`
+	 */
+	uint64_t palette[];
+};
+
+
+/**
+ * Get the allocation size for a `struct libquanta_palette`
+ * 
+ * @param   ncolours   The number of colours the palette may hold
+ * @param   nchannels  The number of colour channels
+ * @param   size_out   Output parameter for the allocation size
+ *                     (may be `NULL`)
+ * @return             0 on success, -1 on failure
+ * 
+ * @throws  EOVERFLOW  The palette is too large to be allocated
+ * @throws  EINVAL     `ncolours` or `nchannels` is 0
+ */
+int libquanta_palette_size(size_t ncolours, size_t nchannels, size_t *size_out);
+
+/**
+ * Allocate and initialise a colour palette
+ * 
+ * @param   ncolours   The number of colours the palette may hold, must be at least 1
+ * @param   nchannels  The number of colour channels, must be at least 1
+ * @return             Colour palette object with `.size` set to `ncolours`,
+ *                     allocated using malloc(3) (caller is responsible for
+ *                     deallocating it using free(3); `NULL` on failure
+ * 
+ * @throws  ENOMEM  The palette is too large to be allocated,
+ *                  or there is not enough memory to allocate it
+ * @throws  EINVAL  `ncolours` or `nchannels` is 0
+ */
+struct libquanta_palette *libquanta_malloc_palette(size_t ncolours, size_t nchannels);
+
+
+/**
+ * Colour-quantise an image using some library-selected colour quantiser
+ * 
+ * @param   palette  Colour palette to fill, `.size` must be set with the
+ *                   maximum number of colours, and will upon successful
+ *                   completion be set to the number of used colours
+ * @param   args     One `const struct libquanta_channel *` per colour
+ *                   channel, providing the image for each colour channel;
+ *                   the list shall be terminated by `NULL`
+ * @return           Colour-quantised image, allocated using malloc(3)
+ *                   (caller is responsible for deallocating it using
+ *                   free(3); `NULL` on failure
+ * 
+ * @throws  EINVAL  `palette` is `NULL`, `palette->size` is 0,
+ *                  or no colour channels are provided (the first argument
+ *                  stored in `args` is `NULL`)
+ * @throws  ENOMEM  Not enough memory available
+ *
+ * The function may write beyond to and beyond `&palette->palette[palette->size]`
+ * for `palette->size` as set after the function returns, but not as set when
+ * the function is called. The function may also have written to `palette->palette`
+ * even if the function fails.
+ */
+struct libquanta_image *libquanta_vquantise(struct libquanta_palette *palette, va_list args);
+
+/**
+ * Colour-quantise an image using some library-selected colour quantiser
+ *
+ * @param   palette  Colour palette to fill, `.size` must be set with the
+ *                   maximum number of colours, and will upon successful
+ *                   completion be set to the number of used colours
+ * @param   ...      One `const struct libquanta_channel *` per colour
+ *                   channel, providing the image for each colour channel;
+ *                   the list shall be terminated by `NULL`
+ * @return           Colour-quantised image, allocated using malloc(3)
+ *                   (caller is responsible for deallocating it using
+ *                   free(3); `NULL` on failure
+ * 
+ * @throws  EINVAL  `palette` is `NULL`, `palette->size` is 0,
+ *                  or no colour channels are provided (second argument is `NULL`)
+ * @throws  EINVAL  One of the colour channel descriptions has an invalid configuration
+ *                  or the colour channel descriptions do not have the same image size
+ *                  configured
+ * @throws  ENOMEM  Not enough memory available
+ *
+ * The function may write beyond to and beyond `&palette->palette[palette->size]`
+ * for `palette->size` as set after the function returns, but not as set when
+ * the function is called. The function may also have written to `palette->palette`
+ * even if the function fails.
+ */
+struct libquanta_image *libquanta_quantise(struct libquanta_palette *palette, ...);
+
+
+/**
+ * Colour-quantise an image using Wu's Colour Quantiser
+ * 
+ * @param   palette  Colour palette to fill, `.size` must be set with the
+ *                   maximum number of colours, and will upon successful
+ *                   completion be set to the number of used colours
+ * @param   args     One `const struct libquanta_channel *` per colour
+ *                   channel, providing the image for each colour channel;
+ *                   the list shall be terminated by `NULL`
+ * @return           Colour-quantised image, allocated using malloc(3)
+ *                   (caller is responsible for deallocating it using
+ *                   free(3); `NULL` on failure
+ * 
+ * @throws  EINVAL  `palette` is `NULL`, `palette->size` is 0,
+ *                  or no colour channels are provided (the first argument
+ *                  stored in `args` is `NULL`)
+ * @throws  EINVAL  One of the colour channel descriptions has an invalid configuration
+ *                  or the colour channel descriptions do not have the same image size
+ *                  configured
+ * @throws  ENOMEM  Not enough memory available
+ *
+ * The function may write beyond to and beyond `&palette->palette[palette->size]`
+ * for `palette->size` as set after the function returns, but not as set when
+ * the function is called. The function may also have written to `palette->palette`
+ * even if the function fails.
+ */
+struct libquanta_image *libquanta_vquantise_wu(struct libquanta_palette *palette, va_list args);
+
+/**
+ * Colour-quantise an image using Wu's Colour Quantiser
+ *
+ * @param   palette  Colour palette to fill, `.size` must be set with the
+ *                   maximum number of colours, and will upon successful
+ *                   completion be set to the number of used colours
+ * @param   ...      One `const struct libquanta_channel *` per colour
+ *                   channel, providing the image for each colour channel;
+ *                   the list shall be terminated by `NULL`
+ * @return           Colour-quantised image, allocated using malloc(3)
+ *                   (caller is responsible for deallocating it using
+ *                   free(3); `NULL` on failure
+ * 
+ * @throws  EINVAL  `palette` is `NULL`, `palette->size` is 0,
+ *                  or no colour channels are provided (second argument is `NULL`)
+ * @throws  EINVAL  One of the colour channel descriptions has an invalid configuration
+ *                  or the colour channel descriptions do not have the same image size
+ *                  configured
+ * @throws  ENOMEM  Not enough memory available
+ *
+ * The function may write beyond to and beyond `&palette->palette[palette->size]`
+ * for `palette->size` as set after the function returns, but not as set when
+ * the function is called. The function may also have written to `palette->palette`
+ * even if the function fails.
+ */
+struct libquanta_image *libquanta_quantise_wu(struct libquanta_palette *palette, ...);
+
+
+#endif
diff --git a/libquanta_bigint_divmod_small__.c b/libquanta_bigint_divmod_small__.c
new file mode 100644
index 0000000..a2c345f
--- /dev/null
+++ b/libquanta_bigint_divmod_small__.c
@@ -0,0 +1,34 @@
+/* See LICENSE file for copyright and license details. */
+#include "common.h"
+
+
+uintmax_t
+libquanta_bigint_divmod_small__(struct bigint *big, uintmax_t small)
+{
+	uintmax_t q = 0, hi, lo;
+	int e = 8 * (int)sizeof(small);
+
+#if 0 /* this would overflow (undefined behaviour) and not fit in the result */
+	q = big->high / small;
+	q <<= 8 * (int)sizeof(small);
+#endif
+	big->high %= small;
+
+	while (big->high && --e) {
+		hi = small >> (8 * (int)sizeof(small) - e);
+		lo = small << e;
+
+		if (hi > big->high)
+			continue;
+		if (hi == big->high && lo > big->low)
+			continue;
+
+		q |= (uintmax_t)1 << e;
+		bigint_sub_small(big, lo);
+		big->high -= hi;
+	}
+
+	q += big->low / small;
+	big->low %= small;
+	return q;
+}
diff --git a/libquanta_channel.3 b/libquanta_channel.3
new file mode 120000
index 0000000..baf5586
--- /dev/null
+++ b/libquanta_channel.3
@@ -0,0 +1 @@
+struct_libquanta_channel.3
\ No newline at end of file
diff --git a/libquanta_image.3 b/libquanta_image.3
new file mode 120000
index 0000000..e5d16c4
--- /dev/null
+++ b/libquanta_image.3
@@ -0,0 +1 @@
+struct_libquanta_image.3
\ No newline at end of file
diff --git a/libquanta_malloc_palette.3 b/libquanta_malloc_palette.3
new file mode 100644
index 0000000..5523b4e
--- /dev/null
+++ b/libquanta_malloc_palette.3
@@ -0,0 +1,73 @@
+.TH LIBQUANTA_MALLOC_PALETTE 3 libquanta
+.SH NAME
+libquanta_malloc_palette \- Allocation a palette
+
+.SH SYNPOSIS
+.nf
+#include <libquanta.h>
+
+\fBstruct libquanta_palette\fP {
+	size_t \fIsize\fP;
+	uint64_t \fIpalette\fP[];
+};
+
+struct libquanta_palette *\fBlibquanta_malloc_palette\fP(size_t \fIncolours\fP, size_t \fInchannels\fP);
+.fi
+.PP
+Link with
+.IR -lquanta .
+
+.SH DESCRIPTION
+The
+.BR libquanta_malloc_palette ()
+function allocates a new palette holding a maximum of
+.I ncolours
+colours, each with
+.I nchannels
+primary colours.
+.PP
+.I .size
+in the structure refered to by the returned pointer will
+be set to
+.IR ncolours .
+.PP
+.I ncolours
+shall be the maximum number of colours the palette may hold.
+Must be at least 1.
+.PP
+.I nchannels
+shall be the number of colour channels in the image.
+Must be at least 1.
+
+.SH RETURN VALUE
+Upon successful completion, the
+.BR libquanta_malloc_palette ()
+function returns a pointer to newly palette.
+The caller is responsable for deallocating the
+memory, that the returned value points to, using the
+.BR free (3)
+function.
+On failure, the function return
+.I NULL
+and set
+.I errno
+to indicate the error.
+
+.SH ERRORS
+The
+.BR libquanta_malloc_palette ()
+function may fail if:
+.TP
+.B EINVAL
+.I ncolours
+or
+.I nchannels
+is 0.
+.TP
+.B ENOMEM
+Storage space available is insufficient.
+
+.SH SEE ALSO
+.BR libquanta (7),
+.BR libquanta_palette_size (3),
+.BR libquanta_quantise (3)
diff --git a/libquanta_malloc_palette.c b/libquanta_malloc_palette.c
new file mode 100644
index 0000000..816fb89
--- /dev/null
+++ b/libquanta_malloc_palette.c
@@ -0,0 +1,21 @@
+/* See LICENSE file for copyright and license details. */
+#include "common.h"
+
+
+struct libquanta_palette *
+libquanta_malloc_palette(size_t ncolours, size_t nchannels)
+{
+	struct libquanta_palette *ret;
+	size_t size;
+
+	if (libquanta_palette_size(ncolours, nchannels, &size)) {
+		if (errno == EOVERFLOW)
+			errno = ENOMEM;
+		return NULL;
+	}
+
+	ret = malloc(size);
+	if (ret)
+		ret->size = ncolours;
+	return ret;
+}
diff --git a/libquanta_palette.3 b/libquanta_palette.3
new file mode 120000
index 0000000..3029ff5
--- /dev/null
+++ b/libquanta_palette.3
@@ -0,0 +1 @@
+struct_libquanta_palette.3
\ No newline at end of file
diff --git a/libquanta_palette_size.3 b/libquanta_palette_size.3
new file mode 100644
index 0000000..146bde3
--- /dev/null
+++ b/libquanta_palette_size.3
@@ -0,0 +1,66 @@
+.TH LIBQUANTA_PALETTE_SIZE 3 libquanta
+.SH NAME
+libquanta_palette_size \- Calculate allocation size of palette
+
+.SH SYNPOSIS
+.nf
+#include <libquanta.h>
+
+int \fBlibquanta_palette_size\fP(size_t \fIncolours\fP, size_t \fInchannels\fP, size_t *\fIsize_out\fP);
+.fi
+.PP
+Link with
+.IR -lquanta .
+
+.SH DESCRIPTION
+The
+.BR libquanta_palette_size ()
+function calculates the allocation size for a
+.I struct libquanta_palette
+whuch may hold
+.I ncolours
+colours, each with
+.I nchannels
+primary colours. The size is stored in
+.I *size_out
+unless
+.I size_out
+is
+.IR NULL .
+.PP
+.I ncolours
+shall be the maximum number of colours the palette may hold.
+Must be at least 1.
+.PP
+.I nchannels
+shall be the number of colour channels in the image.
+Must be at least 1.
+
+.SH RETURN VALUE
+Upon successful completion, the
+.BR libquanta_palette_size ()
+function returns 0.
+On failure, the function return -1 and set
+.I errno
+to indicate the error.
+
+.SH ERRORS
+The
+.BR libquanta_palette_size ()
+function may fail if:
+.TP
+.B EINVAL
+.I ncolours
+or
+.I nchannels
+is 0.
+.TP
+.B EOVERFLOW
+The palette is too large to be allocated
+because the require size exceeds the maximum
+allocation size for the platform.
+
+.SH SEE ALSO
+.BR libquanta (7),
+.BR libquanta_malloc_palette (3),
+.BR libquanta_quantise (3)
diff --git a/libquanta_palette_size.c b/libquanta_palette_size.c
new file mode 100644
index 0000000..e8b95eb
--- /dev/null
+++ b/libquanta_palette_size.c
@@ -0,0 +1,22 @@
+/* See LICENSE file for copyright and license details. */
+#include "common.h"
+
+
+int
+libquanta_palette_size(size_t ncolours, size_t nchannels, size_t *size_out)
+{
+	if (!ncolours || !nchannels) {
+		errno = EINVAL;
+		return -1;
+	}
+
+	if (ncolours > (SIZE_MAX - PALETTE_BASE_SIZE) / PALETTE_VALUE_SIZE / nchannels) {
+		errno = EOVERFLOW;
+		return -1;
+	}
+
+	if (size_out)
+		*size_out = PALETTE_BASE_SIZE + ncolours * nchannels * PALETTE_VALUE_SIZE;
+
+	return 0;
+}
diff --git a/libquanta_quantise.3 b/libquanta_quantise.3
new file mode 100644
index 0000000..3a5db49
--- /dev/null
+++ b/libquanta_quantise.3
@@ -0,0 +1,246 @@
+.TH LIBQUANTA_QUANTISE 3 libquanta
+.SH NAME
+libquanta_quantise \- Colour-quantise an image
+
+.SH SYNPOSIS
+.nf
+#include <libquanta.h>
+
+\fBstruct libquanta_image\fP {
+	size_t \fIwidth\fP;
+	size_t \fIheight\fP;
+	size_t \fIimage\fP[];
+};
+
+\fBstruct libquanta_palette\fP {
+	size_t \fIsize\fP;
+	uint64_t \fIpalette\fP[];
+};
+
+\fBstruct libquanta_channel\fP {
+	size_t \fIbits_in\fP;
+	size_t \fIbits_out\fP;
+	size_t \fIimage_width\fP;
+	size_t \fIimage_height\fP;
+	size_t \fIimage_row_size\fP;
+	size_t \fIimage_cell_size\fP;
+	const void *\fIimage\fP;
+};
+
+struct libquanta_image *\fBlibquanta_quantise\fP(struct libquanta_palette *\fIpalette\fP, ...);
+struct libquanta_image *\fBlibquanta_vquantise\fP(struct libquanta_palette *\fIpalette\fP, va_list \fIargs\fP);
+.fi
+.PP
+Link with
+.I -lquanta
+.IR -lm .
+
+.SH DESCRIPTION
+The
+.BR libquanta_quantise ()
+creates and returnes a colour-quantised image from
+input image and desired palette size.
+.PP
+.I palette->size
+shall the desired palette size: the maximum number
+of colours. This value must be at least 1.
+.I *palette
+may be otherwise uninitialised.
+.PP
+Upon successful completion,
+.I palette->size
+will be set to the actual number of colours stored
+in the palette, which will not exceed the input
+value of
+.IR palette->size.
+Additionally,
+.I palette->palette
+will be filled with the colours in the colour
+palette choosen by the function. For a colour index
+.I i
+(integer within [0,
+.IR palette->size )
+(using the output value of
+.IR palette->size ))
+and channel index
+.I j
+(integer within [0,
+.IR n ),
+where
+.I n
+is the number of colour channels (the number of
+arguments between
+.I palette
+and the
+.I NULL
+at the end of the argument list),
+.I palette->palette[i*n+j]
+is the value of the
+.IR j -th
+primary colour of the
+.IR i -th
+palette colour.
+On failure, arbitrary values may be stored in
+.IR palette->palette .
+.PP
+After
+.I palette
+the caller shall provide one
+.B const struct libquanta_channel *
+per primary colour (at least one) in the image's
+colour model. This list shall be terminated by a
+.I NULL
+as the final argument.
+For each such argument,
+.I .bits_in
+shall be the number of bits used to stored colour
+values for that colour channel in the input image.
+This value must be at least 1, but no greater than 64.
+.I .bits_out
+shall be the number of bits used to stored colour
+values for that colour channel in the output colour
+palette. This value must be at least 1, but no
+greater than
+.IR .bits_in .
+.I .image_width
+shall be the number of pixels the image is wide.
+.I .image_height
+shall be the number of pixels the image is tall.
+.IR .image ,
+.IR .image_row_size ,
+and
+.IR image_cell_size
+shall be configured so that
+.I &((const char *).image)[y * .image_row_size + x * .image_cell_size]
+points to the value of pixel
+.RI ( x ,
+.IR y )
+for the colour channel.
+Although the naming of
+.I .image_row_size
+and
+.I .image_cell_size
+imply row-major layout, column-major layout
+is also supported. If
+.I .bits_in
+is between 1 and 8 (inclusively),
+values in
+.I .image
+shall use the type
+.BR uint8_t .
+If
+.I .bits_in
+is between 9 and 16,
+.B uint16_t
+is used. If
+.I .bits_in
+is between 17 and 32,
+.B uint32_t
+is used. If
+.I .bits_in
+is between 33 and 64,
+.B uint64_t
+is used.
+.PP
+For the returned
+.BR "struct libquanta_image" ,
+.I width
+will be set to the image width
+.RI ( .image_width
+used for the input
+.IR "struct libquanta_channel" s),
+.I height
+will be set to the image height
+.RI ( .image_height
+used for the input
+.IR "struct libquanta_channel" s),
+and
+.I .image[y*.width+x]
+for each
+.I y
+in [0,
+.IR .height )
+and each
+.I x
+in [0,
+.IR .width )
+will be set to the index of the
+palette colour selected for pixel
+.RI ( x ,
+.IR y ),
+will be less than the value assigned to
+.I palette->size
+when the function returns.
+.PP
+Unused elements in
+.I palette->palette
+may remain uninitialised or contain
+arbitrary values when the function
+returns.
+.PP
+The
+.BR libquanta_vquantise ()
+function is a variant of the
+.BR libquanta_quantise ()
+function, that uses
+.I va_list args
+to hold variadic arguments.
+
+.SH RETURN VALUE
+Upon successful completion, the
+.BR libquanta_quantise (3)
+and
+.BR libquanta_vquantise (3)
+functions return a pointer to a newly
+allocated colour-quantised image. The
+caller is responsible for deallocating
+the image by calling the
+.BR free (3)
+with the returned pointer.
+On failure, the functions return
+.I NULL
+and set
+.I errno
+to indicate the error.
+
+.SH ERRORS
+The
+.BR libquanta_quantise ()
+and
+.BR libquanta_vquantise ()
+functions will fail if:
+.TP
+.B EINVAL
+.I palette
+is
+.IR NULL .
+.TP
+.B EINVAL
+.I palette->size
+is 0.
+.TP
+.B EINVAL
+No colour cannel is provided (the first
+variadic argument is
+.IR NULL ).
+.TP
+.B EINVAL
+At least one of the colour channel descriptions
+have an invalid configuration
+.RI ( .bits_in
+is less than 1 or greater than 64, or
+.I .bits_out
+is less than 1 or greater than
+.IR .bits_in ).
+.TP
+.B EINVAL
+The colour channel descriptions do not have the
+same image dimensions configured.
+.TP
+.B ENOMEM
+Storage space available is insufficient.
+
+.SH SEE ALSO
+.BR libquanta (7),
+.BR libquanta_malloc_palette (3),
+.BR libquanta_quantise_wu (3)
diff --git a/libquanta_quantise.c b/libquanta_quantise.c
new file mode 100644
index 0000000..8d7c3e5
--- /dev/null
+++ b/libquanta_quantise.c
@@ -0,0 +1,12 @@
+/* See LICENSE file for copyright and license details. */
+#include "common.h"
+
+
+struct libquanta_image *
+libquanta_quantise(struct libquanta_palette *palette, ...)
+{
+	va_list args;
+	va_start(args, palette);
+	return libquanta_vquantise(palette, args);
+	va_end(args);
+}
diff --git a/libquanta_quantise_wu.3 b/libquanta_quantise_wu.3
new file mode 100644
index 0000000..d7673e8
--- /dev/null
+++ b/libquanta_quantise_wu.3
@@ -0,0 +1,50 @@
+.TH LIBQUANTA_QUANTISE_WU 3 libquanta
+.SH NAME
+libquanta_quantise_wu \- Colour-quantise using Wu's Colour Quantiser
+
+.SH SYNPOSIS
+.nf
+#include <libquanta.h>
+
+struct libquanta_image *\fBlibquanta_quantise_wu\fP(struct libquanta_palette *\fIpalette\fP, ...);
+struct libquanta_image *\fBlibquanta_vquantise_wu\fP(struct libquanta_palette *\fIpalette\fP, va_list \fIargs\fP);
+.fi
+.PP
+Link with
+.I -lquanta
+.IR -lm .
+
+.SH DESCRIPTION
+The
+.BR libquanta_quantise_wu ()
+and
+.BR libquanta_vquantise_wu ()
+functions are variants of the
+.BR libquanta_quantise (3)
+and
+.BR libquanta_vquantise (3)
+functions. These variants always use
+Wu's Colour Quantiser.
+.PP
+See
+.BR libquanta_quantise (3)
+for more information.
+
+.SH RETURN VALUE
+See
+.BR libquanta_quantise (3).
+
+.SH ERRORS
+The
+.BR libquanta_quantise_wu ()
+and
+.BR libquanta_vquantise_wu ()
+functions may fail for any response
+specified for the
+.BR libquanta_quantise (3)
+function.
+
+.SH SEE ALSO
+.BR libquanta (7),
+.BR libquanta_quantise (3),
+.BR libquanta_malloc_palette (3)
diff --git a/libquanta_quantise_wu.c b/libquanta_quantise_wu.c
new file mode 100644
index 0000000..c7fb483
--- /dev/null
+++ b/libquanta_quantise_wu.c
@@ -0,0 +1,12 @@
+/* See LICENSE file for copyright and license details. */
+#include "common.h"
+
+
+struct libquanta_image *
+libquanta_quantise_wu(struct libquanta_palette *palette, ...)
+{
+	va_list args;
+	va_start(args, palette);
+	return libquanta_vquantise_wu(palette, args);
+	va_end(args);
+}
diff --git a/libquanta_vquantise.3 b/libquanta_vquantise.3
new file mode 120000
index 0000000..77c4960
--- /dev/null
+++ b/libquanta_vquantise.3
@@ -0,0 +1 @@
+libquanta_quantise.3
\ No newline at end of file
diff --git a/libquanta_vquantise.c b/libquanta_vquantise.c
new file mode 100644
index 0000000..eb97fac
--- /dev/null
+++ b/libquanta_vquantise.c
@@ -0,0 +1,9 @@
+/* See LICENSE file for copyright and license details. */
+#include "common.h"
+
+
+struct libquanta_image *
+libquanta_vquantise(struct libquanta_palette *palette, va_list args)
+{
+	return libquanta_vquantise_wu(palette, args);
+}
diff --git a/libquanta_vquantise_wu.3 b/libquanta_vquantise_wu.3
new file mode 120000
index 0000000..56b70b8
--- /dev/null
+++ b/libquanta_vquantise_wu.3
@@ -0,0 +1 @@
+libquanta_quantise_wu.3
\ No newline at end of file
diff --git a/libquanta_vquantise_wu.c b/libquanta_vquantise_wu.c
new file mode 100644
index 0000000..6e1ce94
--- /dev/null
+++ b/libquanta_vquantise_wu.c
@@ -0,0 +1,713 @@
+s/* See LICENSE file for copyright and license details. */
+#include "common.h"
+
+/* TODO limitation: 32 values per channel, limits palette size */
+#define WORKING_BITS 5U
+
+
+struct channel_data {
+	const struct libquanta_channel *ch;
+	size_t subindex_multiplier;
+	size_t colour_shift;
+	struct bigint *histogram;
+	struct bigint whole;
+	struct bigint half;
+	struct bigint base;
+	long double max;
+	size_t cut;
+	int has_cut;
+};
+
+struct histogram {
+	size_t occurrences;
+	long double square_sum;
+};
+
+struct cube_boundary {
+	size_t min; /* exclusive */
+	size_t max; /* inclusive */
+};
+
+struct cube {
+	size_t volume;
+	struct cube_boundary bounds[];
+};
+
+
+static void
+compute_histogram(struct libquanta_image *image, struct histogram *histogram, struct channel_data *channels, size_t nchannels)
+{
+	size_t x, width  = channels[0].ch->image_width;
+	size_t y, height = channels[0].ch->image_height;
+	size_t ch, i, index, pixel;
+	long double square_subsum;
+	const char *p;
+	uint64_t v;
+
+	for (y = 0, pixel = 0; y < height; y++) {
+		for (x = 0; x < width; x++, pixel++) {
+			index = 0;
+			square_subsum = 0;
+			for (ch = nchannels; ch--;) {
+				i  = y * channels[ch].ch->image_row_size;
+				i += x * channels[ch].ch->image_cell_size;
+				p = &((const char *)channels[ch].ch->image)[i];
+				if (channels[ch].ch->bits_in <= 8)
+					v = *(const uint8_t *)p;
+				else if (channels[ch].ch->bits_in <= 16)
+					v = *(const uint16_t *)p;
+				else if (channels[ch].ch->bits_in <= 32)
+					v = *(const uint32_t *)p;
+				else if (channels[ch].ch->bits_in <= 64)
+					v = *(const uint64_t *)p;
+				else
+					abort();
+				square_subsum += (long double)v * (long double)v;
+				v >>= channels[ch].colour_shift;
+				v += 1U;
+				index += v * channels[ch].subindex_multiplier;
+				bigint_add_small(&channels[ch].histogram[index], v);
+			}
+			histogram[index].occurrences += 1U;
+			histogram[index].square_sum += square_subsum;
+			image->image[pixel] = index;
+		}
+	}
+}
+
+
+static void
+compute_moments(struct histogram *histogram, struct channel_data *channels, size_t nchannels) /* TODO */
+{
+}
+
+
+static void
+volume_over_square_sums_(const struct cube *cube, const struct histogram *data, const struct channel_data *channels, size_t nchannels,
+                         size_t index, int use_addition, long double *result)
+{
+	size_t ch;
+
+	if (!nchannels) {
+		if (use_addition)
+			*result += data[index].square_sum;
+		else
+			*result -= data[index].square_sum;
+		return;
+	}
+
+	ch = nchannels - 1U;
+
+	volume_over_square_sums_(cube, data, channels, ch,
+	                         index + cube->bounds[ch].max * channels[ch].subindex_multiplier,
+	                         use_addition ^ 1, result);
+
+	volume_over_square_sums_(cube, data, channels, ch,
+	                         index + cube->bounds[ch].min * channels[ch].subindex_multiplier,
+	                         use_addition, result);
+}
+
+static long double
+volume_over_square_sums(const struct cube *cube, const struct histogram *data, const struct channel_data *channels, size_t nchannels)
+{
+	long double result = 0;
+	int use_addition = (int)(~nchannels & 1U);
+	volume_over_square_sums_(cube, data, channels, nchannels, 0, use_addition, &result);
+	return result;
+}
+
+
+static void
+volume_over_occurrences_(const struct cube *cube, const struct histogram *data, const struct channel_data *channels, size_t nchannels,
+                         size_t index, int use_addition, uintmax_t *result)
+{
+	size_t ch;
+
+	if (!nchannels) {
+		uintmax_t v = (uintmax_t)data[index].occurrences;
+		if (use_addition)
+			*result += v;
+		else
+			*result -= v;
+		return;
+	}
+
+	ch = nchannels - 1U;
+
+	volume_over_occurrences_(cube, data, channels, ch,
+	                         index + cube->bounds[ch].max * channels[ch].subindex_multiplier,
+	                         use_addition ^ 1, result);
+
+	volume_over_occurrences_(cube, data, channels, ch,
+	                         index + cube->bounds[ch].min * channels[ch].subindex_multiplier,
+	                         use_addition, result);
+}
+
+static uintmax_t
+volume_over_occurrences(const struct cube *cube, const struct histogram *data, const struct channel_data *channels, size_t nchannels)
+{
+	uintmax_t result = 0;
+	int use_addition = (int)(~nchannels & 1U);
+	volume_over_occurrences_(cube, data, channels, nchannels, 0, use_addition, &result);
+	return result;
+}
+
+
+static void
+volume_over_bigints_(const struct cube *cube, const struct bigint *data, const struct channel_data *channels, size_t nchannels,
+                     size_t index, int use_addition, struct bigint *result)
+{
+	size_t ch;
+
+	if (!nchannels) {
+		if (use_addition)
+			bigint_add_big(result, &data[index]);
+		else
+			bigint_sub_big(result, &data[index]);
+		return;
+	}
+
+	ch = nchannels - 1U;
+
+	volume_over_bigints_(cube, data, channels, ch,
+	                     index + cube->bounds[ch].max * channels[ch].subindex_multiplier,
+	                     use_addition ^ 1, result);
+
+	volume_over_bigints_(cube, data, channels, ch,
+	                     index + cube->bounds[ch].min * channels[ch].subindex_multiplier,
+	                     use_addition, result);
+}
+
+static void
+volume_over_bigints(const struct cube *cube, const struct bigint *data, const struct channel_data *channels, size_t nchannels,
+                    struct bigint *result)
+{
+	int use_addition = (int)(~nchannels & 1U);
+	result->high = 0;
+	result->low = 0;
+	volume_over_bigints_(cube, data, channels, nchannels, 0, use_addition, result);
+}
+
+
+static void
+bottom_over_occurrences_(const struct cube *cube, size_t channel, const struct histogram *data, const struct channel_data *channels,
+                         size_t nchannels, size_t index, int use_addition, uintmax_t *result)
+{
+	size_t ch;
+
+	if (!nchannels) {
+		uintmax_t v = (uintmax_t)data[index].occurrences;
+		if (use_addition)
+			*result += v;
+		else
+			*result -= v;
+		return;
+	}
+
+	ch = nchannels - 1U;
+
+	bottom_over_occurrences_(cube, channel, data, channels, ch,
+	                         index + cube->bounds[ch].min * channels[ch].subindex_multiplier,
+	                         use_addition, result);
+
+	if (ch == channel)
+		return;
+
+	bottom_over_occurrences_(cube, channel, data, channels, ch,
+	                         index + cube->bounds[ch].max * channels[ch].subindex_multiplier,
+	                         use_addition ^ 1, result);
+}
+
+static uintmax_t
+bottom_over_occurrences(const struct cube *cube, size_t channel, const struct histogram *data, const struct channel_data *channels,
+                        size_t nchannels)
+{
+	uintmax_t result = 0;
+	int use_addition = (int)(~nchannels & 1U);
+	bottom_over_occurrences_(cube, channel, data, channels, nchannels, 0, use_addition, &result);
+	return result;
+}
+
+
+static void
+bottom_over_bigints_(const struct cube *cube, size_t channel, const struct bigint *data, const struct channel_data *channels,
+                     size_t nchannels, size_t index, int use_addition, struct bigint *result)
+{
+	size_t ch;
+
+	if (!nchannels) {
+		if (use_addition)
+			bigint_add_big(result, &data[index]);
+		else
+			bigint_sub_big(result, &data[index]);
+		return;
+	}
+
+	ch = nchannels - 1U;
+
+	bottom_over_bigints_(cube, channel, data, channels, ch,
+	                     index + cube->bounds[ch].min * channels[ch].subindex_multiplier,
+	                     use_addition, result);
+
+	if (ch == channel)
+		return;
+
+	bottom_over_bigints_(cube, channel, data, channels, ch,
+	                     index + cube->bounds[ch].max * channels[ch].subindex_multiplier,
+	                     use_addition ^ 1, result);
+}
+
+static void
+bottom_over_bigints(const struct cube *cube, size_t channel, const struct bigint *data, const struct channel_data *channels,
+                    size_t nchannels, struct bigint *result)
+{
+	int use_addition = (int)(~nchannels & 1U);
+	result->high = 0;
+	result->low = 0;
+	bottom_over_bigints_(cube, channel, data, channels, nchannels, 0, use_addition, result);
+}
+
+
+static void
+top_over_occurrences_(const struct cube *cube, size_t channel, size_t position, const struct histogram *data,
+                      const struct channel_data *channels, size_t nchannels, size_t index, int use_addition, uintmax_t *result)
+{
+	size_t ch;
+
+	if (!nchannels) {
+		uintmax_t v = (uintmax_t)data[index].occurrences;
+		if (use_addition)
+			*result += v;
+		else
+			*result -= v;
+		return;
+	}
+
+	ch = nchannels - 1U;
+	
+	if (ch == channel) {
+		top_over_occurrences_(cube, channel, position, data, channels, ch,
+		                      index + position * channels[ch].subindex_multiplier,
+		                      use_addition ^ 1, result);
+
+	} else {
+		top_over_occurrences_(cube, channel, position, data, channels, ch,
+		                      index + cube->bounds[ch].min * channels[ch].subindex_multiplier,
+		                      use_addition, result);
+
+		top_over_occurrences_(cube, channel, position, data, channels, ch,
+		                      index + cube->bounds[ch].max * channels[ch].subindex_multiplier,
+		                      use_addition ^ 1, result);
+	}
+}
+
+static uintmax_t
+top_over_occurrences(const struct cube *cube, size_t channel, size_t position, const struct histogram *data,
+                     const struct channel_data *channels, size_t nchannels)
+{
+	uintmax_t result = 0;
+	int use_addition = (int)(~nchannels & 1U);
+	top_over_occurrences_(cube, channel, position, data, channels, nchannels, 0, use_addition, &result);
+	return result;
+}
+
+
+static void
+top_over_bigints_(const struct cube *cube, size_t channel, size_t position, const struct bigint *data,
+                  const struct channel_data *channels, size_t nchannels, size_t index, int use_addition, struct bigint *result)
+{
+	size_t ch;
+
+	if (!nchannels) {
+		if (use_addition)
+			bigint_add_big(result, &data[index]);
+		else
+			bigint_sub_big(result, &data[index]);
+		return;
+	}
+
+	ch = nchannels - 1U;
+	
+	if (ch == channel) {
+		top_over_bigints_(cube, channel, position, data, channels, ch,
+		                  index + position * channels[ch].subindex_multiplier,
+		                  use_addition ^ 1, result);
+
+	} else {
+		top_over_bigints_(cube, channel, position, data, channels, ch,
+		                  index + cube->bounds[ch].min * channels[ch].subindex_multiplier,
+		                  use_addition, result);
+
+		top_over_bigints_(cube, channel, position, data, channels, ch,
+		                  index + cube->bounds[ch].max * channels[ch].subindex_multiplier,
+		                  use_addition ^ 1, result);
+	}
+}
+
+static void
+top_over_bigints(const struct cube *cube, size_t channel, size_t position, const struct bigint *data,
+                 const struct channel_data *channels, size_t nchannels, struct bigint *result)
+{
+	int use_addition = (int)(~nchannels & 1U);
+	result->high = 0;
+	result->low = 0;
+	top_over_bigints_(cube, channel, position, data, channels, nchannels, 0, use_addition, result);
+}
+
+
+static void
+maximise(const struct cube *cube, size_t channel, const struct histogram *histogram, struct channel_data *channels,
+         size_t nchannels, uintmax_t whole_weight)
+{
+	uintmax_t base_weight;
+	uintmax_t half_weight;
+	size_t i, pos;
+	size_t first, last;
+	long double temp, temp2, v;
+
+	channels[channel].max = 0;
+	channels[channel].cut = 0;
+	channels[channel].has_cut = 0;
+
+	for (i = 0; i < nchannels; i++)
+		bottom_over_bigints(cube, channel, channels[i].histogram, channels, nchannels, &channels[i].base);
+	base_weight = bottom_over_occurrences(cube, channel, histogram, channels, nchannels);
+
+	first = cube->bounds[channel].min + 1U;
+	last = cube->bounds[channel].max;
+	for (pos = first; pos < last; pos++) { /* TODO should `last` really be exclusive? */
+		half_weight = top_over_occurrences(cube, channel, pos, histogram, channels, nchannels);
+		half_weight += base_weight;
+		if (!half_weight)
+			continue;
+
+		temp = 0;
+		for (i = 0; i < nchannels; i++) {
+			top_over_bigints(cube, channel, pos, channels[i].histogram, channels, nchannels, &channels[i].half);
+			bigint_add_big(&channels[i].half, &channels[i].base);
+
+			v = (long double)channels[i].half.high;
+			v = ldexpl(v, 8 * sizeof(channels[i].half.low));
+			v += (long double)channels[i].half.low;
+			temp = fmal(v, v, temp);
+		}
+		temp /= (long double)half_weight;
+
+		half_weight = whole_weight - half_weight;
+		if (!half_weight)
+			continue;
+
+		temp2 = 0;
+		for (i = 0; i < nchannels; i++) {
+			bigint_rsub_big(&channels[i].half, &channels[i].whole); /* half = whole - half */
+
+			v = (long double)channels[i].half.high;
+			v = ldexpl(v, 8 * sizeof(channels[i].half.low));
+			v += (long double)channels[i].half.low;
+			temp2 = fmal(v, v, temp2);
+		}
+		temp += temp2 / (long double)half_weight;
+
+		if (temp > channels[channel].max) {
+			channels[channel].max = temp;
+			channels[channel].cut = pos;
+			channels[channel].has_cut = 1;
+		}
+	}
+}
+
+
+static int
+cut(struct cube *cube1, struct cube *cube2, const struct histogram *histogram, struct channel_data *channels, size_t nchannels)
+{
+	size_t i, ch;
+	uintmax_t whole_weight;
+	long double max_max;
+
+	for (i = 0; i < nchannels; i++)
+		volume_over_bigints(cube1, channels[i].histogram, channels, nchannels, &channels[i].whole);
+	whole_weight = volume_over_occurrences(cube1, histogram, channels, nchannels);
+
+	for (i = 0; i < nchannels; i++)
+		maximise(cube1, i, histogram, channels, nchannels, whole_weight);
+
+	max_max = channels[0].max;
+	ch = 0;
+	for (i = 1; i < nchannels; i++) {
+		if (channels[i].max > max_max) {
+			max_max = channels[i].max;
+			ch = i;
+		}
+	}
+
+	if (!channels[ch].has_cut)
+		return 0;
+
+	for (i = 0; i < nchannels; i++) {
+		cube2->bounds[i].max = cube1->bounds[i].max;
+		cube2->bounds[i].min = cube1->bounds[i].min;
+	}
+
+	cube2->bounds[ch].min = cube1->bounds[ch].max = channels[ch].cut;
+
+	cube1->volume = 1U;
+	cube2->volume = 1U;
+	for (i = 0; i < nchannels; i++) {
+		cube1->volume *= cube1->bounds[i].max - cube1->bounds[i].min;
+		cube2->volume *= cube2->bounds[i].max - cube2->bounds[i].min;
+	}
+
+	return 1;
+}
+
+
+static long double
+variance(const struct cube *cube, const struct histogram *histogram, const struct channel_data *channels, size_t nchannels)
+{
+	uintmax_t weight;
+	long double w, x, y = 0;
+	struct bigint v;
+	size_t i;
+
+	x = volume_over_square_sums(cube, histogram, channels, nchannels);
+
+	for (i = 0; i < nchannels; i++) {
+		volume_over_bigints(cube, channels[i].histogram, channels, nchannels, &v);
+		w = (long double)v.high;
+		w = ldexpl(w, 8 * sizeof(v.low));
+		w += (long double)v.low;
+		y = fmal(w, w, y);
+	}
+
+	weight = volume_over_occurrences(cube, histogram, channels, nchannels);
+	y /= (long double)weight;
+
+	return x - y;
+}
+
+
+static void
+mark_(const struct cube *cube, size_t label, size_t *tag, const struct channel_data *channels, size_t nchannels, size_t index)
+{
+	size_t i, ch;
+
+	if (!nchannels) {
+		tag[index] = label;
+		return;
+	}
+
+	ch = nchannels - 1U;
+	for (i = cube->bounds[ch].min; i++ < cube->bounds[ch].max;)
+		mark_(cube, label, tag, channels, ch, index + i * channels[ch].subindex_multiplier);
+}
+
+static void
+mark(const struct cube *cube, size_t label, size_t *tag, const struct channel_data *channels, size_t nchannels)
+{
+	mark_(cube, label, tag, channels, nchannels, 0);
+}
+
+
+static int
+partition(struct cube **cubes, const struct histogram *histogram, struct channel_data *channels, size_t nchannels,
+          size_t ncolours, size_t *ncubes_out)
+{
+	size_t i, j, next;
+	long double temp;
+	long double *v;
+
+	if (ncolours > SIZE_MAX / sizeof(*v)) {
+		errno = ENOMEM;
+		return -1;
+	}
+	v = malloc(ncolours * sizeof(*v));
+	if (!v)
+		return -1;
+
+	for (i = 0; i < nchannels; i++) {
+		cubes[0]->bounds[i].min = 0;
+		cubes[0]->bounds[i].max = (1U << WORKING_BITS);
+	}
+
+	next = 0;
+	for (i = 1; i < ncolours; i++) {
+		if (cut(cubes[next], cubes[i], histogram, channels, nchannels)) {
+			v[next] = cubes[next]->volume > 1U ? variance(cubes[next], histogram, channels, nchannels) : 0U;
+			v[i] = cubes[i]->volume > 1U ? variance(cubes[i], histogram, channels, nchannels) : 0U;
+		} else {
+			v[next] = 0;
+			i--;
+		}
+		next = 0;
+		temp = v[0];
+		for (j = 0; j++ < i;)
+			if (v[i] > temp)
+				temp = v[next = j];
+		if (temp <= 0) {
+			i++;
+			break;
+		}
+	}
+	*ncubes_out = i;
+
+	free(v);
+	return 0;
+}
+
+
+struct libquanta_image *
+libquanta_vquantise_wu(struct libquanta_palette *palette, va_list args)
+{
+	struct channel_data *channels = NULL;
+	const struct libquanta_channel *channel;
+	va_list args2;
+	size_t ncolours, nchannels;
+	size_t width = 0, height = 0;
+	size_t i, j, n, subindex_multiplier;
+	struct histogram *histogram = NULL;
+	struct libquanta_image *ret = NULL;
+	size_t size, image_size, *tag = NULL;
+	uintmax_t weight, mean;
+	struct bigint sum;
+	struct cube **cubes = NULL;
+	size_t cube_size;
+	size_t ncolours_used;
+
+	/* inspect argument `palette`, and get number of output colours */
+	if (!palette || !palette->size)
+		goto einval;
+	ncolours = palette->size;
+
+	/* inspect argument `args`, and get number of input colour channels */
+	va_copy(args2, args);
+	nchannels = 0;
+	for (; (channel = va_arg(args2, const struct libquanta_channel *)); nchannels++) {
+		if (!nchannels) {
+			width = channel->image_width;
+			height = channel->image_height;
+		} else if (channel->image_width != width || channel->image_height != height) {
+			goto einval;
+		}
+		if (!channel->bits_in || channel->bits_in > PALETTE_VALUE_MAX_BITS)
+			goto einval;
+		if (!channel->bits_out || channel->bits_out > channel->bits_in)
+			goto einval;
+		if (!channel->image && (channel->image_width | channel->image_height))
+			goto einval;
+	}
+	if (!nchannels)
+		goto einval;
+	va_end(args2);
+
+	/* allocate working memory for each channel, and save channel argument and precompute data */
+	channels = calloc(nchannels, sizeof(*channels));
+	if (!channels)
+		goto fail;
+	for (i = 0; (channel = va_arg(args, const struct libquanta_channel *)); i++) {
+		channels[i].ch = channel;
+		channels[i].histogram = NULL;
+	}
+	subindex_multiplier = 1U;
+	for (i = nchannels; i--; subindex_multiplier *= ((1U << WORKING_BITS) + 1U)) {
+		channels[i].subindex_multiplier = subindex_multiplier;
+		channels[i].colour_shift = channels[i].ch->bits_in > WORKING_BITS ? channels[i].ch->bits_in - WORKING_BITS : 0U;
+	}
+	for (i = 0; i < nchannels; i++) {
+		channels[i].histogram = calloc(subindex_multiplier, sizeof(*channels[i].histogram));
+		if (!channels[i].histogram)
+			goto fail;
+	}
+
+	/* allocate working memory that is not per channel */
+	histogram = calloc(subindex_multiplier, sizeof(*histogram));
+	if (!histogram)
+		goto fail;
+	if (subindex_multiplier > SIZE_MAX / sizeof(*tag))
+		goto enomem;
+	tag = malloc(subindex_multiplier * sizeof(*tag));
+	if (!tag)
+		goto fail;
+	if (ncolours > SIZE_MAX / sizeof(*cubes))
+		goto enomem;
+	cubes = malloc(ncolours * sizeof(*cubes));
+	if (!cubes)
+		goto fail;
+	for (i = 0; i < ncolours; i++)
+		cubes[i] = NULL;
+	if (nchannels > (SIZE_MAX - offsetof(struct cube, bounds)) / sizeof(struct cube_boundary))
+		goto enomem;
+	cube_size = offsetof(struct cube, bounds);
+	cube_size += nchannels * sizeof(struct cube_boundary);
+	for (i = 0; i < ncolours; i++) {
+		cubes[i] = malloc(cube_size);
+		if (!cubes[i])
+			goto fail;
+	}
+
+	/* allocate and configure colour-quantised image */
+	size = sizeof(*ret->image);
+	image_size = 0;
+	if (width && height) {
+		if (width > (SIZE_MAX - offsetof(struct libquanta_image, image)) / size / height)
+			goto enomem;
+		size *= image_size = width * height;
+	}
+	size += offsetof(struct libquanta_image, image);
+	ret = malloc(offsetof(struct libquanta_image, image) + size);
+	if (!ret)
+		goto fail;
+	ret->width = width;
+	ret->height = height;
+
+	/* analyse image */
+	compute_histogram(ret, histogram, channels, nchannels);
+	compute_moments(histogram, channels, nchannels);
+	if (partition(cubes, histogram, channels, nchannels, ncolours, &ncolours_used))
+		goto fail;
+
+	/* reduce image */
+	n = 0;
+	for (i = 0; i < ncolours_used; i++) {
+		mark(cubes[i], i, tag, channels, nchannels);
+		weight = volume_over_occurrences(cubes[i], histogram, channels, nchannels);
+		if (!weight)
+			continue;
+		for (j = 0; j < nchannels; j++) {
+			volume_over_bigints(cubes[i], channels[j].histogram, channels, nchannels, &sum);
+			mean = bigint_div_small(&sum, weight);
+			if (channels[j].ch->bits_out < channels[j].ch->bits_in)
+				mean >>= channels[j].ch->bits_in - channels[j].ch->bits_out;
+			palette->palette[n * nchannels + j] = (uint64_t)mean;
+		}
+		n++;
+	}
+	palette->size = n;
+	for (i = 0; i < image_size; i++)
+		ret->image[i] = tag[ret->image[i]];
+
+	goto out;
+
+enomem:
+	errno = ENOMEM;
+	goto fail;
+
+einval:
+	errno = EINVAL;
+fail:
+	free(ret);
+	ret = NULL;
+out:
+	free(tag);
+	free(histogram);
+	if (channels) {
+		for (i = 0; i < nchannels; i++)
+			free(channels[i].histogram);
+		free(channels);
+	}
+	if (cubes) {
+		for (i = 0; i < ncolours; i++)
+			free(cubes[i]);
+		free(cubes);
+	}
+	return ret;
+}
diff --git a/mk/linux.mk b/mk/linux.mk
new file mode 100644
index 0000000..ad58f69
--- /dev/null
+++ b/mk/linux.mk
@@ -0,0 +1,6 @@
+LIBEXT      = so
+LIBFLAGS    = -shared -Wl,-soname,lib$(LIB_NAME).$(LIBEXT).$(LIB_MAJOR)
+LIBMAJOREXT = $(LIBEXT).$(LIB_MAJOR)
+LIBMINOREXT = $(LIBEXT).$(LIB_VERSION)
+
+FIX_INSTALL_NAME = :
diff --git a/mk/macos.mk b/mk/macos.mk
new file mode 100644
index 0000000..256ee5c
--- /dev/null
+++ b/mk/macos.mk
@@ -0,0 +1,6 @@
+LIBEXT      = dylib
+LIBFLAGS    = -dynamiclib -Wl,-compatibility_version,$(LIB_MAJOR) -Wl,-current_version,$(LIB_VERSION)
+LIBMAJOREXT = $(LIB_MAJOR).$(LIBEXT)
+LIBMINOREXT = $(LIB_VERSION).$(LIBEXT)
+
+FIX_INSTALL_NAME = install_name_tool -id "$(PREFIX)/lib/libquanta.$(LIBMAJOREXT)"
diff --git a/mk/windows.mk b/mk/windows.mk
new file mode 100644
index 0000000..ed5ec8d
--- /dev/null
+++ b/mk/windows.mk
@@ -0,0 +1,6 @@
+LIBEXT      = dll
+LIBFLAGS    = -shared
+LIBMAJOREXT = $(LIB_MAJOR).$(LIBEXT)
+LIBMINOREXT = $(LIB_VERSION).$(LIBEXT)
+
+FIX_INSTALL_NAME = :
diff --git a/struct_libquanta_channel.3 b/struct_libquanta_channel.3
new file mode 120000
index 0000000..77c4960
--- /dev/null
+++ b/struct_libquanta_channel.3
@@ -0,0 +1 @@
+libquanta_quantise.3
\ No newline at end of file
diff --git a/struct_libquanta_image.3 b/struct_libquanta_image.3
new file mode 120000
index 0000000..77c4960
--- /dev/null
+++ b/struct_libquanta_image.3
@@ -0,0 +1 @@
+libquanta_quantise.3
\ No newline at end of file
diff --git a/struct_libquanta_palette.3 b/struct_libquanta_palette.3
new file mode 120000
index 0000000..77c4960
--- /dev/null
+++ b/struct_libquanta_palette.3
@@ -0,0 +1 @@
+libquanta_quantise.3
\ No newline at end of file