4 replaced http://mathematica.stackexchange.com/ with https://mathematica.stackexchange.com/

## 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:

Needs["ComputationalGeometry"]
{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]


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_}] :=
Module[
{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]}},
nrOfSamplePoints
]


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.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}]


## 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:

Needs["ComputationalGeometry"]
{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]


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_}] :=
Module[
{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]}},
nrOfSamplePoints
]


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}]


## 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:

Needs["ComputationalGeometry"]
{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]


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_}] :=
Module[
{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]}},
nrOfSamplePoints
]


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}]

3 added 93 characters in body

## 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:

Needs["ComputationalGeometry"]
{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]


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_}] :=
Module[
{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]}},
nrOfSamplePoints
]


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}]


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:

Needs["ComputationalGeometry"]
{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]


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_}] :=
Module[
{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]}},
nrOfSamplePoints
]


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:

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}]


## 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:

Needs["ComputationalGeometry"]
{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]


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_}] :=
Module[
{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]}},
nrOfSamplePoints
]


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}]

2 added 568 characters in body

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:

Needs["ComputationalGeometry"]
{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]


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_}] :=
Module[
{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]}},
nrOfSamplePoints
]


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:

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}]


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:

Needs["ComputationalGeometry"]
{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]


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_}] :=
Module[
{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]}},
nrOfSamplePoints
]


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.

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:

Needs["ComputationalGeometry"]
{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]


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_}] :=
Module[
{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]}},
nrOfSamplePoints
]


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:

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}]
`
1