How to use Mathematica to turn a picture into low poly style?

Question

I want to use Mathematica to turn a picture into low poly style. That is, to make the following conversion: Bill to Bill low poly(Using Image Triangulator App).

I think the biggest difficulty is how to extract the points on the boundaries of different colors appropriately. Too much points make the output picture messy, while the picture generated by too little points are not beautiful.

And then, after I have extract the points of the original picture, how to turn the original picture to the low poly style? Is using DelaunayMesh possible? Or, extract the color of a certain point in a triangle(using ImageValue or some functions similar) and make it to be the main color of the triangle?

Edit: No, I don't think that the answers in artistic image vectorization fit my problem. Because you can see that the question is not using Delaunay mesh or other triangular mesh.

Answer 1

My take on this, with some inspiration from shrx's great answer. I pick points randomly, but weighted by the image gradient to try to get more points on edges in the image. Then I use a Delaunay triangulation so that those edges are maintained in the mesh. Finally I colour it by taking the single pixel value at the centre of each triangle (rather than the slower but more accurate averaging approach used by shrx).

i = ExampleData[{"TestImage", "House"}];
n = 1000;

{x, y} = ImageDimensions[i];

pts = Reverse /@ RandomChoice[
    Flatten @ ImageData @ GradientFilter[i, 2] -> Tuples @ {Range[y, 1, -1], Range[x]}, n];

pts = Join[pts, {{0, 0}, {x, 0}, {x, y}, {0, y}}];

m = DelaunayMesh @ pts;

Graphics[With[{col = RGBColor @ ImageValue[i, Mean @@ #]},
    {EdgeForm @ col, col, #}] & /@ MeshPrimitives[m, 2]]

Answer 2

i = ExampleData[{"TestImage", "Mandrill"}];
id = ImageDimensions[i];

Create some keypoints and use them to make a triangular mesh:

xy = ImageKeypoints[i, MaxFeatures -> 50];
m = VoronoiMesh[xy, {0, #} & /@ id];
mt = TriangulateMesh[#, MaxCellMeasure -> ∞, 
     MeshQualityGoal -> "Minimal"] & /@ 
   Map[MeshRegion[MeshCoordinates[m], #] &, MeshCells[m, 2]];
p = Flatten[MeshPrimitives[#, 2] & /@ mt];

Create masks for calculations of mean color of triangles:

bp = Transpose[{Table[Black, {Length[p]}], p}];
masks = ParallelTable[
   Rasterize[
    Graphics[{EdgeForm[Black], ReplacePart[bp, {n, 1} -> White]}, 
     PlotRangePadding -> None], ImageSize -> id], {n, Length[p]}];
c = ParallelMap[
  RGBColor @@ 
    ImageMeasurements[i, "Mean", Masking -> Opening[#, 1]] &, masks];

Put the colored triangles together:

Graphics[Transpose[{c, EdgeForm /@ c, p}], PlotRangePadding -> None]