Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

Sign up
Here's how it works:
  1. Anybody can ask a question
  2. Anybody can answer
  3. The best answers are voted up and rise to the top

After seeing an example of what is called “Push-pin Art”, I assumed that Mathematica would be able to handle the kind of image processing required to convert an image to this format. I am not well versed in the area of imagine processing: its vocabulary or methods so after several attempts, I decided to get help.

What I would like to do:

  1. Take any image and convert it to a Red-Yellow-Blue-White-Black color space with the ability to designate a desired resolution (width and height of the resulting pixel array) - i.e. to find a minimum acceptable image resolution.
  2. Also, generate a list (2D array) of the colors used in the resulting image e.g. the color names: {{Y,B,R,R},{W,Blk,Y,B},{etc.}}.

When I say RYBWBlk color space I mean that each pixel is forced to only one of these colors, where the eye itself does the blending.

Attempts: [Using Mathematica version 8]

  • In Mathematica, it did not seem possible to create my own indexed color scheme.


  • The RYB… color space is not available to convert to.

  • Posterization did not appear to give me control over the colors chosen, only the number of colors.


  • ColorQuantize seemed like a possible help but no success here.

The following YouTube video and link show Exactly what I would like to be able to accomplish using Mathematica alone.

Photoshop image processing example.

Further explanation.

share|improve this question
up vote 16 down vote accepted

Your problem is essentially colour quantization, and to make it look good you need some form of dithering. Floyd-Steinberg dithering is a classic algorithm that also happens to have nice pseudocode in the Wikipedia article. (Rojo's answer appears to do something like random dithering.)

Note: This functionality probably exists somewhere inside Mathematica already. I'm talking about the ColorQuantize[image, n] function, which has to first (i) compute a good palette with n colours for the image, and then (ii) apply dithering to the image using this palette. Perhaps in a future version, it might be useful to allow the desired palette itself to be passed into ColorQuantize.

Here's a completely direct, extremely inefficient translation of said pseudocode to Mathematica. It can probably be made much faster using Compile and directly poking at the ImageData. If anyone wants to improve it, be my guest: I've marked this Community Wiki.

I'll run it on a small version of the image, because my code is slow.

i = ImageResize[
F7VUQOEH0OJ3L5X.LARGE.jpg"], 100]

Module[{w, h, listArgMin, closestPaletteColor, x, y, oldPixel, 
  newPixel, quantError},
 j = i;
 {w, h} = ImageDimensions[i];
 listArgMin[f_, L_List] := L~Extract~Ordering[f /@ L, 1];
 closestPaletteColor[c_] := 
   SquaredEuclideanDistance[c, #] &, {{1, 0, 0}, {1, 1, 0}, {0, 0, 
     1}, {1, 1, 1}, {0, 0, 0}}];
 For[y = 1, y <= h, y++,
  For[x = 1, x <= w, x++,
   oldPixel = PixelValue[j, {x, y}];
   newPixel = closestPaletteColor[oldPixel];
   quantError = oldPixel - newPixel;
   j = ReplacePixelValue[j, {{x, y} -> newPixel,
      {x + 1, y} -> PixelValue[j, {x + 1, y}] + 7/16 quantError,
      {x - 1, y + 1} -> 
       PixelValue[j, {x - 1, y + 1}] + 3/16 quantError,
      {x + 0, y + 1} -> 
       PixelValue[j, {x + 0, y + 1}] + 5/16 quantError,
      {x + 1, y + 1} -> 
       PixelValue[j, {x + 1, y + 1}] + 1/16 quantError}];]]]

{i, j}

enter image description here

P.S. I adapted Mr.Wizard's listMaxArg function from this previous question.

Here is a compiled version (snappier):

