The approach I would use, and indeed have used, would be to
- run principal component analysis on your image data to generate "the" main axis
- Project your pixels onto that axis (i.e. compute the dot product of each colour value against the axis)
- Using the dot prods, sort the the colour values from smallest to largest and
- Find the "partitioning" point, i.e. that splits smaller values in one set, larger in the other, that minimises, say, the mean squared error. i.e. each partition is represented by the average of its members and you compute the total error.
This last step, i.e. finding the best partition, can be done in linear time by starting with the "least" item in partition "A", all the rest in partition "B", and progressively shifting the next smaller over from "B" to "A, one at a time. You just need to keep separate sums of the values and sums of squares of the values for each partition and then you can compute the total MSE relatively cheaply.
Note that the function should follow a quadratic curve so it (should) begin decreasing until it reaches the minimum. Once it starts increasing again you can stop the search.
A cheap way to find the principal axis is to first generate the 3x3 covariance matrix and then use the "repeated matrix vector multiplication" trick to find the principal eigenvector (i.e. the principal axis)