#AmReading Wikipedia about Expectation Maximization Algorithm for Gaussian Mixture https://en.wikipedia.org/wiki/Expectation%E2%80%93maximization_algorithm#Gaussian_mixture . I think the 4 update equations for T, \tau, \mu, \Sigma in the example section should be sufficient for me to implement it. Generalizing from 2 to arbitrary count should be simple enough too.

I want to implement this as part of a colour quantization algorithm, the other part of the problem I'll try using the Jump Method for determining optimal number of colours, the pseudo code on Wikipedia is fairly clear: https://en.wikipedia.org/wiki/Determining_the_number_of_clusters_in_a_data_set#An_information%E2%80%93theoretic_approach

Worked (almost) first try! (after fixing compilation errors, one stupid typo in a variable name, and having to set ulimit -s unlimited in my shell because I have Big arrays on the stack).

Show thread

Some strange colour choices though, the dark green on the waves seems too dark, and merged with the black for boundaries. Not sure why. Maybe just bad luck in the random number generation found a poor local optimum.

Hmm it's very sensitive, these are all with the same input, only thing changing is the pseudo-random number generator seed... various numbers of colours are chosen as optimal, ranging from 2 to 11.

Probably I should minimize over some randoms in an innerer loop instead of outside the whole process.

Show thread

Follow

Implemented median cut colour quantization as described at https://en.wikipedia.org/wiki/Median_cut

The only non-straightforward part was needing to split at the boundary between median and next value, rather than (potentially) in the middle of a joint median region, to avoid weird artifacts where different parts of the image had different colours.

Hardcoded to 8 colours for now, not sure how to do optimal palette size choice with this method. Maybe https://en.wikipedia.org/wiki/Akaike_information_criterion or similar could work. But I'd need to redo the code structure, currently it's simply recursive (binary splitting, depth 3) and I'd need it to be iterative with a set of blocks to be able to make the decisions about when to stop splitting.