palette = {Black, White, Red, Yellow, Blue};
colors = List @@ ColorConvert[#, "RGB"] & /@ palette;
cf = Compile[{{imdata, _Real, 3}, {colors, _Real, 2}},
       Module[{res, in, new, i, j, w, h, err, data, n},
       data = imdata;
       {h, w} = Most@Dimensions@data;
       res = Array[0. &, {h, w, 3}];
           in = data[[i, j]];
           n = First@Ordering[With[{\[Delta] = in - #}, \[Delta].\[Delta]] & /@ colors, 1]; 
           new = colors[[n]];
           res[[i, j]] = new;
           err = (in - new)/16;
           If[j < w, data[[i, j + 1]] += 7 *err];
           If[j > 1 && i < h, data[[i + 1, j - 1]] += 3 *err];
           If[i < h, data[[i + 1, j]] += 5 *err];
           If[i < h && j < w, data[[i + 1, j + 1]] += err],
         {i, h}, {j, w}];
       res], CompilationTarget -> "C"];
Image@cf[ImageData[resized], colors]

Full size result:

enter image description here

share|improve this answer
@Matthias: Sweet, thanks! – Rahul Dec 15 '12 at 2:46
We should probably linearize the colour space before doing any dithering, too... – Rahul Dec 15 '12 at 2:47
Very nice! This looks fantastic – R. M. Dec 15 '12 at 2:50
@RahulNarain probably not necessary, but thanks for the reminder about this other question – Matthias Odisio Dec 15 '12 at 2:52
Very good +1. I was expecting to get a way better answer than mine a few minutes afterwards, and you didn't disappoint – Rojolalalalalalalalalalalalala Dec 15 '12 at 2:55


The image

Let's work on the lady of the link

im = Import[""];
imSmall = im~ImageResize~200

Some palettes

Let's represent a palette as an image in the desired color space.

colors2palette[colors_List, colSpace_?Image`ColorSpaceQ] := 
  Image[{Flatten[List @@ ColorConvert[#, colSpace]] & /@ colors}, 
   ColorSpace -> colSpace];
paletteData = First~Composition~ImageData;


(paletteOP = colors2palette[{Red, Yellow, Blue, White, Black}, "RGB"]) // Framed
(paletteBW = colors2palette[{Black, White}, "Grayscale"]) // Framed

paletteOP is the palette proposed in the question. {Red, Yellow, Blue, White, Black} gives {RGBColor[1, 0, 0], RGBColor[1, 1, 0], RGBColor[0, 0, 1], GrayLevel[1], GrayLevel[0]}. There doesn't seem to be a way to blend up a color with more green than red.


The most basic approach to quantizing an image to a palette is to round each pixel to its closest palette color, with some metric in the color space. This implies a function appliead to each pixel separately, point to point, without knowledge of the neighbouring pixels.

Let's look at the result

(* Convert the image to the proper color space *)

Options[thresholding] = {"DistanceFunction" -> Automatic};

thresholding[im_Image, palette_Image, op : OptionsPattern[]] /; 
   ImageColorSpace[im] =!= ImageColorSpace@palette :=
    thresholding[ColorConvert[im, ImageColorSpace@palette], palette, 

thresholding[im_Image, palette_Image, OptionsPattern[]] :=
    DistanceFunction -> OptionValue["DistanceFunction"]]~ImageApply~im;


thresholding[imSmall, #] & /@ {paletteOP, paletteBW} // FlipView

Mathematica graphics

Ugly: if a region in the image is of a certain color more or less between two palette colors, all the region will be rounded off to the same color.


Now, here comes something weird an interesting. The problems with simple thresholding can be reduced by adding noise to our image.

randomDithering[im_Image, palette_Image, \[Sigma] : (_?NumericQ) : 0.15] :=
  randomDithering[im, palette, NormalDistribution[0, \[Sigma]]];

randomDithering[im_Image, palette_Image, 
        dist_?Statistics`Library`UnivariateDistributionQ] :=
     Reverse@ImageDimensions@im~Append~ImageChannels@im], palette]


randomDithering[imSmall, #] & /@ {paletteOP, paletteBW} // FlipView

Mathematica graphics

Every time you run the code you'll get a slightly different result with similar overall look.


This involves replacing every pixel by a matrix of pixels of given size. So, if the overall resolution is maintained, say by previously reducing the size of the image, then the effective resolution is reduced.

A quite slow implementation based on linear programming follows

getCoordsFromColor[xtra_, paletteData : {_?NumericQ ..}] :=
  getCoordsFromColor[xtra, List /@ paletteData];

g_getCoordsFromColor[Except[_List, i_]] := g[{i}];

getCoordsFromColor[xtra_Integer?Positive, paletteData_][color_] :=
  Append[paletteData\[Transpose], 1~ConstantArray~Length@paletteData],
  Transpose@{(2 xtra + 1)^2 color~Append~1, 
    1~ConstantArray~Length@color~Append~0}, 0`, Integers]

getCoordsFromImage[im_Image, xtra_Integer?Positive, palette_][im_Image] := 
 Map[getCoordsFromColor[xtra, palette], ImageData[im], {2}]

This takes that matrix and rebuilds it,randomly sorting the pixels in each submatrix

rebuild[paletteData_][coords_] := Block[{RGBColor = List,
     GrayLevel = ConstantArray[#, 3] &,
     side = Sqrt@Total@coords[[1, 1]]},
         ConstantArray @@@ Transpose@{paletteData, coord}, 1], side]],
     coords, {2}]] // ArrayFlatten // Image

