Artistic image vectorization



The question is how can we use Mathematica to create vectorized versions of low-resolution images? The goal is to get an image suitable for quality printing at any resolution.

Since "true" vectorization performed by various specialized software is a tough problem, I suggest to consider "artistic" approaches, which produce inexact, but beautiful version of the original. Colored implementations are highly appreciated.


Posted 2012-07-19T22:33:47.067

Reputation: 12 161

1This is also built into Inkscape: load the rose image (e.g.), go to Path > Trace Bitmap > Mode > Multiple scans and select (e.g.) color. – Jens – 2012-07-20T00:07:42.780

By "vectorization" he probably meant that in Mathematica image data will be put into graphics primitives. I think the keyword here is "Artistic" though ;-) - meaning you can play with those graphics primitives. – Vitaliy Kaurov – 2012-07-20T01:37:16.000


A small note:

Now in Mathematica 11.1, ImageGraphics can convert an Image object into a scalable vector Graphics.

Documentation is Here.

– Wjx – 2017-03-26T08:46:49.593



This vectorisation attempts to represent the image with coloured triangles. The code selects a user defined number of sample points, with the selection weighted according to the image gradient, to obtain finer sampling in more detailed regions of the image. I use ListPlot3D to triangulate the sample points into a set of polygons - there is probably a neater way. The output from ListPlot3D is stripped of the third dimension and VertexColors are applied to the polygons based on the image colour at the sample points.



