I think your algorithm of generating (galaxy) stars is basically not correct, apart from their HR color distribution. There are several better attempts already in the net. In the first place, one has to take the basic inhomogeneity in star density into account that corresponds to the galaxy arms. Here is a close-up of the central quadrant of your M31 galaxy image.
Your star distribution is locally almost homogeneous, while the galaxy "arms" are still consisting of my nebula sprites blended in.
In contrast, here is e.g. an interesting site introducing a generator for (spiral) galaxies
in terms of "stars only":http://dmytry.com/galaxy/index.html
(Dmytry Lavrov's, a mathematically oriented software dev)
The galaxy generator may be downloaded freely. Unfortunately, it's not open source. A simple example thereof illustrates pretty clearly what I mean:
It's straightforward to derive a 2D probability distribution with an arm structure, such that the stars distributed according to it follows the right pattern. It can be easily coded with Von Neuman's well-known acceptance-rejection method
that I have also heavily used in my galaxy (and globular) code for Celestia. A nice summary of that method is here, for example,http://pdg.lbl.gov/2011/reviews/rpp2011 ... niques.pdf
Your star distribution might be useful, however, for modelling elliptical
galaxies as well as nearby spheroidal dwarfs that all are lacking any arms. Also never forget to consider distributing your stars according to a Perlin noise function...
The best approach I can see is a blended
one, depending on the actual screen resolution
or -- equivalently-- the observer distance from the galaxy. As long as the observer is far away, a sprite-based rendering is fast and perfectly sufficient. Upon a closer approach, when the resolution improves a lot, a blended transition to individual stars should be done. A big advantage is that at close distance from the galaxy, one only needs to render a tiny portion
thereof in terms of stars, since the rest is outside the field of view! With clever culling, this fact might make a real-time rendering possible...
Such a two-component rendering has to be set up in the core code, since it is clearly not feasible via add-ons (at least in Celestia).