Convert 32 bit PNGs to high quality 8 bit PNGs with C# / by Matt Wrock

I’ve been working on a C# HTTPModule called RequestReduce that sprites CSS background images on the fly as well as merges and Minifes all CSS resources in the document head. For those who may not know, image spriting is a technique that takes individual images in separate image files and merges them into a single file. By tracking the offsets of where each image begins, you can use CSS syntax to limit a css background image to just one of the images in the file like so:

background: url(/RequestReduceContent/94b703a23cf3a2644577ab68fa465844-e2d1ced7efd8cd883516e2da4e51a9d5.png) -241 0 no-repeat

By using this technique effectively, page load times will be faster since there will be fewer HTTP connections opened.

One of the items on my backlog has been image optimization. I had envisioned a feature that would either shell out to one of the popular PNG compressors or allow users a means of plugging in a compressor of their choice. About a week or two before releasing this project at work where I work on the Microsoft EPX Galleries team where we produce among many other galleries, the Visual Studio Gallery and the MSDN Code Samples Gallery, I was in a meeting discussing image optimization and someone was demonstrating the benefits of optimizing some of the images that get published in the MSDN dev centers. The impact on size was enormous and on the level of producing images several times smaller than their unoptimized originals.

So I figured how difficult could it be to add this feature to RequestReduce? My plan was to shell out to optipng.exe - a very popular and effective C Library that losslessly compresses PNG images. The answer was not difficult at all and the code looked something like this:

private void InvokeExecutable(string arguments, string executable){    using(var process = new Process())    {        process.StartInfo = new ProcessStartInfo()        {            UseShellExecute = false,            RedirectStandardOutput = true,            CreateNoWindow = true,            FileName = executable,            Arguments = arguments        };        process.Start();        process.StandardOutput.ReadToEnd();        process.WaitForExit(10000);        if(!process.HasExited)        {            process.Kill();            throw new OptimizationException                (string.Format("Unable to optimize image using executable {0} with arguments {1}",                 executable, arguments));        }    }}

However, I was seeing image size reductions that were much less dramatic than the ones demoed in the meeting I attended. The reason was that I work with Web Devs who are typically good about optimizing images and the presenter in the meeting was referring to “site managers” who might know little to nothing about image optimization and sometimes upload huge images unaware of the consequences.

So I was a bit disappointed but decided to move forward with the feature. I was after all seeing some improvement but more on the order of 5%.

A couple days later I plugged RequestReduce into this blog. While it was very effective in reducing the number of HTTP requests, it appeared as though the total count of downloaded bytes was greater that before implementing RequestReduce. How could this be? With the minification of the CSS, it would certainly be less right? Apparently not. The main culprit was a JPG image that was originally abut  5k and was consuming well over this amount in my sprite.I figured this was a flaw I had to fix and I knew very little about the details of image optimization.

My first inclination was that I was creating the image incorrectly. I assumed there was some encoder setting I was failing to set correctly that was forcing the generation of 32 bit PNGs. I spent hours googling only to find there were countless others with the same inclination and no solution. Previously I had reviewed the code written for the Asp.Net Sprite Framework. This framework does something similar to RequestReduce but takes a very different approach that requires some intentional organization of the local image files before it can work. It had code, and eventually I had code that looked like this when finally generating the final PNG bitmap:

public byte[] GetBytes(string mimeType){    using (var spriteEncoderParameters = new EncoderParameters(1))    {        spriteEncoderParameters.Param[0] = new EncoderParameter(Encoder.Quality, 90);        using (var stream = new MemoryStream())        {            SpriteImage.Save(stream, ImageCodecInfo.GetImageEncoders()            .First(x => x.MimeType == mimeType), spriteEncoderParameters);            return stream.GetBuffer();        }    }}