patterning[im_Image, palette_Image, side_Integer?Positive] := 
 With[{smallerIm = 
    ImageResize[im, Scaled[1/(2 side + 1)]]~ColorConvert~
    getCoordsFromImage[smallerIm, side, paletteData@palette]][smallerIm]


patterning[imSmall, #, 1] & /@ {paletteOP, paletteBW} // FlipView

Mathematica graphics


The general idea of these algorithms is best illustrated with a black and white palette.

  • Create a threshold matrix, of, say, 4*4, with values between 0 and 1

  • Partition the image into small matrices of the same 4*4 size

  • Compare each image partition's data values with the threshold matrix. If it is higher, the output pixel will be white; lower, black

The name of this method comes from the fact that given an image of a plain color, pixels are turned on one by one as the color gets darker and surpasses the different thresholds. Take a look at this demonstration.

This is a very smart and neat idea since it breaks a little bit with the compromise between resolution and number of colors you can represent that patterning and regular thresholding seemed to impose. At the same time, the output is deterministic, and the computation is fast because in the end each pixel is being thresholded on its own.

Options[orderedDithering] = {"ThresholdMatrixGenerator" -> 
    sparseThresholdMatrix, "Grayscale" -> True};

orderedDithering[im_Image, side_Integer, op : OptionsPattern[]] :=
  orderedDitheringPvt[im, side, op];

orderedDitheringPvt[im_Image, side_Integer, op : OptionsPattern[]] /; 
   OptionValue[orderedDithering, {op}, "Grayscale"] :=
        With[{gim = ColorConvert[im, "Grayscale"]},
   UnitStep[ImageData@gim - 
      OptionValue[orderedDithering, {op}, "ThresholdMatrixGenerator"][
       side, ImageDimensions@gim]] // Image];

The thresholding matrices can be built according to your requirements. Old printers needing to use only black and white to print images, may have preferred to put the ink all together, and "turn the pixels on" in clusters.

tileMatrix[mat_?MatrixQ, {dimensions : Repeated[_Integer?Positive, {2}]}] := 
   ConstantArray[mat, Ceiling[{dimensions}/Dimensions@mat]]~Take~dimensions

This defines a clustered thresholding matrix generator

   dimensions : {_Integer?Positive, _Integer?Positive}] := 
  tileMatrix[SpiralMatrix[Reverse@Range[side^2]] // N // Rescale,

(* Folds the list in a spiral to build a matrix *)
SpiralMatrix[l_List] /; OddQ@Sqrt@Length@l := spiralMatrix[l];
SpiralMatrix[l_List] /; EvenQ@Sqrt@Length@l := 
  spiralMatrix[l]~Drop~Sequence[-1, 1]; 

(* Thanks @Xerxes for the code and thanks a lot @belisarius for the untimely help *)
spiralMatrix[l_List] := 
 Module[{dim = 2 Ceiling[(Sqrt[Length[l]] - 1)/2] + 1, x0}, 
  x0 = Floor[dim/2]; # /. x_?NumericQ :>
      PadRight[l, dim^2][[Position[Sort[Flatten[#]], x][[1, 1]]]] &@
   Array[2 Norm[{##}, \[Infinity]] + 
      Mod[(ArcTan[##] - ArcTan[x0, x0 + 1])/\[Pi], 2] &, {dim, 
     dim}, {-x0 - 1.*^-6, -x0}]]

We can plot it

plotThresholdMatrix[mat_?MatrixQ] := 
 ArrayPlot[mat, ImageSize -> Small]


plotThresholdMatrix@clusteredThresholdMatrix[8, 8 {1, 1}]

Mathematica graphics

However, for other applications you may prefer the visuals of a sparser matrix.

This defines a sparse thresholding matrix generator

    dimensions : {_Integer?Positive, _Integer?Positive}] /;
   Log2@side \[Element] Integers :=
  tileMatrix[Nest[expandMatrix, {{1}}, Log2@side] // N // Rescale,

expandMatrix[mat_] := With[{numElms = Times @@ Dimensions[mat],
   rotationSequence = {{0, 0}, {1, 1}, {0, 1}, {1, 0}}},
  Table[Upsample[mat + i numElms, {2, 2}], {i, 0, 3}] /. mats_ :>
    Total@MapThread[RotateRight, {mats, rotationSequence}]

plotThresholdMatrix@sparseThresholdMatrix[8, 8 {1, 1}]

Mathematica graphics

Let's try them

orderedDithering[imSmall, 4, 
    "ThresholdMatrixGenerator" -> #] & /@ {sparseThresholdMatrix, 
   clusteredThresholdMatrix} // FlipView

Mathematica graphics

Extending this to arbitrary palettes does not seem automatic to me. But we can see some colors if we apply it to each RGB channel separately

orderedDitheringPvt[im_Image, side_Integer, op : OptionsPattern[]] :=
  orderedDitheringPvt[#, side, "Grayscale" -> True, op] & /@ 
   ColorSeparate@ColorConvert[im, "RGBColor"], "RGBColor"]


orderedDithering[imSmall, 4, "ThresholdMatrixGenerator" -> #, 
    "Grayscale" -> False] & /@ {sparseThresholdMatrix, 
   clusteredThresholdMatrix} // FlipView

Mathematica graphics

This would be using the following palette

Image[{Tuples[{1, 0}~ConstantArray~3]}, ColorSpace -> "RGB"] // Framed

Mathematica graphics

Error diffusion

These algorithms are procedural. The idea is to sweep the image, say, in "reading direction" (left to right, top to down), thresholding each pixel and spreading the error among the neighbouring unprocessed pixels.

The Floyd-Steinberg algorithm implemented in the other great answer is an example of this. In particular, it uses the following matrix to propagate the error. The asterisk references the current pixel being processed.

$$M=\frac{1}{16}\left( \begin{array}{cc} & * & 7 \\ 3 & 5 & 1\\ \end{array} \right)$$

For completeness, here goes the stolen code from the other answer

cfFloydSteinberg = 
  Compile[{{imdata, _Real, 3}, {colors, _Real, 2}}, 
   Module[{res, in, new, i, j, w, h, err, data, n}, data = imdata;
    {h, w} = Most@Dimensions@data;
    res = Array[0. &, {h, w, 3}];
    Do[in = data[[i, j]];
     n = First@Ordering[
        With[{\[Delta] = in - #}, \[Delta].\[Delta]] & /@ colors, 1];
     new = colors[[n]];
     res[[i, j]] = new;
     err = (in - new)/16;
     If[j < w, data[[i, j + 1]] += 7*err];
     If[j > 1 && i < h, data[[i + 1, j - 1]] += 3*err];
     If[i < h, data[[i + 1, j]] += 5*err];
     If[i < h && j < w, data[[i + 1, j + 1]] += err], {i, h}, {j, w}];

The interface

floyedSteinbergDithering[im_Image, palette_] := 
 cfFloydSteinberg[ImageData@im, paletteData@palette] // Image

And this would give the output

floyedSteinbergDithering[imSmall, paletteOP]
share|improve this answer

Your Answer


By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.