img = ImageResize[ExampleData[{"TestImage", "Lena"}], 200];
vectorise[img, #] & /@ {100, 1000, 10000}

enter image description here

Having vectorised the image we can do silly things with the graphics:

Export["fragmentlena.gif",Table[MapAt[#+rand (1-Cos[t])^2&,pic,{1,1,1}],{t,2\[Pi]/40,2\[Pi],2\[Pi]/40}],ImageSize->200,"DisplayDurations"->0.05,AnimationRepetitions->Infinity]

enter image description here

Simon Woods

Posted 2012-07-19T22:33:47.067

Reputation: 81 905

1A big +1! Gradient polygon coloring gives an extraordinary result, and the speed is impressive. This is the most profound approach so far, IMO. – faleichik – 2012-07-22T18:09:26.960


Let's get a low-res image:

enter image description here

And put in in gray-scale mode:

gimg = ColorConvert[ImageResize[
        Import[""], 300], "Grayscale"];

Now extract the image data (pixel values) together with pixel indexes:

data = MapIndexed[Append[#2, #1] &, ImageData[gimg], {2}];

I, of course, couldn't pass on Voronoi styling. We will add random noise to perfectly integer pixel coordinates and make a mosaic animation:

Table[Rotate[ListDensityPlot[(MapThread[Append, {3 RandomReal[{-1, 1}, {Length[#], 2}], 
      ConstantArray[0, Length[#]]}] + #) &@Flatten[Transpose@data, 1][[1 ;; -1 ;; 15]], 
    InterpolationOrder -> 0, ColorFunction -> "GrayTones", 
    BoundaryStyle -> Directive[Black, Opacity[.2]], Frame -> False, 
    PlotRangePadding -> 0, AspectRatio -> Automatic, ImageSize -> 600], -Pi/2], {10}];

Export["test.gif", %]

enter image description here

And various outlandish coloring

           Append, {3 RandomReal[{-1, 1}, {Length[#], 2}], 
            ConstantArray[0, Length[#]]}] + #) &@
       Flatten[Transpose@data, 1][[1 ;; -1 ;; 45]], 
      InterpolationOrder -> 0, ColorFunction -> #, 
      BoundaryStyle -> Directive[Black, Opacity[.2]], Frame -> False, 
      PlotRangePadding -> 0, AspectRatio -> Automatic, 
      ImageSize -> 300], -Pi/2] & /@ {"CherryTones", "CoffeeTones", 
    "DarkRainbow", "DeepSeaColors", "PlumColors", "Rainbow", 
    "StarryNightColors", "SunsetColors", "ValentineTones"}, 3], 
 Spacings -> 0] 

enter image description here

If we fix the noise sampling with SeedRandom and change only magnitude of the noise, we can create a sort of order-from-chaos appearance effect:

id = ParallelTable[Rotate[ListDensityPlot[(MapThread[
          Append, {SeedRandom[1]; 
           200 (1 - st^(1/8)) RandomReal[{-1, 1}, {Length[#], 2}], 
           ConstantArray[0, Length[#]]}] + #) &@
      Flatten[Transpose@data, 1][[1 ;; -1 ;; 15]], 
     InterpolationOrder -> 0, ColorFunction -> GrayLevel, 
     BoundaryStyle -> Opacity[.1], Frame -> False, 
     PlotRangePadding -> 0, AspectRatio -> Automatic, 
     ImageSize -> 350], -Pi/2], {st, 0.2, 1, .05}];

idd = id~Join~Table[id[[-1]], {7}];

Export["appear.gif", idd, ImageSize -> 350]

enter image description here

Vitaliy Kaurov

Posted 2012-07-19T22:33:47.067

Reputation: 66 672

3This is a viable alternative to my interpretation of the question. +1 – Mr.Wizard – 2012-07-20T06:07:32.010

3love your Voronoism! – cormullion – 2012-07-20T08:25:25.077

Great job! Vitaliy, is it possible to take cell colors from the original color image?... Also, I would suggest to use data = MapIndexed[#2~Join~{#1} &, ImageData[gimg], {2}]; – faleichik – 2012-07-21T09:50:34.953

@faleichik Thanks! Yes it is possible "to take cell colors from the original color image". ListDensityPlot produces graphics based on Polygon all of which can be extracted and put together with RGBColor function that was mapped on original ImageData. It's a bit tedious though ;) And good tip on Join - I use it myself when speed counts. See this on details how to extract polygons:

– Vitaliy Kaurov – 2012-07-24T08:22:34.450

2@VitaliyKaurov very nice! Suggestion: For copy&paste&evaluate purposes it would be more comfortable if your image were imported via some URL. – Yves Klett – 2012-07-24T09:04:31.013

@YvesKlett thank you! I am not sure I understand what do you mean by "image were imported via some URL" and what it got to do with Copy/Paste - could you please explain? All images on MSE have their URLs and are downloadable. – Vitaliy Kaurov – 2012-07-24T20:41:42.137

@VitaliyKaurov Meaning acutally importing it in the code like in e.g. @Szabolcs answer (img = Import[""]) – Yves Klett – 2012-07-24T21:01:43.750

@YvesKlett The code should work with any image. Are you interested in the used image in particular? - I do not remember where I got that image. – Vitaliy Kaurov – 2012-07-24T21:30:11.040

The last one is sort of creepy. :) – rcollyer – 2012-07-25T01:27:25.253

Sorry, for leaving the code so late but that ParallelTable construct you use to extract pixel and pixel-indices is, aehmm the horror;-) It takes forever (23sec) here while the obvious MapIndexed[Append[#2, #1] &, ImageData[gimg], {2}] takes 0.1sec! – halirutan – 2012-10-02T16:54:43.610


Let me start from an approach which I believe has its own name, but unfortunately I don't know it. The idea is to generate circles whith radii depending on the intensity of corresponding and surrounding pixels.

img = Import["ExampleData/rose.gif"];
bubbled[img_, r_, delta_, rmax_] := 
  Block[{ker, data, radii, thresh = 0.99},
   ker = N[#/Total@Total@#] &@DiskMatrix[r];
   data = Map[Mean, ImageData[ImageConvolve[img, ker]], {2}];
   radii = Partition[data, {1, 1}, delta] /. {{a_?NumberQ}} -> a;
    MapIndexed[If[#1 < thresh, Disk[#2, rmax Max[1 - #1, 0]]] &, 
     Reverse /@ (Transpose@radii), {2}]
is = 360;
 Row[Show[#, ImageSize -> is] & /@ {img, bubbled[img, r, d, rmax]}],
 {{r, 1, "Smooth radius"}, 1, 10, 1, ControlType -> Setter},
 {{d, 3, "Offset"}, 1, 10, 1, ControlType -> Setter},
 {{rmax, Sqrt[2.]/2, "Maximum radius"}, 0.5, 1, 0.05}

enter image description here

Some more examples:

enter image description here

enter image description here

One can then save the image for printing:

gfx = bubbled[img, 2, 3, 1]
Export[NotebookDirectory[] <> "gfx.pdf", gfx]


Posted 2012-07-19T22:33:47.067

Reputation: 12 161

One thing which I don't like is radii = Partition[data, {1, 1}, delta] /. {{a_?NumberQ}} -> a; I guess Flatten must be used here, but levelspecs are still a mystery for me. – faleichik – 2012-07-19T22:53:11.190

9 – Brett Champion – 2012-07-19T22:53:54.043

1@faleichik nice answer! Flatten[Partition[data, {1, 1}, delta], {{1}, {2, 3, 4}}] is what you're looking for, I think. – Oleksandr R. – 2012-07-20T12:53:40.960

Brett, thank you, I didn't suspect it to be so simple :-) Many thanks to Oleksandr as well! I would never come to this myself... Need to make myself study the corresponding thread here on – faleichik – 2012-07-21T09:23:52.083


Take an image:


Grayscale, flip, and smooth it:

j = CurvatureFlowFilter[ImageReflect[ColorConvert[i, "Grayscale"]]]


Then let ListContourPlot do all the work!

lcp = ListContourPlot[ImageData[j], Frame -> False, Contours -> Table[p, {p, 0, 1, 0.025}],
  ColorFunction -> "GrayTones", ContourStyle -> None, PlotRange -> {0, 1}]


Make it as big as you please:

Export["sw.png", lcp, ImageSize -> 1280]


Posted 2012-07-19T22:33:47.067

Reputation: 13 711

1You can use the DataReversed option to ImageData and avoid reflecting the image. (+1 BTW) – Szabolcs – 2012-07-20T23:52:17.647


The following approach was suggested by cormullion in his answer, namely it is "spiral engraving". The idea is to approximate the image by short lines arranged in a spiral. The main issue is the lines lengths: they should depend on the spiral curvature and on the image details. I decided to implement this by extracting points from ParametricPlot and then spliting the lines which are longer than maxlen. Then each line is drawn with fradient colors from the corresponding pixels of the original image.

img=enter image description here

ptsSpiral[d_, tmax_, maxlen_] := Block[{pts, prev},
   d - the spiral "density",
   tmax - determines spiral length,
   maxlen - maximum line length in resulting splitting 
   pts = ParametricPlot[d # {Sin@#, Cos@#} &[t], {t, 0, tmax}, 
      PlotPoints -> 100][[1, 1, 3, 2, 1]];
   prev = pts[[1]];
        If[(len = Norm[# - prev]) <= maxlen, Sow@#,
         Sow /@ 
            Range[prev, #, (# - prev)/Ceiling[len/maxlen]]]
         ]; prev = #) &
      , pts]
     ][[2, 1]]

spiralify[img_, pts_, th_] := Block[{c},
   c = ImageDimensions@img/2;
     Line[pts, VertexColors -> ImageValue[img, c + # & /@ pts]]}]

img = ImageResize[img, {Automatic, 200}]
pts = ptsSpiral[0.25, 380, 2];
Show[spiralify[img, pts, 2], ImageSize -> 400]

enter image description here

Next thing we can try is to play with thickness rather than color: we draw black line with thickness depending on the intensity of corresponding pixels.

engrave[img_, pts_, maxth_] := Block[{c},
   c = (ImageDimensions@img)/2;
      "Round"], {AbsoluteThickness[
         maxth (1 - Mean@ImageValue[img, #])], Line@#} & /@ 
      Partition[# + c & /@ pts, 2, 1]}]
pts = ptsSpiral[0.25, 399, 1];
engrave[ColorConvert[img, "Grayscale"], pts, 4]

enter image description here

And one more image I love the most:

enter image description here

It was made from this image.

And finally I must apologize for the second self-answer here...


Posted 2012-07-19T22:33:47.067

Reputation: 12 161

4Why apologize for a self-answer, especially one inspired by another answer? SE explicitly encourages self-answers, so don't worry about it. +1 – rcollyer – 2012-07-25T01:35:02.000

This is lovely! – Yves Klett – 2012-07-25T06:30:39.847

That's really excellent - the last one is especially effective! – cormullion – 2012-07-25T06:33:43.603

1In Mathematica 12.1.1, ParametricPlot in ptsSpiral needs to be extracted with [[1, 1, 1, 3, 1, 2, 1]] instead. – Taiki – 2020-10-19T07:42:48.043

1In addition, ImageValue in engrave needs to be replaced by Map[First]@ImageValue. – Taiki – 2020-10-19T08:31:36.000


This is from an old notebook of mine. The girl from Ipanema, composed of 35,000 points. (Note true pointillism, perhaps someone can do that one as well.) Sorry about the messy code and the slow processing.

The idea is to randomly scatter disks in a rectangle, and colour them according to the corresponding points of the photo. More points with a higher density are used in regions of high detail or sharp transitions, fewer points elsewhere (to keep their number down).

img = Import[""]

The Girl from Ipanema

etf = EntropyFilter[img, 12] // ImageAdjust
sdf = ColorConvert[StandardDeviationFilter[img, 5], "GrayScale"] // ImageAdjust
map = ImageAdd[sdf, etf] // ImageAdjust

Point size map

mapdata = ImageData[map];
data = ImageData[img];

{w, h} = ImageDimensions[img];

ch = RandomChoice[
         (Flatten[mapdata] + 0.1)^1.7 -> Join @@ Table[{i, j}, {i, h}, {j, w}], 

spots = {data[[#1, #2]], {#2, -#1}, 15 (1.1 - mapdata[[#1, #2]])^1.8} & @@@ ch;
spots = Reverse@SortBy[spots, Last];

Graphics[{RGBColor[#1], Disk[#2, #3]} & @@@ spots,
 Background -> GrayLevel[0.75],
 PlotRange -> {{1, w}, {1, -h}}

Mathematica graphics

Another example (unfortunately I lost the original source photo and it's not as easy to google up as the other one above):

Mathematica graphics


Posted 2012-07-19T22:33:47.067

Reputation: 213 047

@Szabolcs - BTW, the word is search by image. There's a Firefox plug-in which I use almost daily.

– stevenvh – 2015-04-19T07:30:07.177

7Your Google-fu is weak. Or maybe it only finds bikini... – Mr.Wizard – 2012-07-20T09:01:28.043

@Mr.Wizard Did you search for the image itself? I forgot completely about that. I did this some 3 years ago, and I remember trying to find the big cat picture again, but I didn't manage at that time. – Szabolcs – 2012-07-20T09:18:37.450

7I think I chose the wrong image for my answer ... :) – cormullion – 2012-07-20T10:01:51.760

+1 for that summer feeling (and now I have to get the tune out of my head: na-nanana-na-na-nanana...) – Yves Klett – 2012-07-20T10:50:17.307

1@cormullion Not too late to change it though:) – Ajasja – 2012-07-20T11:46:53.333

1@cormullion I´d second the motion to change ;-) – Yves Klett – 2012-07-20T12:15:07.717


I enjoy representing colour images with tilings, preferably non-periodic ones. Alternatively, one can use overlapping tiles to give a slight 3D effect. Here is some example code.

TranslateObject[p_, {x_, y_}] := Map[{x, y} + # &, p, {2}]

HouseHexPolygon[s_, 0] := 
   Polygon[s*{{1/2,-1}, {3/2,0}, {3/2,1}, {-1/2,1}, {-3/2,0}, {-3/2,-1}}]
HouseHexPolygon[s_, 1] := 
   Polygon[s*{{1,1/2}, {0, 3/2}, {-1,3/2}, {-1,-1/2}, {0,-3/2}, {1,-3/2}}]

HouseHexGrid[s_, imax_, jmax_] :=
   Block[{imod, jmod, k, m},
         imod = Mod[2 i, 5];
         jmod = Mod[2 i + 3, 5];
         k = Range[0, Floor[(jmax - imod)/5]];
         m = Range[0, Floor[(jmax - jmod)/5]];
         {i, 0, imax}], 2]]

HouseHexGridPerturbed[s_, h_, v_, r_] :=
   Block[{poly=Map[Round[#,10.^-10]&, HouseHexGrid[N[s],h,v], {2}], pts, len, rules},
      pts = DeleteDuplicates[Flatten[poly[[All, 1]], 1]];
     len = Length[pts];
     rules = Dispatch[Thread[pts -> Range[len]]];
        pts + RandomReal[{-r, r}, {len, 2}],
        RandomSample[poly] /. rules]

HouseHexColourCentroid[image_Image, s_, t_, opts___] :=
   Block[{d, dim, g, centroids, colours, rr, gg, bb, greys},
      d = Reverse[ImageData[image]]; (* Reverse to get bottom of image on bottom *)
      dim = Dimensions[d]; (* dim[[1]] is number of rows, dim[[2]] is number of cols *)
      g = HouseHexGridPerturbed[s, Floor[dim[[2]]/s-3], Floor[dim[[1]]/s-3], 0];
      centroids = Apply[Mean[g[[1, #]]] &, g[[2]], 1];
      colours = Apply[RGBColor, Extract[d, Map[Reverse, Round[centroids+s/2]]], 1];
      greys = colours /. RGBColor[rr_,gg_,bb_] -> 0.299 rr + 0.587 gg + 0.114 bb;
      g[[2]] = Transpose[{colours, g[[2]]}][[Ordering[t*greys]]];
      Graphics[{EdgeForm[{Thickness[0.0003], Black}], 
                GraphicsComplex[g[[1]], g[[2]]]},

In the following (low-resolution) example of the Helix nebula, the brighter filaments appear in front of the fainter ones, and the bright stars appear layered.

   Import[""], 3, 1]

Helix Nebula


Posted 2012-07-19T22:33:47.067

Reputation: 14 269

Nice, art lies in concealing art. – Narasimham – 2016-06-23T08:47:58.830


Well, this is really, really inefficient, but there you go:

comp = WatershedComponents[img, Method -> "Rainfall"];

imdat = ImageData[img];

colorrectangle[n_] := 
 With[{pos = Position[comp, n]}, {RGBColor[
    Mean[imdat[[Sequence @@ #]] & /@ pos]], 
   Rectangle[{#[[2]], -#[[1]]}] & /@ pos}]

newimg = Graphics[
   colorrectangle /@ 
    Range[Length[comp // Flatten // DeleteDuplicates]]];

Show[#, ImageSize -> 200, PlotRangePadding -> 0] & /@ {img, newimg}

Mathematica graphics

WatershedComponents returns the segements which are brute-forced to Rectangles and colored with the mean of the segment pixels in the image.

More effective renderingwise:

colorep[n_] := 
 With[{pos = Position[comp, n]}, 
  n -> RGBColor[Mean[imdat[[Sequence @@ #]] & /@ pos]]]

reps = Flatten[
    colorep /@ Range[Length[comp // Flatten // DeleteDuplicates]]] // 

ArrayPlot[comp, ColorRules -> reps, PlotRangePadding -> 0]

Mathematica graphics

It would be even better to get a vector outline of the segments and represent them as polygons. Then you´d get away from the 8-bit Nintendo look...

Yves Klett

Posted 2012-07-19T22:33:47.067

Reputation: 14 743

Do you know about this?

– Vitaliy Kaurov – 2012-07-20T14:49:42.607

@VitaliyKaurov I knew it was somewhere ;-) Started on a vectorization idea but did not get to the end (yet). Still thought it might be worth putting it here in (potentially) intermediate state. – Yves Klett – 2012-07-22T18:13:15.580


In this attempt, each row of the image is drawn as a polygon, like a mountain range, with the peaks and troughs controlled by the pixel values of the image data. These polygons are (clumsily) constructed and stacked one of top of the other. The idea is echoing the old artist engravers of earlier centuries, who used the varying thickness of engraved lines to suggest tone. Somewhere I've got an abandoned attempt at a spiral engraving thing like the work of some guy called Mellan.

i = enter image description here

im = ImagePad[i, 1, White]
{width, height} = ImageDimensions[im];
pixelData = Reverse[ImageData[im]];

rowToLine[r_, c_, amplitude_] := {
   firstPos = pixelData[[r]][[1]];
   (* add last point to complete polygon *)
       (* x coords are 1 to width, 
          y coords are row height + pixel value * amplitude *)
        Range[1, c], r + (amplitude * pixelData[[r]])], 2]],
     {c, r + (amplitude * firstPos)}, 
     {1, -1}]}};

Manipulate[g = Graphics[{Opacity[opacity],
    Table[{rowToLine[row, width, x]}, {row, 1, height, step} ]}, 
   ImageSize -> 500],
  {x, 0, 5, 1}, 
 {{thickness, 0.1}, 0.1, 2},
 {step, 1, 12, 1}, 
 {{opacity, .5}, 0.1, 1}, 
 {face, Blue}]


Apologies if this image has copyright implications...


Posted 2012-07-19T22:33:47.067

Reputation: 23 565

19This face is copyrighted. You'll have to replace it with a non-hyperbolic heptagon to avoid being sued. – rm -rf – 2012-07-20T14:12:48.660

1Could you explain the algorithm? – rcollyer – 2012-07-21T20:00:00.477

@rcollyer good thinking, added context – cormullion – 2012-07-22T09:38:55.903

1My, +1, then. I'd like to see if you could get the spiral engraving working, too. – rcollyer – 2012-07-22T12:55:21.407


Sorry for bumping a two-year-old question, but I asked about it and people said it was okay!

Since, as the question states, "true" vectorization of the whole image is a tough problem, how about we break up the image into little blocks and vectorize those?

I'll take one of the Kodak test images:

img = Import[""];

enter image description here

blocks = ImagePartition[img, 32];
fromImg[p_] := {1, -1} Reverse@p (* image coordinates to graphics coordinates *)
drawBlocks[f_, blocks_] := 
 Graphics[MapIndexed[Translate[f[#1], fromImg@#2] &, blocks, {2}], ImageSize -> Large]

For each block, let's do the simplest vectorization possible: make a straight-line slice through the block and draw the two pieces with different colours. Now the question is where and in what direction should we slice, and what should the two colours be? Instead of cooking up some heuristics, we'll let Mathematica find the slice and colours that best match the pixels of the original block. This is a bit involved, though...

lerp[a_, b_, t_] := a + (b - a) t (* linear interpolation *)
clip[poly_, f_] := (* Sutherland-Hodgman polygon clipping with a single clipping edge *)
 MapThread[{If[f[#1] >= 0, #1],
   With[{f1 = f[#1], f2 = f[#2]}, If[f1 f2 < 0, lerp[#1, #2, f1/(f1 - f2)]]]} &,
  {poly, RotateLeft@poly}] // Flatten[#, 1] & // DeleteCases[#, Null] &
slice[block_] := 
 Module[{m, n, params, interp, line, fun, error, sol, rect},
  {m, n} = ImageDimensions[block];
  params = {θ, d, {r1, g1, b1}, {r2, g2, b2}}; (* direction, position, colour1, colour2 *)
  smoothStep[t_] := (ArcTan[t] + π/2)/π;
  line[{x_, y_}, {θ_, d_}] := ({x, y} - {m, n}/2).{Cos@θ, Sin@θ} - d;
  fun[{i_, j_}, {θ_, d_, c1 : {r1_, g1_, b1_}, c2 : {r2_, g2_, b2_}}] := 
   lerp[c1, c2, smoothStep@line[{i - 1/2, j - 1/2}, {θ, d}]];
  error[params_] := Total[MapIndexed[Norm[#1 - fun[#2, params]]^2 &, ImageData@block, {2}], ∞];
  sol = Last@FindMinimum[error[params], Flatten@params]; (* this is where the magic happens *)
  rect = {{0, 0}, {m, 0}, {m, n}, {0, n}};
  {RGBColor[r1, g1, b1] /. sol, 
   Polygon[fromImg[#/{m, n}] & /@ clip[rect, -line[#, {θ, d} /. sol] &]], 
   RGBColor[r2, g2, b2] /. sol, 
   Polygon[fromImg[#/{m, n}] & /@ clip[rect, line[#, {θ, d} /. sol] &]]}]
drawBlocks[slice, blocks] (* takes a while *)

It looks odd in Mathematica because of antialiasing artifacts, but here's a screenshot of an exported PDF. I quite like the chunky sort of look, though the eyes are a little wonky.

enter image description here

Improved version with $16\times16$ blocks, NMinimize instead of FindMinimum, and a very very very long time (click for bigger):

enter image description here

Also, everyone loves GIFs:

enter image description here


Posted 2012-07-19T22:33:47.067



It's not only "okay" it's encouraged. +1 on your way toward that badge. As an extension consider drawing smaller squares (or rectangles?) in more detailed areas of the image.

– Mr.Wizard – 2014-07-02T20:22:53.080

@Mr.Wizard: Thanks! If one chose to draw smaller squares in more detailed areas one could skip the slicing part entirely. :)

– None – 2014-07-02T20:28:24.700

2Personally I'd prefer to see a combination of the two. – Mr.Wizard – 2014-07-02T20:34:49.060


Original version

It was stated in the comments to Vitaliy's answer that it would be tedious to rewrite it so that the colors of the regions correspond to the colors in the original image. For that reason I would like to supply a different way of "vectorizing" with a Voronoi diagram, which is quite straightforward. For a set of sample points pts:

{w, h} = ImageDimensions[img];
{coords, polys} = BoundedDiagram[{{0, 0}, {0, h}, {w, 0}, {w, h}}, pts];
colorList = MapIndexed[First@#2 -> RGBColor @@ Extract[ImageData[img, DataReversed -> True], Reverse@#] &, pts];
Graphics[{First@# /. colorList, GraphicsComplex[coords, Polygon@Last@#]} & /@ polys]

voronoi example

The sample points can be selected in many different ways. The easiest way is to select points at random, and the best way is a Poisson-disc distribution. The latter is actually how our eyes sample light. This distribution selects sample points so that each Voronoi region is approximately the same size. This ensures that the same quality is achieved in every part of the final image. The poisson-disc distribution is a bit involved, but there is an approximation which is pretty good which is known as the best-candidate distribution. I used that to generate the sample points for the image above. If you are interested you can read more here. The code is given below:

findBestCandidate[samplePoints_, nrOfCandidates_, {w_, h_}] := 
  {nearestFunction, candidates},
  nearestFunction = Nearest[samplePoints];
  candidates = Transpose[{
     RandomInteger[{1, w}, nrOfCandidates],
     RandomInteger[{1, h}, nrOfCandidates]
  Last@SortBy[candidates, Norm[# - First@nearestFunction@#] &]

sample[nrOfSamplePoints_, nrOfCandidates_, {w_, h_}] := Nest[
  #~Append~findBestCandidate[#, nrOfCandidates, {w, h}] &,
  {{RandomInteger[w], RandomInteger[h]}},

The more candidates, the better distribution. The more sample points the better the quality of the final image. The vast majority of the computing time, in the end, goes to calculating the bounded diagram and plotting it.

For Mathematica version 10 and above (written by R.M.):

The builtin VoronoiMesh is much faster than BoundedDiagram and the above code can be written using new functions as:

{w, h} = ImageDimensions[img];
pts = sample[2000, 20, {w, h}];
mesh = VoronoiMesh[pts, {{0, w}, {0, h}}];
coords = MeshPrimitives[mesh, 0];
polys = MeshPrimitives[mesh, 2];
colorList = RGBColor @@ Extract[ImageData[img, DataReversed -> True], Reverse@#] & /@ 
    (Mean @@@ List @@@ polys /. x_?NumericQ :> Round[x]);
Graphics[Transpose@{colorList, polys}]

C. E.

Posted 2012-07-19T22:33:47.067

Reputation: 67 448

Using VoronoiMesh in v10 makes this much faster than BoundedDiagram and allows you to sample more points and still finish in under 20-30 secs. (It does require changing some of the structure of colorList and the call to Graphics.) Probably too long to post as a comment, but I can add it as an edit it if you'd like. – rm -rf – 2016-01-26T15:21:41.377

@R.M. I would appreciate it very much if you did. – C. E. – 2016-01-26T16:11:04.050

Please take a look at the edit when you can :) – rm -rf – 2016-01-26T16:30:17.483

@R.M. It works great! – C. E. – 2016-01-26T16:52:52.093


As of version 11.1 we can use ImageGraphics.

The image from Vitaliy's answer:

im = ImageResize[Import[""], 300]

enter image description here

g = ImageGraphics[im, MinColorDistance -> 0]

enter image description here

We can get rid of the boundary artifacts by adding an EdgeForm:

g /. c_?ColorQ :> Directive[EdgeForm[c], c]

enter image description here

We can watch the detail grow by specifying an increasing amount of colors to quantize with:

frames = DeleteDuplicates @ Table[
  ImageGraphics[im, n, MinColorDistance -> 0], 
  {n, 2, 100}
] /. c_?ColorQ :> Directive[EdgeForm[c], c];

enter image description here

Chip Hurst

Posted 2012-07-19T22:33:47.067

Reputation: 29 735