Wikipedia about Expectation Maximization Algorithm for Gaussian Mixture en.wikipedia.org/wiki/Expectat . 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: en.wikipedia.org/wiki/Determin

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).

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.

Oh. I think the model is inappropriate for this use case - the mixture of Gaussians means some clusters can be very wide and flat which means the center is not a good approximation of the colour for most of the region... no wonder my results were a bit all over the place...

Implemented median cut colour quantization as described at en.wikipedia.org/wiki/Median_c

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 en.wikipedia.org/wiki/Akaike_i 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.

Welcome to post.lurk.org, an instance for discussions around cultural freedom, experimental, new media art, net and computational culture, and things like that.