If the number of clusters is known (like here)
You may use Lloyd's clustering [1]
The idea is as follows:
it optimizes a set of cluster centers $p_i$:
Initialize the p_i's with an initial guess, or randomly
For each iteration:
Compute the cluster associated with each p_i,
(the cluster is the set of points nearer to p_i than to the other p_j's)
Move each p_i to the weighted centroid of its cluster
For an image, the iteration can be implented as follows, computing the mass m_i and the centroid g_i of each cluster:
For each i
m_i = 0
g_i = (0,0)
For each pixel (x,y) of the image
let i denote the index of the center p_i nearest to (x,y)
m_i = m_i + pixel_intensity(x,y)
g_i = g_i + pixel_intensity(x,y) * (x,y)
For each i
p_i = (1/m_i)*g_i
Since the number of clusters is small, you can find the nearest p_i using a simple loop. If you have a higher number of sites, you may either use a kd-tree, or compute the Voronoi diagram of the sites and iterate on the pixels of each Voronoi cell.
I used this algorithm to cluster the colors of a rubics cube acquired by a lego color sensor, and it works reasonably well while being very easy to implement [3]
If the number of clusters is unknown
then the problem is much more difficult.
You may use "mean shift clustering" [2], that will apply a filter-like operation to the image, and make the "modes" appear. It acts like the inverse of a smoothing filter.
[1] https://en.wikipedia.org/wiki/Lloyd%27s_algorithm
[2] https://en.wikipedia.org/wiki/Mean_shift
[3] http://alice.loria.fr/WIKI/index.php/Graphite/Lego