One of the effects I needed for a previous game was creating nebula for outerspace. But as a pixel art game, the nebula needed to look… pixel-artsy?
Anyway, my solution was to generate nebula via smooth gradient noise, and then apply the following method to it. The results made the nebula look like they were made specifically for pixel art.
Here’s our example image, and the result with 4 colours and dithering.
In this article I’ll take you through how it’s done. If you’re not familiar with the Median Cut algorithm, or with dithering, you should read about them here, and here.
Define a Pixel and Comparisons
You can probably get away with using whatever existing definition of a pixel you have, along with making some changes to my code; but for illustration purposes, this is a pixel for our algorithm:
struct P
{
std::uint8_t b;
std::uint8_t g;
std::uint8_t r;
std::uint8_t a;
};Code language: C++ (cpp)This works fairly well for OpenGL being in BGRA order. Each channel is 8 bits wide.
Next we need some comparison functions that’ll be used later on. They’re all straightforward.
bool operator==(const P& lhs, const P& rhs)
{
return lhs.r == rhs.r && lhs.g == rhs.g && lhs.b == rhs.b && lhs.a == rhs.a;
}
bool red_comp(const P& a, const P& b)
{
return a.r < b.r;
}
bool green_comp(const P& a, const P& b)
{
return a.g < b.g;
}
bool blue_comp(const P& a, const P& b)
{
return a.b < b.b;
}
bool alpha_comp(const P& a, const P& b)
{
return a.a < b.a;
}Code language: C++ (cpp)Primary Channel and Sort
As a part of the median cut algorithm, you need to be able to determine the primary colour channel being used in a box/bucket. You then sort the colours from highest to lowest according to that channel.
void determine_primary_color_and_sort_box(std::vector<P>& box, std::uint8_t& redRange, std::uint8_t& greenRange, std::uint8_t& blueRange, std::uint8_t& alphaRange)
{
redRange = (*std::max_element(box.begin(), box.end(), red_comp)).r - (*std::min_element(box.begin(), box.end(), red_comp)).r;
greenRange = (*std::max_element(box.begin(), box.end(), green_comp)).g - (*std::min_element(box.begin(), box.end(), green_comp)).g;
blueRange = (*std::max_element(box.begin(), box.end(), blue_comp)).b - (*std::min_element(box.begin(), box.end(), blue_comp)).b;
alphaRange = (*std::max_element(box.begin(), box.end(), alpha_comp)).a - (*std::min_element(box.begin(), box.end(), alpha_comp)).a;
if(redRange >= greenRange && redRange >= blueRange && redRange >= alphaRange)
{
std::sort(box.begin(), box.end(), red_comp);
}
else if(greenRange >= redRange && greenRange >= blueRange && greenRange >= alphaRange)
{
std::sort(box.begin(), box.end(), green_comp);
}
else if(blueRange >= redRange && blueRange >= greenRange && blueRange >= alphaRange)
{
std::sort(box.begin(), box.end(), blue_comp);
}
else
{
std::sort(box.begin(), box.end(), alpha_comp);
}
}Code language: C++ (cpp)To start with we determine what the range, i.e. the highest minus the lowest, value of each channel is. Based on which channel is predominent across the box, we sort it. We pass the ranges back out to the caller.
Generating the Palette, i.e. the Median Cut Algorithm Part 1
Here we define what is essentially the first part of the median-cut algorithm. It’s provided with our source pixels and the number of colours we want to reduce down to. We are returned the palette to use, or the numColors most descriptive colours of the source.
std::vector<P> median_cut_generate_palette(const std::vector<P>& source, const std::uint_fast32_t numColors)
{
typedef std::vector<P> Box;
typedef std::pair<std::uint8_t,Box> RangeBox;
std::vector<RangeBox> boxes;
Box init = source;
boxes.push_back(RangeBox(0,init));
while(boxes.size() < numColors)
{
/* for each box, sort the boxes pixels according to the colour it has the most range in */
for(RangeBox& boxData : boxes)
{
std::uint8_t redRange;
std::uint8_t greenRange;
std::uint8_t blueRange;
std::uint8_t alphaRange;
if(std::get<0>(boxData) == 0)
{
determine_primary_color_and_sort_box(std::get<1>(boxData), redRange, greenRange, blueRange, alphaRange);
if(redRange >= greenRange && redRange >= blueRange && redRange >= alphaRange)
{
std::get<0>(boxData) = redRange;
}
else if(greenRange >= redRange && greenRange >= blueRange && greenRange >= alphaRange)
{
std::get<0>(boxData) = greenRange;
}
else if(blueRange >= redRange && blueRange >= greenRange && blueRange >= alphaRange)
{
std::get<0>(boxData) = blueRange;
}
else
{
std::get<0>(boxData) = alphaRange;
}
}
}
std::sort(boxes.begin(), boxes.end(), [](const RangeBox& a, const RangeBox& b)
{
return std::get<0>(a) < std::get<0>(b);
});
std::vector<RangeBox>::iterator itr = std::prev(boxes.end());
Box biggestBox = std::get<1>(*itr);
boxes.erase(itr);
// the box is sorted already, so split at median
Box splitA(biggestBox.begin(), biggestBox.begin() + biggestBox.size() / 2);
Box splitB(biggestBox.begin() + biggestBox.size() / 2, biggestBox.end());
boxes.push_back(RangeBox(0, splitA));
boxes.push_back(RangeBox(0, splitB));
}
// each box in boxes can be averaged to determine the colour
std::vector<P> palette;
for(const RangeBox& boxData : boxes)
{
Box box = std::get<1>(boxData);
std::uint_fast32_t redAccum = 0;
std::uint_fast32_t greenAccum = 0;
std::uint_fast32_t blueAccum = 0;
std::uint_fast32_t alphaAccum = 0;
std::for_each(box.begin(),box.end(),[&](const P& p)
{
redAccum += p.r;
greenAccum += p.g;
blueAccum += p.b;
alphaAccum += p.a;
});
redAccum /= static_cast<std::uint_fast32_t>(box.size());
greenAccum /= static_cast<std::uint_fast32_t>(box.size());
blueAccum /= static_cast<std::uint_fast32_t>(box.size());
alphaAccum /= static_cast<std::uint_fast32_t>(box.size());
palette.push_back(
{
static_cast<std::uint8_t>(std::min(blueAccum, 255u)),
static_cast<std::uint8_t>(std::min(greenAccum, 255u)),
static_cast<std::uint8_t>(std::min(redAccum, 255u)),
static_cast<std::uint8_t>(std::min(alphaAccum, 255u))
});
}
return palette;
}Code language: C++ (cpp)What we’re doing is generating a new box until we have the number of boxes needed. The number of boxes needed is the number of colours we want.
But when splitting a box we split whichever box is covering the largest range of colours. What this does is that it ensures that if a box contains pixels sorted by red, and that redness covers 75% of the channel range, then it’s split before a blue-sorted box that only covers 10% of the blue-channel range.
By doing this we get the colours that best represent the image.
Median-Cut Part 2
Now we get to the final part.
void birunji::algorithm::median_cut(std::vector<P>& source,
std::vector<P>& destination,
std::uint_fast32_t numColors,
const bool differ,
std::size_t width)
{
std::vector<P> palette = median_cut_generate_palette(source, numColors);
destination = source;
if(!differ)
{
for(P& p : destination)
{
p = *std::min_element(palette.begin(), palette.end(), [p](const P& a, const P& b)
{
return std::sqrt(std::pow(p.r-a.r, 2) +
std::pow(p.g-a.g, 2) +
std::pow(p.b-a.b, 2) +
std::pow(p.a-a.a, 2))
< std::sqrt(std::pow(p.r-b.r, 2) +
std::pow(p.g-b.g, 2) +
std::pow(p.b-b.b, 2) +
std::pow(p.a-b.a, 2));
});
}
}
else
{
/* dither */
}
}Code language: C++ (cpp)We’ll get to the dithering part next. The non-dithering part is a lot simpler, we just use the euclidean distance to determine which colour in the palette we generated, best represents each source pixel.
There are other methods to determine the distance. Euclidean works for my purpoese, by other perceptual methods can produce more accurate results (accurate to what the eyes see, rather than what the data says).
Dithering the Results
If we go down the dithering path above, we would use the following code parts.
First we’ll define a float-based pixel. This is just easier to work with as it gives us more accuracy.
struct PX
{
float r;
float g;
float b;
float a;
};Code language: C++ (cpp)Now we’ll transform the destination into this float range (remember we set the destination to be a copy of the source before).
std::vector<PX> largerRange;
for(const P& p : destination)
{
largerRange.push_back(
{
static_cast<float>(p.r)/255.f,
static_cast<float>(p.g)/255.f,
static_cast<float>(p.b)/255.f,
static_cast<float>(p.a)/255.f
});
}Code language: C++ (cpp)Then we’ll go through each pixel top to bottom, left to right and determine which pixel in the palette best represents the source pixel. Again, same euclidean measurement as before.
for(std::size_t y = 0; y < destination.size() / width; ++y)
{
for(std::size_t x = 0; x < width; ++x)
{
PX oldPixel = largerRange[x + y * width];
P newPixel = *std::min_element(palette.begin(), palette.end(), [oldPixel](const P& a, const P& b)
{
return std::sqrt(std::pow((oldPixel.r) - static_cast<float>(a.r) / 255.f, 2) +
std::pow((oldPixel.g) - static_cast<float>(a.g) / 255.f, 2) +
std::pow((oldPixel.b) - static_cast<float>(a.b) / 255.f, 2) +
std::pow((oldPixel.a) - static_cast<float>(a.a) / 255.f, 2))
< std::sqrt(std::pow((oldPixel.r) - static_cast<float>(b.r) / 255.f, 2) +
std::pow((oldPixel.g) - static_cast<float>(b.g) / 255.f, 2) +
std::pow((oldPixel.b) - static_cast<float>(b.b) / 255.f, 2) +
std::pow((oldPixel.a) - static_cast<float>(b.a) / 255.f, 2));
});
destination[x + y * width] = newPixel;Code language: C++ (cpp)Next we calculate the amount of error between the original pixel and the newly chosen colour from the palette:
float errorR = oldPixel.r - static_cast<float>(newPixel.r) / 255.f;
float errorG = oldPixel.g - static_cast<float>(newPixel.g) / 255.f;
float errorB = oldPixel.b - static_cast<float>(newPixel.b) / 255.f;
float errorA = oldPixel.a - static_cast<float>(newPixel.a) / 255.f;
errorR *= 0.75;
errorG *= 0.75;
errorB *= 0.75;
errorA *= 0.75;Code language: C++ (cpp)Notice that we only carry other 75% of the error. You can adjust this as needed, probably set it as a variable.
Finally we carry the error other according to the Floyd-Steinberg matrix. We’ve complicated things with some bounds checking as well.
if(x < width -1)
{
largerRange[x+1 + y * width].r = std::min(std::max(largerRange[x+1 + y * width].r + errorR / 32 * 5,0.f),1.f);
largerRange[x+1 + y * width].g = std::min(std::max(largerRange[x+1 + y * width].g + errorG / 32 * 5,0.f),1.f);
largerRange[x+1 + y * width].b = std::min(std::max(largerRange[x+1 + y * width].b + errorB / 32 * 5,0.f),1.f);
largerRange[x+1 + y * width].a = std::min(std::max(largerRange[x+1 + y * width].a + errorA / 32 * 5,0.f),1.f);
if(y < destination.size() / width - 1)
{
largerRange[x+1 + (y+1) * width].r = std::min(std::max(largerRange[x+1 + (y+1) * width].r + errorR / 32 * 4,0.f),1.f);
largerRange[x+1 + (y+1) * width].g = std::min(std::max(largerRange[x+1 + (y+1) * width].g + errorG / 32 * 4,0.f),1.f);
largerRange[x+1 + (y+1) * width].b = std::min(std::max(largerRange[x+1 + (y+1) * width].b + errorB / 32 * 4,0.f),1.f);
largerRange[x+1 + (y+1) * width].a = std::min(std::max(largerRange[x+1 + (y+1) * width].a + errorA / 32 * 4,0.f),1.f);
}
if(y < destination.size() / width - 2)
{
largerRange[x+1 + (y+2) * width].r = std::min(std::max(largerRange[x+1 + (y+2) * width].r + errorR / 32 * 2,0.f),1.f);
largerRange[x+1 + (y+2) * width].g = std::min(std::max(largerRange[x+1 + (y+2) * width].g + errorG / 32 * 2,0.f),1.f);
largerRange[x+1 + (y+2) * width].b = std::min(std::max(largerRange[x+1 + (y+2) * width].b + errorB / 32 * 2,0.f),1.f);
largerRange[x+1 + (y+2) * width].a = std::min(std::max(largerRange[x+1 + (y+2) * width].a + errorA / 32 * 2,0.f),1.f);
}
}
if(x < width - 2)
{
largerRange[x+2 + y * width].r = std::min(std::max(largerRange[x+2 + y * width].r + errorR / 32 * 3,0.f),1.f);
largerRange[x+2 + y * width].g = std::min(std::max(largerRange[x+2 + y * width].g + errorG / 32 * 3,0.f),1.f);
largerRange[x+2 + y * width].b = std::min(std::max(largerRange[x+2 + y * width].b + errorB / 32 * 3,0.f),1.f);
largerRange[x+2 + y * width].a = std::min(std::max(largerRange[x+2 + y * width].a + errorA / 32 * 3,0.f),1.f);
if(y < destination.size() / width - 1)
{
largerRange[x+2 + (y+1) * width].r = std::min(std::max(largerRange[x+2 + (y+1) * width].r + errorR / 32 * 2,0.f),1.f);
largerRange[x+2 + (y+1) * width].g = std::min(std::max(largerRange[x+2 + (y+1) * width].g + errorG / 32 * 2,0.f),1.f);
largerRange[x+2 + (y+1) * width].b = std::min(std::max(largerRange[x+2 + (y+1) * width].b + errorB / 32 * 2,0.f),1.f);
largerRange[x+2 + (y+1) * width].a = std::min(std::max(largerRange[x+2 + (y+1) * width].a + errorA / 32 * 2,0.f),1.f);
}
}
if(y < destination.size() / width - 1)
{
largerRange[x + (y+1) * width].r = std::min(std::max(largerRange[x + (y+1) * width].r + errorR / 32 * 5,0.f),1.f);
largerRange[x + (y+1) * width].g = std::min(std::max(largerRange[x + (y+1) * width].g + errorG / 32 * 5,0.f),1.f);
largerRange[x + (y+1) * width].b = std::min(std::max(largerRange[x + (y+1) * width].b + errorB / 32 * 5,0.f),1.f);
largerRange[x + (y+1) * width].a = std::min(std::max(largerRange[x + (y+1) * width].a + errorA / 32 * 5,0.f),1.f);
if(x > 0.f)
{
largerRange[x-1 + (y+1) * width].r = std::min(std::max(largerRange[x-1 + (y+1) * width].r + errorR / 32 * 4,0.f),1.f);
largerRange[x-1 + (y+1) * width].g = std::min(std::max(largerRange[x-1 + (y+1) * width].g + errorG / 32 * 4,0.f),1.f);
largerRange[x-1 + (y+1) * width].b = std::min(std::max(largerRange[x-1 + (y+1) * width].b + errorB / 32 * 4,0.f),1.f);
largerRange[x-1 + (y+1) * width].a = std::min(std::max(largerRange[x-1 + (y+1) * width].a + errorA / 32 * 4,0.f),1.f);
}
if(x > 1)
{
largerRange[x-2 + (y+1) * width].r = std::min(std::max(largerRange[x-2 + (y+1) * width].r + errorR / 32 * 2,0.f),1.f);
largerRange[x-2 + (y+1) * width].g = std::min(std::max(largerRange[x-2 + (y+1) * width].g + errorG / 32 * 2,0.f),1.f);
largerRange[x-2 + (y+1) * width].b = std::min(std::max(largerRange[x-2 + (y+1) * width].b + errorB / 32 * 2,0.f),1.f);
largerRange[x-2 + (y+1) * width].a = std::min(std::max(largerRange[x-2 + (y+1) * width].a + errorA / 32 * 2,0.f),1.f);
}
}
if(y < destination.size() / width - 2)
{
largerRange[x + (y+2) * width].r = std::min(std::max(largerRange[x + (y+2) * width].r + errorR / 32 * 3,0.f),1.f);
largerRange[x + (y+2) * width].g = std::min(std::max(largerRange[x + (y+2) * width].g + errorG / 32 * 3,0.f),1.f);
largerRange[x + (y+2) * width].b = std::min(std::max(largerRange[x + (y+2) * width].b + errorB / 32 * 3,0.f),1.f);
largerRange[x + (y+2) * width].a = std::min(std::max(largerRange[x + (y+2) * width].a + errorA / 32 * 3,0.f),1.f);
if(x > 0.f)
{
largerRange[x-1 + (y+2) * width].r = std::min(std::max(largerRange[x-1 + (y+2) * width].r + errorR / 32 * 2,0.f),1.f);
largerRange[x-1 + (y+2) * width].g = std::min(std::max(largerRange[x-1 + (y+2) * width].g + errorG / 32 * 2,0.f),1.f);
largerRange[x-1 + (y+2) * width].b = std::min(std::max(largerRange[x-1 + (y+2) * width].b + errorB / 32 * 2,0.f),1.f);
largerRange[x-1 + (y+2) * width].a = std::min(std::max(largerRange[x-1 + (y+2) * width].a + errorA / 32 * 2,0.f),1.f);
}
}Code language: C++ (cpp)And that’s it.
Here’s a few different examples:
Here’s a section of a colourful parrot to try and push the boundaries of median cuts ability to represent the original image:
And here it is with a few different levels of cut:
And now 8 colours, with and without dithering.
I hope someone has found this useful to reduce the colours used in their own games.
Hi. Thank you for this article. The algorithm is useful but its implementation has a bug: if you try to reduce colors of an image to a number of colors bigger than the image itself has (for instance, a plain yellow image to 256 colors) the program will crash. Could you please fix it?