Let me just set the record straight here because there are oodles of forum posts and stackoverflow answers concerning this and many give wrong or misleading answers. The Encoder.Quality parameter does absolutely nothing in the default PNG encoder as well as most of the other encoder parameters. I had thought that this was where the bulk of my problem lied. It wasn’t.

Eventually it became clear that this was not going to be solved by some hidden property in the GDI+ api which is the API exposed in .net BCL for image manipulation. I stumbled upon a concept called image quantization. The process of reducing the number of colors in an image. This is a key concept in getting a 32 bit PNG to convert to a 8 bit PNG because an 8 bit PNG can have no more than 256 colors. It has what is called an Indexed Color palette. A part of the PNG file structure holds pointers to 256 colors and then each pixel in the image gets its color from one of those pointers. Thus each pixel only consumes one bytem its 0-255 value pointing to its color on the palette. On the other hand, a 32 bit PNG is 4 bytes per pixel and each pixel can represent a different ARGB color value. A 32 bit PNG also has a lot of wasted space since in may have multiple pixels that use the same color but they each hold their own copy.

So with this information it makes sense why my sprites had to be 32 bit PNGs. As I added images to the sprites, even if the individual images being added were 8 bit PNGs, I would inevitably accrue more than 256 colors, once that happens, an 8 bit PNG can no longer be used without some data loss. That said, data loss does not necessarily need to translate to “perceptible” quality loss. This is where color quantizers come into play.

Typically image optimization for the web takes one to two paths of optimization: lossless and lossy. Lossless guaeantees no quality loss and implements various compression strategies to shrink the size of a PNG. Popular tools like Optipng or PngCrush utilize this method, The savings can sometimes be significant and other times not so much. However, the savings is always “free” of quality loss and is therefore highly advisable.

