image processing - Is there any faster implementation of DominantColors?

The function DominantColors is a simple function but runs incredibly slow when the number n of dominant colors to find is large - it is currently the biggest bottleneck in my code:

tst = ExampleData[{"TestImage", "Lena"}];

gr = ListDensityPlot[SeedRandom[1]; RandomReal[{}, {4000, 3}], 
  InterpolationOrder -> 0, Frame -> False, Mesh -> All, 
  ImageSize -> 475]
polys = Cases[Normal@gr, _Polygon, \[Infinity]];
vpolys = Polygon[#[[1]], VertexTextureCoordinates -> #[[1]]] & /@ 
   polys;
pieces = Table[Graphics[{Texture[tst], p}], {p, vpolys}];
colors = Table[
   RandomChoice[DominantColors[p, 5][[2 ;;]]], {p, pieces}];
cpolys = Map[{#[[2]], Polygon[#[[1, 1]]]} &, Thread@{polys, colors}];

gcp = Graphics[cpolys]

Has anyone any ideas to reimplement it to be faster?

Motivation & Clarification

I use mathematica to do explorative art. I'm using this subroutine DominantColors as a tool in my various software art projects for example: I don't need someone to implement a specific image effect, rather I need a faster dominant-n-colors-finder. I'm thinking of using cudalink...

Answer

The function domCol below is about 100 times faster than DominantColors.

Basic plan: Create enough color bins throughout the color space occupied by the image; count the number of pixels in each bin; return the sorted colors. The function works in the LAB color space so we can use EuclideanDistance for the distance between colors. The centers of the bins are jiggled by a random offset each call, so that the return value is not deterministic. The option "RandomSeed" can be used to give a reproducible result.

Examples:

domCol[ExampleData[{"TestImage", "Mandrill"}], 20] // AbsoluteTiming

domCol[ExampleData[{"TestImage", "Lena"}], 20] // AbsoluteTiming

It seems to find similar colors to the ones found by DominantColors, although in a different order here and there.

Code:

Suppose DominantColors spends some time fiddling with constructing an optimal set of color bins, which seems a reasonable hypothesis. This approach does not, but uses a simple heuristic to get a binning that's likely to do a pretty good job. That seems the probable reason for the difference in performance. The function domCol uses a face-center cubic close-packing of spheres, with slightly larger spheres that overlap a little (radius 0.55 dx where dx is the diameter determined by the heuristic). Some of the space is not covered. (Technically, they're not really bins.) The distance can be set explicitly by the option MinColorDistance, which in this case will actually determine the exact color distance between adjacent bins, but the option exists already (for DominantColors). One could make the 0.55 fudge factor an option, too, for more control. With this approach Nearest is a convenient and efficient way to do the bin counts.

ClearAll[domCol];
Options[domCol] = {"RandomSeed" -> Automatic, MinColorDistance -> Automatic};
domCol[img_, n_, OptionsPattern[]] := 
 Module[{idata, nf, bounds, basis, points, colorbins, intensity},
  idata = Flatten[ImageData@ColorConvert[img, "LAB"], 1];
  nf = Nearest@idata;
  bounds = MinMax /@ Transpose@idata;
  With[{dx = 
     OptionValue[MinColorDistance] /.
      Automatic -> (Times @@ Flatten[Differences /@ bounds]/n)^(1/3)/4},

   basis = LatticeData["FaceCenteredCubic", "Basis"];
   If[IntegerQ@OptionValue["RandomSeed"] || 
      StringQ@OptionValue["RandomSeed"],
    SeedRandom[OptionValue["RandomSeed"]]
    ];
   points = Tuples[Range[# - dx + RandomReal[{-dx, dx}/2], #2 + dx, dx] & @@@
      (Sort /@ (Inverse@Transpose[basis].bounds))].basis;
   colorbins = Select[points, Length@nf[#, {1, dx}] >= 1 &];
   intensity = Length@nf[#, {All, 0.55 dx}] & /@ colorbins
   ];

  LABColor /@ colorbins[[Ordering[intensity, -Min[n, Length@colorbins]]]] // Reverse
  ]

The use of the "FaceCenteredCubic" lattice is based in part on this answer by s0rce: Using LatticeData to fill a space with spheres in a face-centered cubic (fcc) lattice packing arrangement