Preliminaries
The image
Let's work on the lady of the link
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_?Image`ColorSpaceQ] :=
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[1], GrayLevel[0]}. 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_?Statistics`Library`UnivariateDistributionQ] :=
thresholding[
im~ImageAdd~Image@RandomVariate[dist,
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 :>
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]
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_ :>
Total@MapThread[RotateRight, {mats, rotationSequence}]
]
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]