Lossy compression does certainly run the risk of unacceptable quality loss but can often render images that are dramatically smaller than their 32 bit cousins and their quality loss is practically imperceptible. There is inevitably loss, but if the loss can either not be perceived at all or is so slight that the savings far outweigh the color loss, taking on the loss may be advisable. Often the deciding factor is the number of colors in the original image, The closer that number is to 256, the less quality loss there will be and therefore the less likely that any of this loss will be perceptible. RequestReduce’s default color limit to quantize is 5000. However, this rule is not concrete. There are images where noticable color loss may occur beyond 3000 colors and other times (especially dealing with photos, where 10000 images is fine to quantize.

So I only managed to find one open source .net based library that provided quantization and I was not comfortable with its license. There are a couple of C based EXEs that perform fast and effective quantization. The ones I used were PNGQuant and Pngnq. The difference is the algorithm they employ for the quantization. I found that these produced decent quality images but one of my images (a rating star) was always off color. and a couple sprites had images with gradients that did look rather awful after the quantization. I eventially tried tweaking my code to limit the number of colors stored in a single sprite. This meant producing up to 3 sprite images for a page which seemed to me to be compromising the benefit of the sprite.

So I looked further into other quantization algorithms. I found an old sample on MSDN that I also found in an older version of Paint.net source code. My images looked terrible using this sample. Then I stumbled upon this article on Code Project. It was a C# library that contained several different quantization algorithms. They all looked quite subpar except for one: an algorithm developed by Xiaolin Wu. My images looked perfect. I could not notice any loss. Furthermore, their sizes were even smaller than the images generated by the C exe libraries. There was just one catch: the algorithm assumed RGB values pixels with no alpha opacity value. Now this Code Project sample included an “alpha blending” technique that made the transparency “appear” to be visible but you had to know what the background color is to blend the transparency with. RequestReduce works with sites that it has no idea what the background color is.

Thus began a three week journey to tweak Xiaolin Wu’s algorithm to include the fourth alpha layer. Let me just preface this with the fact that I have no math or statistics background to speak of. Sadly, someone with such a background probably would have spent a lot less time on this.But it was literally almost all I could think of for three weeks. It was the type of problem that I constantly felt I was about to solve within an hour but each breakthrough simply exposed a new problem that became more complex than I anticipated as I got closer.Before really jumping in, I had thought that adding the fourth Alpha dimension to RGB would be rather trivial. Just look for lines like:

halfDistance = halfRed * halfRed + halfGreen * halfGreen + halfBlue * halfBlue;

and change it to:

halfDistance = halfRed * halfRed + halfGreen * halfGreen + halfBlue * halfBlue + halfAlpha * halfAlpha;

I can do that and a lot of it was just that. Except for a few key parts. I’m not going to agonize the details of this discovery process. While perhaps theraputic to myself, it would have negative value for any reader. Suffice it to say that in the end after lots of scratch pads and notpad.exe windows full of random numbers and sequences, I eventually had a quantized image with 256 colors that looked really good, better and smaller than the ones produced by pngquant or pngnq. However they were not perfect.

As one can assume, adding an extra byte for alpha will translate to a source image .with potentially much more colors than a fully opaque image. This is because, technically, the same color with a different alpha value, is a different value and competes for one of the 256 slots in the 8 bit PNG indexed palette. In fact, there is the potential for 256x more colors. So I was not surprised to find that these quantized images were of slightly lower quality that the ones produced by the RGB quantizer.

Fortunately, in reality, the number of values used in the alpha channel are usually far less than those used in the RGB channels.That being the case, I took some liberties with the alpha values that increased the image quality while sacrificing the number of unique alpha values used. I used two strategies here:

1. I created an AlphaThreshold value. This is a constant that I set but should be configurable if this were production code. Any alpha value equal to or below the threshold gets translated to 0. I set this threshold to 30, finding that values with an alpha of 30 or less are for the most part invisible.

2. I normalized the source alpha values to their closest multiple of 70. Why 70? Because that is the number transmitted to me from the great master Astra who is the giver of justice in the curious singing forest. Gotta love Astra, A heart of gold that one..This essentially limits the number of values that the alpha channel can occupy but continues to allow for decent gradients. I’m sure there are images where this would not work well but it looks good on all the images I have processed so far.

So lets dive into the code, You can find a complete working VS Solution on the MSDN Samples Gallery. The process of quantization can be split into the following parts:

  1. Read the unoptimized source image and create a histogram of its colors. Which colors are used and how often.
  2. Break down this data into several cumulative arrays that will make it much more efficient to slice and dice.
  3. Create 256 clusters of colors from which we will pick the resulting palette. This step is what sets Wu’s algorithm apart from others. It is a process of perpetually splitting the color space into pairs of boxes to find clusters of colors with minimal variance.
  4. Iterate over the source image again to find find the actual color closest to the mean of each cluster (these will be our palette colors) and map each pixel to its appropriate cluster.
  5. Generate the new quantized image from the data in step 4.

Again, what lends to the accuracy of this algorithm is the clustering of colors of minimal variance. It is looking for the core colors in the image and not simply the colors most often used. With lower quality algorithms, it is far more likely for parts of an image to appear discolored as a more predominant color in the image. The color may be somewhat similar but quite noticeably different. From my tests, the Wu algorithm seems to eliminate these inaccuracies. For a deeper analysis from the author himself see  Graphics Gems vol. II, pp. 126-133.

Building the Histogram


private static ColorData BuildHistogram(Bitmap sourceImage){    var data = sourceImage.LockBits(Rectangle.FromLTRB(0, 0, sourceImage.Width, sourceImage.Height),                                    ImageLockMode.ReadOnly, sourceImage.PixelFormat);    var colorData = new ColorData(MaxSideIndex);

    try    {        var byteLength = data.Stride < 0 ? -data.Stride : data.Stride;        var byteCount = Math.Max(1, BitDepth >> 3);        var offset = 0;        var buffer = new Byte[byteLength * sourceImage.Height];        var value = new Byte[byteCount];

        Marshal.Copy(data.Scan0, buffer, 0, buffer.Length);        for (var y = 0; y < sourceImage.Height; y++)        {            var index = 0;            for (var x = 0; x < sourceImage.Width; x++)            {                var indexOffset = index >> 3;

                for (var valueIndex = 0; valueIndex < byteCount; valueIndex++)                    value[valueIndex] = buffer[offset + valueIndex + indexOffset];

                var indexAlpha = (byte)((value[Alpha] >> 3) + 1);                var indexRed = (byte)((value[Red] >> 3) + 1);                var indexGreen = (byte)((value[Green] >> 3) + 1);                var indexBlue = (byte)((value[Blue] >> 3) + 1);

                if (value[Alpha] > AlphaThreshold)                {                    if (value[Alpha] < 255)                    {                        var alpha = value[Alpha] + (value[Alpha] % 70);                        value[Alpha] = (byte)(alpha > 255 ? 255 : alpha);                        indexAlpha = (byte)((value[Alpha] >> 3) + 1);                    }

                    colorData.Weights[indexAlpha, indexRed, indexGreen, indexBlue]++;                    colorData.MomentsRed[indexAlpha, indexRed, indexGreen, indexBlue] 

                     += value[Red];                    colorData.MomentsGreen[indexAlpha, indexRed, indexGreen, indexBlue] 

                     += value[Green];                    colorData.MomentsBlue[indexAlpha, indexRed, indexGreen, indexBlue] 

                     += value[Blue];                    colorData.MomentsAlpha[indexAlpha, indexRed, indexGreen, indexBlue] 

                     += value[Alpha];                    colorData.Moments[indexAlpha, indexRed, indexGreen, indexBlue] 

                     += (value[Alpha]*value[Alpha]) +                                                              

                      (value[Red]*value[Red]) +                      (value[Green]*value[Green]) +                      (value[Blue]*value[Blue]);                }

                colorData.QuantizedPixels.Add(BitConverter.ToInt32(new[] { indexAlpha, indexRed, 

                 indexGreen, indexBlue }, 0));                colorData.Pixels.Add(new Pixel(value[Alpha], value[Red], value[Green], 

                 value[Blue]));                index += BitDepth;            }            offset += byteLength;        }    }    finally    {        sourceImage.UnlockBits(data);    }    return colorData;}

This is rather straight forward with just two things to note. First, note that there are a few ways to iterate over the pixels of a bit map. Three that I can think of:

  1. The simplest but by far the most inefficient is to simply create a for loop for x and y and call GetPixel(x,y) which returns the color of that pixel.
  2. Use lockbits and unlock bits (as seen above) but with the difference of using pointers to iterate over the locked memory space of the image. This is much faster than technique number one but requires you to compile the code as unsafe.
  3. Use lockbits and unlock but use Marshal.Copy to feed the memory into a buffer.

Second, note that we drop the three rightmost bit of each color dimension.This reduces the granularity of the color maps produced in the next step, making it much faster, This allows us to create our clusters using color values from 1 to 32 instead of 256.

Creating the Color Arrays

These are the data structures that all of the upcoming analysis is performed upon. If this step is off, everything else will be off. This is one of the steps that really led me down a rat hole. For the longest time I thought this was correct and was looking for bugs elsewhere.


private static ColorData CalculateMoments(ColorData data){    for (var alphaIndex = 1; alphaIndex <= MaxSideIndex; ++alphaIndex)    {        var xarea = new long[SideSize, SideSize, SideSize];        var xareaAlpha = new long[SideSize, SideSize, SideSize];        var xareaRed = new long[SideSize, SideSize, SideSize];        var xareaGreen = new long[SideSize, SideSize, SideSize];        var xareaBlue = new long[SideSize, SideSize, SideSize];        var xarea2 = new float[SideSize, SideSize, SideSize];        for (var redIndex = 1; redIndex <= MaxSideIndex; ++redIndex)        {            var area = new long[SideSize];            var areaAlpha = new long[SideSize];            var areaRed = new long[SideSize];            var areaGreen = new long[SideSize];            var areaBlue = new long[SideSize];            var area2 = new float[SideSize];            for (var greenIndex = 1; greenIndex <= MaxSideIndex; ++greenIndex)            {                long line = 0;                long lineAlpha = 0;                long lineRed = 0;                long lineGreen = 0;                long lineBlue = 0;                var line2 = 0.0f;                for (var blueIndex = 1; blueIndex <= MaxSideIndex; ++blueIndex)                {                    line += data.Weights[alphaIndex, redIndex, greenIndex, blueIndex];                    lineAlpha += data.MomentsAlpha[alphaIndex, redIndex, greenIndex, blueIndex];                    lineRed += data.MomentsRed[alphaIndex, redIndex, greenIndex, blueIndex];                    lineGreen += data.MomentsGreen[alphaIndex, redIndex, greenIndex, blueIndex];                    lineBlue += data.MomentsBlue[alphaIndex, redIndex, greenIndex, blueIndex];                    line2 += data.Moments[alphaIndex, redIndex, greenIndex, blueIndex];

                    area[blueIndex] += line;                    areaAlpha[blueIndex] += lineAlpha;                    areaRed[blueIndex] += lineRed;                    areaGreen[blueIndex] += lineGreen;                    areaBlue[blueIndex] += lineBlue;                    area2[blueIndex] += line2;

                    xarea[redIndex, greenIndex, blueIndex] = xarea[redIndex - 1, greenIndex, 

                     blueIndex] + area[blueIndex];                    xareaAlpha[redIndex, greenIndex, blueIndex] = xareaAlpha[redIndex - 1, 

                      greenIndex, blueIndex] + areaAlpha[blueIndex];                    xareaRed[redIndex, greenIndex, blueIndex] = xareaRed[redIndex - 1, 

                      greenIndex, blueIndex] + areaRed[blueIndex];                    xareaGreen[redIndex, greenIndex, blueIndex] = xareaGreen[redIndex - 1, 

                      greenIndex, blueIndex] + areaGreen[blueIndex];                    xareaBlue[redIndex, greenIndex, blueIndex] = xareaBlue[redIndex - 1, 

                      greenIndex, blueIndex] + areaBlue[blueIndex];                    xarea2[redIndex, greenIndex, blueIndex] = xarea2[redIndex - 1, 

                      greenIndex, blueIndex] + area2[blueIndex];

                    data.Weights[alphaIndex, redIndex, greenIndex, blueIndex] = 

                      data.Weights[alphaIndex - 1, redIndex, greenIndex, blueIndex] + 

                      xarea[redIndex, greenIndex, blueIndex];                    data.MomentsAlpha[alphaIndex, redIndex, greenIndex, blueIndex] = 

                      data.MomentsAlpha[alphaIndex - 1, redIndex, greenIndex, blueIndex] + 

                      xareaAlpha[redIndex, greenIndex, blueIndex];                    data.MomentsRed[alphaIndex, redIndex, greenIndex, blueIndex] = 

                      data.MomentsRed[alphaIndex - 1, redIndex, greenIndex, blueIndex] + 

                      xareaRed[redIndex, greenIndex, blueIndex];                    data.MomentsGreen[alphaIndex, redIndex, greenIndex, blueIndex] = 

                      data.MomentsGreen[alphaIndex - 1, redIndex, greenIndex, blueIndex] + 

                      xareaGreen[redIndex, greenIndex, blueIndex];                    data.MomentsBlue[alphaIndex, redIndex, greenIndex, blueIndex] = 

                      data.MomentsBlue[alphaIndex - 1, redIndex, greenIndex, blueIndex] + 

                      xareaBlue[redIndex, greenIndex, blueIndex];                    data.Moments[alphaIndex, redIndex, greenIndex, blueIndex] = 

                      data.Moments[alphaIndex - 1, redIndex, greenIndex, blueIndex] + 

                      xarea2[redIndex, greenIndex, blueIndex];                }            }        }    }    return data;}

In the end you should have a set of multi dimensional arrays that cumulatively increase both vertically down each “line” (or most finely grained dimension – blue here) and across to the right from A to R to G to B. The very last value in the array should be the total sum of the original array.

 
Arranging the data n this way allows you to make very efficient measurements on regions of the color space from just a few calculations instead of having to iterate over every point.
 

Clustering the Data

There is a fair amount of code involved here so I will not include all of it. Here is the entry point into this step:


private static IList<Box> SplitData(ref int colorCount, ColorData data){    --colorCount;    var next = 0;    var volumeVariance = new float[MaxColor];    var cubes = new Box[MaxColor];    cubes[0].AlphaMaximum = MaxSideIndex;    cubes[0].RedMaximum = MaxSideIndex;    cubes[0].GreenMaximum = MaxSideIndex;    cubes[0].BlueMaximum = MaxSideIndex;    for (var cubeIndex = 1; cubeIndex < colorCount; ++cubeIndex)    {        if (Cut(data, ref cubes[next], ref cubes[cubeIndex]))        {            volumeVariance[next] = 

              cubes[next].Size > 1 ? CalculateVariance(data, cubes[next]) : 0.0f;            volumeVariance[cubeIndex] = 

              cubes[cubeIndex].Size > 1 ? CalculateVariance(data, cubes[cubeIndex]) : 0.0f;        }        else        {            volumeVariance[next] = 0.0f;            cubeIndex--;        }

        next = 0;        var temp = volumeVariance[0];

        for (var index = 1; index <= cubeIndex; ++index)        {            if (volumeVariance[index] <= temp) continue;            temp = volumeVariance[index];            next = index;        }

        if (temp > 0.0) continue;        colorCount = cubeIndex + 1;        break;    }    return cubes.Take(colorCount).ToList();}

To dive deeper into what is going on, follow the Cut method. Note that this is looking to build 256 clusters of the smallest possible variance. One interesting piece of code a bit deeper down that I would like to point out is the calculation to obtain the volume of each cube. This demonstrates how the cumulative array allows you to make this calculation without iterating each element in the array:

private static long Volume(Box cube, long[,,,] moment){    return (moment[cube.AlphaMaximum, cube.RedMaximum, cube.GreenMaximum, cube.BlueMaximum] -            moment[cube.AlphaMaximum, cube.RedMaximum, cube.GreenMinimum, cube.BlueMaximum] -            moment[cube.AlphaMaximum, cube.RedMinimum, cube.GreenMaximum, cube.BlueMaximum] +            moment[cube.AlphaMaximum, cube.RedMinimum, cube.GreenMinimum, cube.BlueMaximum] -            moment[cube.AlphaMinimum, cube.RedMaximum, cube.GreenMaximum, cube.BlueMaximum] +            moment[cube.AlphaMinimum, cube.RedMaximum, cube.GreenMinimum, cube.BlueMaximum] +            moment[cube.AlphaMinimum, cube.RedMinimum, cube.GreenMaximum, cube.BlueMaximum] -            moment[cube.AlphaMinimum, cube.RedMinimum, cube.GreenMinimum, cube.BlueMaximum]) -

            (moment[cube.AlphaMaximum, cube.RedMaximum, cube.GreenMaximum, cube.BlueMinimum] -            moment[cube.AlphaMinimum, cube.RedMaximum, cube.GreenMaximum, cube.BlueMinimum] -            moment[cube.AlphaMaximum, cube.RedMaximum, cube.GreenMinimum, cube.BlueMinimum] +            moment[cube.AlphaMinimum, cube.RedMaximum, cube.GreenMinimum, cube.BlueMinimum] -            moment[cube.AlphaMaximum, cube.RedMinimum, cube.GreenMaximum, cube.BlueMinimum] +            moment[cube.AlphaMinimum, cube.RedMinimum, cube.GreenMaximum, cube.BlueMinimum] +            moment[cube.AlphaMaximum, cube.RedMinimum, cube.GreenMinimum, cube.BlueMinimum] -            moment[cube.AlphaMinimum, cube.RedMinimum, cube.GreenMinimum, cube.BlueMinimum]);}

Just 15 simple arithmetic operations instead of 32768. This is another formula that I agonized over for days. The trick is figuring out where to put the pluses and minuses. The original RGB method is:

private static Int64 Volume(WuColorCube cube, Int64[, ,] moment){    return moment[cube.RedMaximum, cube.GreenMaximum, cube.BlueMaximum] -           moment[cube.RedMaximum, cube.GreenMaximum, cube.BlueMinimum] -           moment[cube.RedMaximum, cube.GreenMinimum, cube.BlueMaximum] +           moment[cube.RedMaximum, cube.GreenMinimum, cube.BlueMinimum] -           moment[cube.RedMinimum, cube.GreenMaximum, cube.BlueMaximum] +           moment[cube.RedMinimum, cube.GreenMaximum, cube.BlueMinimum] +           moment[cube.RedMinimum, cube.GreenMinimum, cube.BlueMaximum] -           moment[cube.RedMinimum, cube.GreenMinimum, cube.BlueMinimum];}

The similarity is obvious. The key is building the arrays correctly. With incorrect data, this formula will never come out right. I also did not expect the addition of alpha to require parentheses.

 

Generating the Final Palette

This builds the final data structure we need to create out quantized image:


private static QuantizedPalette GetQuantizedPalette(int colorCount, ColorData data, 

                                                    IEnumerable<Box> cubes){    var imageSize = data.Pixels.Count;    var lookups = BuildLookups(cubes, data);

    for(var index = 0; index < imageSize; ++index)    {        var indexParts = BitConverter.GetBytes(data.QuantizedPixels[index]);        data.QuantizedPixels[index] = lookups.Tags[indexParts[Alpha],         indexParts[Red], indexParts[Green], indexParts[Blue]];    }

    var alphas = new int[colorCount + 1];    var reds = new int[colorCount + 1];    var greens = new int[colorCount + 1];    var blues = new int[colorCount + 1];    var sums = new int[colorCount + 1];    var palette = new QuantizedPalette(imageSize);    var index2 = -1;

    foreach (var pixel in data.Pixels)    {        palette.PixelIndex[++index2] = -1;        if(pixel.Alpha <= AlphaThreshold)            continue;

        var match = data.QuantizedPixels[index2];        var bestMatch = match;        var bestDistance = 100000000;        var index = -1;

        foreach (var lookup in lookups.Lookups)        {            ++index;            var deltaAlpha = pixel.Alpha - lookup.Alpha;            var deltaRed = pixel.Red - lookup.Red;            var deltaGreen = pixel.Green - lookup.Green;            var deltaBlue = pixel.Blue - lookup.Blue;

            var distance = deltaAlpha * deltaAlpha + deltaRed * deltaRed             + deltaGreen * deltaGreen + deltaBlue * deltaBlue;

            if (distance >= bestDistance) continue;

            bestDistance = distance;            bestMatch = index;        }

        alphas[bestMatch] += pixel.Alpha;        reds[bestMatch] += pixel.Red;        greens[bestMatch] += pixel.Green;        blues[bestMatch] += pixel.Blue;        sums[bestMatch]++;

        palette.PixelIndex[index2] = bestMatch;    }

    for (var paletteIndex = 0; paletteIndex < colorCount; paletteIndex++)    {        if (sums[paletteIndex] > 0)        {            alphas[paletteIndex] /= sums[paletteIndex];            reds[paletteIndex] /= sums[paletteIndex];            greens[paletteIndex] /= sums[paletteIndex];            blues[paletteIndex] /= sums[paletteIndex];        }

        var color = Color.FromArgb(alphas[paletteIndex], reds[paletteIndex],             greens[paletteIndex], blues[paletteIndex]);        palette.Colors.Add(color);    }    palette.Colors.Add(Color.FromArgb(0, 0, 0, 0));

    return palette;}

private static LookupData BuildLookups(IEnumerable<Box> cubes, ColorData data){    var lookups = new LookupData(SideSize);

    foreach (var cube in cubes)    {        for (var alphaIndex = (byte)(cube.AlphaMinimum + 1); alphaIndex <= cube.AlphaMaximum; ++alphaIndex)        {            for (var redIndex = (byte)(cube.RedMinimum + 1); redIndex <= cube.RedMaximum; ++redIndex)            {                for (var greenIndex = (byte)(cube.GreenMinimum + 1); greenIndex <= cube.GreenMaximum; ++greenIndex)                {                    for (var blueIndex = (byte)(cube.BlueMinimum + 1); blueIndex <= cube.BlueMaximum; ++blueIndex)                        lookups.Tags[alphaIndex, redIndex, greenIndex, blueIndex] = lookups.Lookups.Count;                }            }        }

        var weight = Volume(cube, data.Weights);

        if (weight <= 0) continue;

        var lookup = new Lookup                         {                             Alpha = (int) (Volume(cube, data.MomentsAlpha)/weight),                             Red = (int) (Volume(cube, data.MomentsRed)/weight),                             Green = (int) (Volume(cube, data.MomentsGreen)/weight),                             Blue = (int) (Volume(cube, data.MomentsBlue)/weight)                         };        lookups.Lookups.Add(lookup);    }    return lookups;}

Building the Image

private static Bitmap ProcessImagePixels(Image sourceImage, QuantizedPalette palette){    var result = new Bitmap(sourceImage.Width, sourceImage.Height, PixelFormat.Format8bppIndexed);    var newPalette = result.Palette;    for (var index = 0; index < palette.Colors.Count; index++)        newPalette.Entries[index] = palette.Colors[index];    result.Palette = newPalette;

    BitmapData targetData = null;    try    {        targetData = result.LockBits(Rectangle.FromLTRB(0, 0, result.Width, result.Height),             ImageLockMode.WriteOnly, result.PixelFormat);        const byte targetBitDepth = 8;        var targetByteLength = targetData.Stride < 0 ? -targetData.Stride : targetData.Stride;        var targetByteCount = Math.Max(1, targetBitDepth >> 3);        var targetSize = targetByteLength * result.Height;        var targetOffset = 0;        var targetBuffer = new byte[targetSize];        var targetValue = new byte[targetByteCount];        var pixelIndex = 0;

        for (var y = 0; y < result.Height; y++)        {            var targetIndex = 0;            for (var x = 0; x < result.Width; x++)            {                var targetIndexOffset = targetIndex >> 3;                targetValue[0] = (byte)(palette.PixelIndex[pixelIndex] == -1 ? palette.Colors.Count - 1 :                     palette.PixelIndex[pixelIndex]);                pixelIndex++;

                for (var valueIndex = 0; valueIndex < targetByteCount; valueIndex++)                    targetBuffer[targetOffset + valueIndex + targetIndexOffset] = targetValue[valueIndex];

                targetIndex += targetBitDepth;            }

            targetOffset += targetByteLength;        }

        Marshal.Copy(targetBuffer, 0, targetData.Scan0, targetSize);    }    finally    {        if(targetData != null)            result.UnlockBits(targetData);    }

    return result;}

Here we simply use our generated palette and fill in the pixels.

If you would like to download a full copy of this code along with a .EXE wrapper to test the quality of the quantizations, you can download it from the MSDN Code Samples Gallery.