# Median cut algorithm (J)

**Other implementations**: C++ |**J***If you're looking to find the median of a list of elements, see the nth_element standard library function or Category:Selection_algorithm.*

The median cut algorithm is a popular algorithm for color quantization, but can be applied to virtually any point clustering problem. In outline, it works as follows:

- Let B be a set of boxes containing points, initially containing only a single box containing all points.
- While the number of boxes in B is less than desired number of clusters...
- Find the largest side length of any side of any box.
- Cut that box into two boxes along its largest side in such a way that half the contained points fall into each new box.
- Shrink the two new boxes so that they are just large enough to contain their points.

## [edit] Theory

We follow the above algorithm nearly exactly — although instead of *shrinking* bounding boxes, we just recalculate them on each pass.

*Exercise: Our code is specialized for 3-dimensional color data; which definition can be easily generalized for n-dimensional clustering?*

## [edit] Practice

Because each step of the algorithm produces exactly one more box,
and we start with one, we just need to **step**
one fewer times than the number (**x**) of boxes desired.
In J, this is easily accomplished
via the iteration operator (**^:**).
We prepare the initial box from an unboxed array (**<**),
and, having taken the **mean** of each box,
**>** removes the boxes afterwards
to produce an unboxed array of the final palette entries.

<<mediancut.ijs>>=mediancut =: dyad : '> mean each step^:(x-1) <y' mean =: +/%#

At each step,
we must find the longest side among all our boxes (**fmax**)
and use that information to split the largest box (**smax**).

<<mediancut.ijs>>=step =: monad : '(fmax y) smax y'

To find the longest side,
we first compute the bounding box (**bbox**)
by subtracting the minimum value (**<./**) from the maximum value (**>./**)
along each dimension and for each box.
We then collapse these lengths to a 1D list
using the boxed equivalent of ravel (**;**),
sort them (**/:**)
and take the largest (**{:**).

<<mediancut.ijs>>=fmax =: monad : '{:/: ; bbox each y' bbox =: monad : '(>./y)-(<./y)'

By rotating (**|.**) the box list (**y**)
according to the box with the maximum dimension (**<.x%3**),
we reduce the problem of splitting the largest box (**smax**)
to that of splitting the first box (**sfst**) along its largest dimension (**3|x**).

In turn, the problem of splitting the first box
can be reduced to splitting a single box (**sbox**),
operating on the head (**{.y**) of the box list,
and appending the tail (**}.y**) of the box list to the result.

Finally, we split a single box by **rank**ing its entries
according to their values along the desired dimension (**x**),
**chop**ping the entries into
those below and those above the median,
and using the arrax indexing operator (**{**)
to select the appropriate values from the original list (**y**).

<<mediancut.ijs>>=sbox =: dyad : '({&y) each chop x rank y' sfst =: dyad : '(x sbox >{.y),(}.y)' smax =: dyad : '(3|x) sfst (<.x%3) |. y'

**rank**ing the entries is easy:
we generate the necessary permutation (**/:**)
from the values along the desired dimension (**x {"1 y**).

**chop**ping is also straightforward:
we divide (**%**) the number of entries (**#**) in half,
and use the resulting index to produce two boxes:
one with entries before (**{.**)
the other with entries after (**{.**)
the chopping point.
These boxes are joined (**;**)
to produce a single return value.

<<mediancut.ijs>>=rank =: dyad : '/: x {"1 y' chop =: monad : '(<.2%~#y) ({.;}.) y'

## [edit] Wrapping up

We still need to describe file I/O, but assuming one has the raw rgb values in the array **pixeldata**, the following command generates a 16-color palette. (**<.** is the *floor* function, converting floating-point values to integers for the palette)

<. 16 mediancut pixeldata

The output will be:

127 126 112 135 139 139 167 119 112 172 184 101 188 212 112 180 153 158 211 207 190 243 246 220 15 15 13 45 25 32 44 59 26 88 63 46 98 108 54 76 64 89 126 79 102 125 133 84

I reduced color depth in Gimp using this custom palette. With Gimp's dither, the resulting image is a little nicer than the Photoshop/C++ color reduction:

Download code |