# Push-pin Art – Convert image to R-Y-B-W-Blk color space

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.

Table[ColorData[i,"Image"],{i,ColorData["Indexed"]}]//TableForm

• 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.

ImageEffect[Import["https://i.stack.imgur.com/wtgxH.jpg"],{"Posterization",8}]

• 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.

http://www.instructables.com/id/How-to-Make-a-Push-Pin-Portrait/

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[
Import["http://cdn.instructables.com/F7V/UQOE/H0OJ3L5X/\
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_] :=
listArgMin[
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} 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}];
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}];
res], CompilationTarget -> "C"];
Image@cf[ImageData[resized], colors]


Full size result: • @Matthias: Sweet, thanks!
– user484
Dec 15, 2012 at 2:46
• We should probably linearize the colour space before doing any dithering, too...
– user484
Dec 15, 2012 at 2:47
• Very nice! This looks fantastic
– rm -rf
Dec 15, 2012 at 2:50
• @RahulNarain probably not necessary, but thanks for the reminder about this other question Dec 15, 2012 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
– Rojo
Dec 15, 2012 at 2:55

# Preliminaries

## The image

im = Import["http://cdn.instructables.com/F7V/UQOE/H0OJ3L5X/F7VUQOEH0OJ3L5X.LARGE.jpg"];
imSmall = im~ImageResize~200 ## Some palettes

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

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


So

(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, GrayLevel}. There doesn't seem to be a way to blend up a color with more green than red.

# Thresholding

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,
op];

thresholding[im_Image, palette_Image, OptionsPattern[]] :=
First~Composition~Nearest[paletteData@palette,
DistanceFunction -> OptionValue["DistanceFunction"]]~ImageApply~im;


So

thresholding[imSmall, #] & /@ {paletteOP, paletteBW} // FlipView 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.

# Random

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_?StatisticsLibraryUnivariateDistributionQ] :=
thresholding[
Reverse@ImageDimensions@im~Append~ImageChannels@im], palette]


So

randomDithering[imSmall, #] & /@ {paletteOP, paletteBW} // FlipView Every time you run the code you'll get a slightly different result with similar overall look.

# Patterning

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_] :=
LinearProgramming[paletteData~Total~{2},
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]]},
Map[Function[coord,
Partition[RandomSample@Flatten[
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~
ImageColorSpace@palette},
Composition[rebuild[paletteData@palette],
getCoordsFromImage[smallerIm, side, paletteData@palette]][smallerIm]
]


So

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

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}]}] :=
ArrayFlatten@
ConstantArray[mat, Ceiling[{dimensions}/Dimensions@mat]]~Take~dimensions


This defines a clustered thresholding matrix generator

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

(* 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 :>
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]


So

plotThresholdMatrix@clusteredThresholdMatrix[8, 8 {1, 1}] However, for other applications you may prefer the visuals of a sparser matrix.

This defines a sparse thresholding matrix generator

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

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_ :>
]

plotThresholdMatrix@sparseThresholdMatrix[8, 8 {1, 1}] Let's try them

orderedDithering[imSmall, 4,
"ThresholdMatrixGenerator" -> #] & /@ {sparseThresholdMatrix,
clusteredThresholdMatrix} // FlipView 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[]] :=
ColorCombine[
orderedDitheringPvt[#, side, "Grayscale" -> True, op] & /@
ColorSeparate@ColorConvert[im, "RGBColor"], "RGBColor"]


So

orderedDithering[imSmall, 4, "ThresholdMatrixGenerator" -> #,
"Grayscale" -> False] & /@ {sparseThresholdMatrix,
clusteredThresholdMatrix} // FlipView This would be using the following palette

Image[{Tuples[{1, 0}~ConstantArray~3]}, ColorSpace -> "RGB"] // Framed # 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}];
res]];


The interface

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


And this would give the output

floyedSteinbergDithering[imSmall, paletteOP]