text:
Draw an Image as a Voronoi Map
Credits to Calvin's Hobbies for nudging my challenge idea in the right direction.
Consider a set of points in the plane, which we will call sites, and associate a colour with each site. Now you can paint the entire plane by colouring each point with the colour of the closest site. This is called a Voronoi map (or Voronoi diagram). In principle, Voronoi maps can be defined for any distance metric, but we'll simply use the usual Euclidean distance r = √(x² + y²). (Note: You do not necessarily have to know how to compute and render one of these to compete in this challenge.)
Here is an example with 100 sites:

If you look at any cell, then all points within that cell are closer to the corresponding site than to any other site.
Your task is to approximate a given image with such a Voronoi map. You're given the image in any convenient raster graphics format, as well as an integer N. You should then produce up to N sites, and a colour for each site, such that the Voronoi map based on these sites resembles the input image as closely as possible.
You can use the Stack Snippet at the bottom of this challenge to render a Voronoi map from your output, or you can render it yourself if you prefer.
You may use built-in or third-party functions to compute a Voronoi map from a set of sites (if you need to).
This is a popularity contest, so the answer with the most net votes wins. Voters are encouraged to judge answers by

how well the original images and their colours are approximated.
how well the algorithm works on different kinds of images.
how well the algorithm works for small N.
whether the algorithm adaptively clusters points in regions of the image that require more detail.

Test Images
Here are a few images to test your algorithm on (some of our usual suspects, some new ones). Click the pictures for larger versions.

The beach in the first row was drawn by Olivia Bell, and included with her permission.
If you want an extra challenge, try Yoshi with a white background and get his belly line right.
You can find all of these test images in this imgur gallery where you can download all of them as a zip file. The album also contains a random Voronoi diagram as another test. For reference, here is the data that generated it.
Please include example diagrams for a variety of different images and N, e.g. 100, 300, 1000, 3000 (as well as pastebins to some of the corresponding cell specifications). You may use or omit black edges between the cells as you see fit (this may look better on some images than on others). Do not include the sites though (except in a separate example maybe if you want to explain how your site placement works, of course).
If you want to show a large number of results, you can create a gallery over on imgur.com, to keep the size of the answers reasonable. Alternatively, put thumbnails in your post and make them links to larger images, like I did in my reference answer. You can get the small thumbnails by appending s to the file name in the imgur.com link (e.g. I3XrT.png -> I3XrTs.png). Also, feel free to use other test images, if you find something nice.
Renderer
Paste your output into the following Stack Snippet to render your results. The exact list format is irrelevant, as long as each cell is specified by 5 floating point numbers in the order x y r g b, where x and y are the coordinates of the cell's site, and r g b are the red, green and blue colour channels in the range 0 ≤ r, g, b ≤ 1.
The snippet provides options to specify a line width of the cell edges, and whether or not the cell sites should be shown (the latter mostly for debugging purposes). But note that the output is only re-rendered when the cell specifications change - so if you modify some of the other options, add a space to the cells or something.

function draw() {
  document.getElementById("output").innerHTML = svg
}

function drawMap() {
  var cells = document.getElementById("cells").value;
  var showSites = document.getElementById("showSites").checked;
  var showCells = document.getElementById("showCells").checked;
  var lineWidth = parseFloat(document.getElementById("linewidth").value);
  var width = parseInt(document.getElementById("width").value);
  var height = parseInt(document.getElementById("height").value);
  var e = prefix.replace('{{WIDTH}}', width)
                .replace('{{HEIGHT}}', height);
  cells = cells.match(/(?:\d*\.\d+|\d+\.\d*|\d+)(?:e-?\d+)?/ig);
  var sitesDom = '';
  var sites = []
  for (i = 0; i < cells.length; i+=5)
  {
    x = parseFloat(cells[i]);
    y = parseFloat(cells[i+1]);
    r = Math.floor(parseFloat(cells[i+2])*255);
    g = Math.floor(parseFloat(cells[i+3])*255);
    b = Math.floor(parseFloat(cells[i+4])*255);
    sitesDom += '<circle cx="' + x + '" + cy="' + y + '" r="1" fill="black"/>';
    sites.push({ x: x, y: y, r: r, g: g, b: b });
  }

  if (showCells)
  {
    var bbox = { xl: 0, xr: width, yt: 0, yb: height };
    var voronoi = new Voronoi();

    var diagram = voronoi.compute(sites, bbox);
    for (i = 0; i < diagram.cells.length; ++i)
    {
      var cell = diagram.cells[i];
      var site = cell.site;
      var coords = '';
      for (var j = 0; j < cell.halfedges.length; ++j)
      {
        var vertex = cell.halfedges[j].getStartpoint();
        coords += ' ' + vertex.x + ',' + vertex.y;
      }
      var color = 'rgb('+site.r+','+site.g+','+site.b+')';
      e += '<polygon fill="'+color+'" points="' + coords + '" stroke="' + (lineWidth ? 'black' : color) + '" stroke-width="'+Math.max(0.6,lineWidth)+'"/>';
    }
  }
  if (showSites)
    e += sitesDom;
  e += suffix;
  e += '
 Using '+sites.length+' cells.';
  svg = e;
  draw();
}
var prefix = '<svg width="{{WIDTH}}" height="{{HEIGHT}}">',
  suffix = "</svg>",
  scale = 0.95,
  offset = 225,
  radius = scale*offset,
  svg = "";
svg {
  position: relative;
}

Line width: <input id="linewidth" type="text" size="1" value="0" />

Show sites: <input id="showSites" type="checkbox" />

Show cells: <input id="showCells" type="checkbox" checked="checked" />

Dimensions: <input id="width" type="text" size="1" value="400" /> x <input id="height" type="text" size="1" value="300" />

Paste cell specifications here

<textarea id="cells" cols="60" rows="10" oninput='drawMap()'></textarea>

Massive credits to Raymond Hill for writing this really nice JS Voronoi library.
Related Challenges

Paint by Numbers
Photomosaics or: How Many Programmers Does it Take to Replace a Light Bulb?

How do we verify we have done it the best way?

@frogeyedpeas By looking at the votes you get. ;) This is a popularity contest. There is not necessarily the best way to do. The idea is that you try to do it as well as you can, and the votes will reflect whether people agree that you've done a good job. There is a certain amount of subjectivity in these, admittedly. Have a look at the related challenges I linked, or at this one. You'll see that there is usually a wide variety of approaches but the voting system helps the better solutions bubble up to the top and decide on a winner.

Olivia approves of the approximations of her beach submitted thus far.

@AlexA. Devon approves of some of the approximations of his face submitted thus far. He's not a big fan of any of the n=100 versions ;)

@Geobits: He'll understand when he's older.

Here's a page about a centroidal voronoi-based stippling technique. A good source of inspiration (the related master thesis has a nice discussion of possible improvements to the algorithm).

Did I make a comment about the use of the term "distance measure" and how it should be "metric" or am I going mental? I can't find the comment!

@AlecTeal You did, and I fixed it, and I replied... and I saw you were online since my reply, so I assumed you'd read it and deleted both comments as they were obsolete. ;)

I had linked to the definition because.... I have a maths background, but it might be the first time someone has seen a formal notion of "distance" and their "metrics" - that's why I included the link. You're using just 3 properties to generate this map. I am also quite interested in seeing what the pictures look like with different metrics.... http://www.maths.kisogo.com/index.php?title=Metric_space is the link BTW

Is this just a puzzle or do you also think this could be an approach for a new type of compressed images? To me it looks like storing the coordinates and colors of the sites would need much less information than storing the image pixel by pixel.

@ThomasWeller I just thought it's a fun challenge that is both challenging and would yield varied and interesting results. But the same thing has been discussed when I posted it in Hacker News, so maybe you can weigh in there.

Python + scipy + scikit-image, weighted Poisson disc sampling
My solution is rather complex. I do some preprocessing on the image to remove noise and get a mapping how 'interesting' each point is (using a combination of local entropy and edge detection):

Then I choose sampling points using Poisson disc sampling with a twist: the distance of the circle is determined by the weight we determined earlier.
Then once I have the sampling points I divide up the image in voronoi segments and assign the L*a*b* average of the color values inside each segment to each segment.
I have a lot of heuristics, and I also must do a bit of math to make sure the number of sample points is close to N. I get N exactly by overshooting slightly, and then dropping some points with an heuristic.
In terms of runtime, this filter isn't cheap, but no image below took more than 5 seconds to make.
Without further ado:
import math
import random
import collections
import os
import sys
import functools
import operator as op
import numpy as np
import warnings

from scipy.spatial import cKDTree as KDTree
from skimage.filters.rank import entropy
from skimage.morphology import disk, dilation
from skimage.util import img_as_ubyte
from skimage.io import imread, imsave
from skimage.color import rgb2gray, rgb2lab, lab2rgb
from skimage.filters import sobel, gaussian_filter
from skimage.restoration import denoise_bilateral
from skimage.transform import downscale_local_mean

# Returns a random real number in half-open range [0, x).
def rand(x):
    r = x
    while r == x:
        r = random.uniform(0, x)
    return r

def poisson_disc(img, n, k=30):
    h, w = img.shape[:2]

    nimg = denoise_bilateral(img, sigma_range=0.15, sigma_spatial=15)
    img_gray = rgb2gray(nimg)
    img_lab = rgb2lab(nimg)

    entropy_weight = 2**(entropy(img_as_ubyte(img_gray), disk(15)))
    entropy_weight /= np.amax(entropy_weight)
    entropy_weight = gaussian_filter(dilation(entropy_weight, disk(15)), 5)

    color = [sobel(img_lab[:, :, channel])**2 for channel in range(1, 3)]
    edge_weight = functools.reduce(op.add, color) ** (1/2) / 75
    edge_weight = dilation(edge_weight, disk(5))

    weight = (0.3*entropy_weight + 0.7*edge_weight)
    weight /= np.mean(weight)
    weight = weight

    max_dist = min(h, w) / 4
    avg_dist = math.sqrt(w * h / (n * math.pi * 0.5) ** (1.05))
    min_dist = avg_dist / 4

    dists = np.clip(avg_dist / weight, min_dist, max_dist)

    def gen_rand_point_around(point):
        radius = random.uniform(dists[point], max_dist)
        angle = rand(2 * math.pi)
        offset = np.array([radius * math.sin(angle), radius * math.cos(angle)])
        return tuple(point + offset)

    def has_neighbours(point):
        point_dist = dists[point]
        distances, idxs = tree.query(point,
                                    len(sample_points) + 1,
                                    distance_upper_bound=max_dist)

        if len(distances) == 0:
            return True

        for dist, idx in zip(distances, idxs):
            if np.isinf(dist):
                break

            if dist < point_dist and dist < dists[tuple(tree.data[idx])]:
                return True

        return False

    # Generate first point randomly.
    first_point = (rand(h), rand(w))
    to_process = [first_point]
    sample_points = [first_point]
    tree = KDTree(sample_points)

    while to_process:
        # Pop a random point.
        point = to_process.pop(random.randrange(len(to_process)))

        for _ in range(k):
            new_point = gen_rand_point_around(point)

            if (0 <= new_point[0] < h and 0 <= new_point[1] < w
                    and not has_neighbours(new_point)):
                to_process.append(new_point)
                sample_points.append(new_point)
                tree = KDTree(sample_points)
                if len(sample_points) % 1000 == 0:
                    print("Generated {} points.".format(len(sample_points)))

    print("Generated {} points.".format(len(sample_points)))

    return sample_points

def sample_colors(img, sample_points, n):
    h, w = img.shape[:2]

    print("Sampling colors...")
    tree = KDTree(np.array(sample_points))
    color_samples = collections.defaultdict(list)
    img_lab = rgb2lab(img)
    xx, yy = np.meshgrid(np.arange(h), np.arange(w))
    pixel_coords = np.c_[xx.ravel(), yy.ravel()]
    nearest = tree.query(pixel_coords)[1]

    i = 0
    for pixel_coord in pixel_coords:
        color_samples[tuple(tree.data[nearest[i]])].append(
            img_lab[tuple(pixel_coord)])
        i += 1

    print("Computing color means...")
    samples = []
    for point, colors in color_samples.items():
        avg_color = np.sum(colors, axis=0) / len(colors)
        samples.append(np.append(point, avg_color))

    if len(samples) > n:
        print("Downsampling {} to {} points...".format(len(samples), n))

    while len(samples) > n:
        tree = KDTree(np.array(samples))
        dists, neighbours = tree.query(np.array(samples), 2)
        dists = dists[:, 1]
        worst_idx = min(range(len(samples)), key=lambda i: dists[i])
        samples[neighbours[worst_idx][1]] += samples[neighbours[worst_idx][0]]
        samples[neighbours[worst_idx][1]] /= 2
        samples.pop(neighbours[worst_idx][0])

    color_samples = []
    for sample in samples:
        color = lab2rgb([[sample[2:]]])[0][0]
        color_samples.append(tuple(sample[:2][::-1]) + tuple(color))

    return color_samples

def render(img, color_samples):
    print("Rendering...")
    h, w = [2*x for x in img.shape[:2]]
    xx, yy = np.meshgrid(np.arange(h), np.arange(w))
    pixel_coords = np.c_[xx.ravel(), yy.ravel()]

    colors = np.empty([h, w, 3])
    coords = []
    for color_sample in color_samples:
        coord = tuple(x*2 for x in color_sample[:2][::-1])
        colors[coord] = color_sample[2:]
        coords.append(coord)

    tree = KDTree(coords)
    idxs = tree.query(pixel_coords)[1]
    data = colors[tuple(tree.data[idxs].astype(int).T)].reshape((w, h, 3))
    data = np.transpose(data, (1, 0, 2))

    return downscale_local_mean(data, (2, 2, 1))

if __name__ == "__main__":
    warnings.simplefilter("ignore")

    img = imread(sys.argv[1])[:, :, :3]

    print("Calibrating...")
    mult = 1.02 * 500 / len(poisson_disc(img, 500))

    for n in (100, 300, 1000, 3000):
        print("Sampling {} for size {}.".format(sys.argv[1], n))

        sample_points = poisson_disc(img, mult * n)
        samples = sample_colors(img, sample_points, n)
        base = os.path.basename(sys.argv[1])
        with open("{}-{}.txt".format(os.path.splitext(base)[0], n), "w") as f:
            for sample in samples:
                f.write(" ".join("{:.3f}".format(x) for x in sample) + "\n")

        imsave("autorenders/{}-{}.png".format(os.path.splitext(base)[0], n),
            render(img, samples))

        print("Done!")

Images
Respectively N is 100, 300, 1000 and 3000:

I like this; it looks a bit like smoked glass.

I messed around with this a bit, and you get better results, particularly for the low triangle images, if you replace the denoise_bilatteral with a denoise_tv_bregman. It generates more even patches in its denoising, which helps.

@LKlevin What weight did you use?

I used 1.0 as the weight.

C++
My approach is quite slow, but I'm very happy with quality of the results that it gives, particularly with respect to preserving edges.  For example, here's Yoshi and the Cornell Box with just 1000 sites each:

There are two main parts that make it tick.  The first, embodied in the evaluate() function takes a set of candidate site locations, sets the optimal colors on them and returns a score for the PSNR of the rendered Voronoi tessellation versus the target image.  The colors for each site are determined by averaging target image pixels covered by the cell around the site.  I use Welford's algorithm to help compute both the best color for each cell and the resulting PSNR using just a single pass over the image by exploiting the relationship between variance, MSE, and PSNR.  This reduces the problem to one of finding the best set of site locations without particular regard to color.
The second part then, embodied in main(), tries to find this set.  It starts by choosing a set of points at random.  Then, in each step it removes one point (going round-robin) and tests a set of random candidate points to replace it.  The one that yields the highest PSNR of the bunch is accepted and kept.  Effectively, this causes the site to jump to a new location and generally improves the image bit-by-bit.  Note that the algorithm intentionally does not retain the original location as a candidate.  Sometimes this means that the jump lowers the overall image quality.  Allowing this to happen helps to avoid getting stuck in local maxima.  It also gives a stopping criteria; the program terminates after going a certain number of steps without improving on the best set of sites found so far.
Note that this implementation is fairly basic and can easily take hours of CPU-core time, especially as the number of sites grows.  It recomputes the complete Voronoi map for every candidate and brute force tests the distance to all sites for each pixel.  Since each operation involves removing one point at a time and adding another, the actual changes to the image at each step are going to be fairly local.  There are algorithms for efficiently incrementally updating a Voronoi diagram and I believe they'd give this algorithm a tremendous speedup.  For this contest, however, I've chosen to keep things simple and brute-force.
Code
#include <cstdlib>
#include <cmath>
#include <string>
#include <vector>
#include <fstream>
#include <istream>
#include <ostream>
#include <iostream>
#include <algorithm>
#include <random>

static auto const decimation = 2;
static auto const candidates = 96;
static auto const termination = 200;

using namespace std;

struct rgb {float red, green, blue;};
struct img {int width, height; vector<rgb> pixels;};
struct site {float x, y; rgb color;};

img read(string const &name) {
    ifstream file{name, ios::in | ios::binary};
    auto result = img{0, 0, {}};
    if (file.get() != 'P' || file.get() != '6')
        return result;
    auto skip = [&](){
        while (file.peek() < '0' || '9' < file.peek())
            if (file.get() == '#')
                while (file.peek() != '\r' && file.peek() != '\n')
                    file.get();
    };
     auto maximum = 0;
     skip(); file >> result.width;
     skip(); file >> result.height;
     skip(); file >> maximum;
     file.get();
     for (auto pixel = 0; pixel < result.width * result.height; ++pixel) {
         auto red = file.get() * 1.0f / maximum;
         auto green = file.get() * 1.0f / maximum;
         auto blue = file.get() * 1.0f / maximum;
         result.pixels.emplace_back(rgb{red, green, blue});
     }
     return result;
 }

 float evaluate(img const &target, vector<site> &sites) {
     auto counts = vector<int>(sites.size());
     auto variance = vector<rgb>(sites.size());
     for (auto &site : sites)
         site.color = rgb{0.0f, 0.0f, 0.0f};
     for (auto y = 0; y < target.height; y += decimation)
         for (auto x = 0; x < target.width; x += decimation) {
             auto best = 0;
             auto closest = 1.0e30f;
             for (auto index = 0; index < sites.size(); ++index) {
                 float distance = ((x - sites[index].x) * (x - sites[index].x) +
                                   (y - sites[index].y) * (y - sites[index].y));
                 if (distance < closest) {
                     best = index;
                     closest = distance;
                 }
             }
             ++counts[best];
             auto &pixel = target.pixels[y * target.width + x];
             auto &color = sites[best].color;
             rgb delta = {pixel.red - color.red,
                          pixel.green - color.green,
                          pixel.blue - color.blue};
             color.red += delta.red / counts[best];
             color.green += delta.green / counts[best];
             color.blue += delta.blue / counts[best];
             variance[best].red += delta.red * (pixel.red - color.red);
             variance[best].green += delta.green * (pixel.green - color.green);
             variance[best].blue += delta.blue * (pixel.blue - color.blue);
         }
     auto error = 0.0f;
     auto count = 0;
     for (auto index = 0; index < sites.size(); ++index) {
         if (!counts[index]) {
             auto x = min(max(static_cast<int>(sites[index].x), 0), target.width - 1);
             auto y = min(max(static_cast<int>(sites[index].y), 0), target.height - 1);
             sites[index].color = target.pixels[y * target.width + x];
         }
         count += counts[index];
         error += variance[index].red + variance[index].green + variance[index].blue;
     }
     return 10.0f * log10f(count * 3 / error);
 }

 void write(string const &name, int const width, int const height, vector<site> const &sites) {
     ofstream file{name, ios::out};
     file << width << " " << height << endl;
     for (auto const &site : sites)
         file << site.x << " " << site.y << " "
              << site.color.red << " "<< site.color.green << " "<< site.color.blue << endl;
 }

 int main(int argc, char **argv) {
     auto rng = mt19937{random_device{}()};
     auto uniform = uniform_real_distribution<float>{0.0f, 1.0f};
     auto target = read(argv[1]);
     auto sites = vector<site>{};
     for (auto point = atoi(argv[2]); point; --point)
         sites.emplace_back(site{
             target.width * uniform(rng),
             target.height * uniform(rng)});
     auto greatest = 0.0f;
     auto remaining = termination;
     for (auto step = 0; remaining; ++step, --remaining) {
         auto best_candidate = sites;
         auto best_psnr = 0.0f;
         #pragma omp parallel for
         for (auto candidate = 0; candidate < candidates; ++candidate) {
             auto trial = sites;
             #pragma omp critical
             {
                 trial[step % sites.size()].x = target.width * (uniform(rng) * 1.2f - 0.1f);
                 trial[step % sites.size()].y = target.height * (uniform(rng) * 1.2f - 0.1f);
             }
             auto psnr = evaluate(target, trial);
             #pragma omp critical
             if (psnr > best_psnr) {
                 best_candidate = trial;
                 best_psnr = psnr;
             }
         }
         sites = best_candidate;
         if (best_psnr > greatest) {
             greatest = best_psnr;
             remaining = termination;
             write(argv[3], target.width, target.height, sites);
         }
         cout << "Step " << step << "/" << remaining
              << ", PSNR = " << best_psnr << endl;
     }
     return 0;
 }

Running
The program is self-contained and has no external dependencies beyond the standard library, but it does require images to be in the binary PPM format.  I use ImageMagick to convert images to PPM, though GIMP and quite a few other programs can do it too.
To compile it, save the program as voronoi.cpp and then run:
g++ -std=c++11 -fopenmp -O3 -o voronoi voronoi.cpp

I expect it will probably work on Windows with recent versions of Visual Studio, though I haven't tried this.  You'll want to make sure that you're compiling with C++11 or better and OpenMP enabled if you do.  OpenMP is not strictly necessary, but it helps a lot in making the execution times more tolerable.
To run it, do something like:
./voronoi cornell.ppm 1000 cornell-1000.txt

The later file will be updated with the site data as it goes.  The first line will have the width and height of the image, followed by lines of x, y, r, g, b values suitable for copying and pasting into the Javascript renderer in the problem description.
The three constants at the top of the program allow you to tune it for speed versus quality.  The decimation factor coarsens the target image when evaluating a set of sites for color and PSNR.  The higher it is, the faster the program will run.  Setting it to 1 uses the full resolution image.  The candidates constant controls how many candidates to test on each step.  Higher gives a better chance of finding a good spot to jump to but makes the program slower.  Finally, termination is how many steps the program can go without improving its output before it quits.  Increasing it may give better results but make it take marginally longer.
Images
N = 100, 300, 1000, and 3000:

This should've won IMO - much better than mine.

@orlp -- Thanks!  To be fair though, you posted yours much sooner and it runs much more quickly.  Speed counts!

Well, mine isn't really a voronoi map answer :) It's a really really good sampling algorithm, but turning sample points into voronoi sites is clearly not optimal.

IDL, adaptive refinement
This method is inspired by Adaptive Mesh Refinement from astronomical simulations, and also Subdivision surface. This is the kind of task that IDL prides itself on, which you'll be able to tell by the large number of builtin functions I was able to use. :D
I've output some of the intermediates for the black-background yoshi test image, with n = 1000.
First, we perform a luminance greyscale on the image (using ct_luminance), and apply a Prewitt filter (prewitt, see wikipedia) for good edge detection:

Then comes the real grunt-work: we subdivide the image into 4, and measure the variance in each quadrant in the filtered image. Our variance is weighted by the size of the subdivision (which in this first step is equal), so that "edgy" regions with high variance don't get subdivided smaller and smaller and smaller. Then, we use the weighted variance to target subdivisions with more detail, and iteratively subdivide each detail-rich section into 4 additional, until we hit our target number of sites (each subdivision contains exactly one site). Since we're adding 3 sites each time we iterate, we end up with n - 2 <= N <= n sites.
I made a .webm of the subdivision process for this image, which I can't embed, but it's here; the color in each subsection is determined by the weighted variance. (I made the same kind of video for the white-background yoshi, for comparison, with the color table reversed so it goes toward white instead of black; it's here.) The final product of the subdivision looks like this:

Once we have our list of subdivisions, we go through each subdivision. The final site location is the position of the minimum of the Prewitt image, i.e., the least "edgy" pixel, and the color of the section is the color of that pixel; here's the original image, with sites marked:

Then, we use the built-in triangulate to calculate the Delaunay triangulation of the sites, and the built-in voronoi to define the vertices of each Voronoi polygon, before drawing each polygon to an image buffer in its respective color. Finally, we save a snapshot of the image buffer.

The code:
function subdivide, image, bounds, vars
  ;subdivide a section into 4, and return the 4 subdivisions and the variance of each
  division = list()
  vars = list()
  nx = bounds[2] - bounds[0]
  ny = bounds[3] - bounds[1]
  for i=0,1 do begin
    for j=0,1 do begin
      x = i * nx/2 + bounds[0]
      y = j * ny/2 + bounds[1]
      sub = image[x:x+nx/2-(~(nx mod 2)),y:y+ny/2-(~(ny mod 2))]
      division.add, [x,y,x+nx/2-(~(nx mod 2)),y+ny/2-(~(ny mod 2))]
      vars.add, variance(sub) * n_elements(sub)
    endfor
  endfor
  return, division
end

pro voro_map, n, image, outfile
  sz = size(image, /dim)
  ;first, convert image to greyscale, and then use a Prewitt filter to pick out edges
  edges = prewitt(reform(ct_luminance(image[0,*,*], image[1,*,*], image[2,*,*])))
  ;next, iteratively subdivide the image into sections, using variance to pick
  ;the next subdivision target (variance -> detail) until we've hit N subdivisions
  subdivisions = subdivide(edges, [0,0,sz[1],sz[2]], variances)
  while subdivisions.count() lt (n - 2) do begin
    !null = max(variances.toarray(),target)
    oldsub = subdivisions.remove(target)
    newsub = subdivide(edges, oldsub, vars)
    if subdivisions.count(newsub[0]) gt 0 or subdivisions.count(newsub[1]) gt 0 or subdivisions.count(newsub[2]) gt 0 or subdivisions.count(newsub[3]) gt 0 then stop
    subdivisions += newsub
    variances.remove, target
    variances += vars
  endwhile
  ;now we find the minimum edge value of each subdivision (we want to pick representative
  ;colors, not edge colors) and use that as the site (with associated color)
  sites = fltarr(2,n)
  colors = lonarr(n)
  foreach sub, subdivisions, i do begin
    slice = edges[sub[0]:sub[2],sub[1]:sub[3]]
    !null = min(slice,target)
    sxy = array_indices(slice, target) + sub[0:1]
    sites[*,i] = sxy
    colors[i] = cgcolor24(image[0:2,sxy[0],sxy[1]])
  endforeach
  ;finally, generate the voronoi map
  old = !d.NAME
  set_plot, 'Z'
  device, set_resolution=sz[1:2], decomposed=1, set_pixel_depth=24
  triangulate, sites[0,*], sites[1,*], tr, connectivity=C
  for i=0,n-1 do begin
    if C[i] eq C[i+1] then continue
    voronoi, sites[0,*], sites[1,*], i, C, xp, yp
    cgpolygon, xp, yp, color=colors[i], /fill, /device
  endfor
  !null = cgsnapshot(file=outfile, /nodialog)
  set_plot, old
end

pro wrapper
  cd, '~/voronoi'
  fs = file_search()
  foreach f,fs do begin
    base = strsplit(f,'.',/extract)
    if base[1] eq 'png' then im = read_png(f) else read_jpeg, f, im
    voro_map,100, im, base[0]+'100.png'
    voro_map,500, im, base[0]+'500.png'
    voro_map,1000,im, base[0]+'1000.png'
  endforeach
end

Call this via voro_map, n, image, output_filename. I included a wrapper procedure as well, which went through each test image and ran with 100, 500, and 1000 sites.
Output collected here, and here are some thumbnails:
n = 100

n = 500

n = 1000

I really like the fact that this solution puts more points in more complex areas, which is I think the intent, and sets it apart from the others at this point.

yeah, the idea of detail-grouped points is what led me to adaptive refinement

Very neat explanation, and the images are impressive! I have a question--It looks like you get much different images when Yoshi is on a white background, where we have a few odd shapes. What might be causing that?

@BrianSteel I imagine the outlines get picked up as high-variance areas and focused on unnecessarily, and then other truly high-detail areas have less points assigned because of that.

@BrainSteel i think doppel is right - there's a strong edge between the black border and the white background, which asks for a lot of detail in the algorithm. i'm not sure if this is something i can (or, more importantly, should) fix...

Python 3 + PIL + SciPy, Fuzzy k-means
from collections import defaultdict
import itertools
import random
import time

from PIL import Image
import numpy as np
from scipy.spatial import KDTree, Delaunay

INFILE = "planet.jpg"
OUTFILE = "voronoi.txt"
N = 3000

DEBUG = True # Outputs extra images to see what's happening
FEATURE_FILE = "features.png"
SAMPLE_FILE = "samples.png"
SAMPLE_POINTS = 20000
ITERATIONS = 10
CLOSE_COLOR_THRESHOLD = 15

"""
Color conversion functions
"""

start_time = time.time()

# http://www.easyrgb.com/?X=MATH
def rgb2xyz(rgb):
  r, g, b = rgb
  r /= 255
  g /= 255
  b /= 255

  r = ((r + 0.055)/1.055)**2.4 if r > 0.04045 else r/12.92
  g = ((g + 0.055)/1.055)**2.4 if g > 0.04045 else g/12.92
  b = ((b + 0.055)/1.055)**2.4 if b > 0.04045 else b/12.92

  r *= 100
  g *= 100
  b *= 100

  x = r*0.4124 + g*0.3576 + b*0.1805
  y = r*0.2126 + g*0.7152 + b*0.0722
  z = r*0.0193 + g*0.1192 + b*0.9505

  return (x, y, z)

def xyz2lab(xyz):
  x, y, z = xyz
  x /= 95.047
  y /= 100
  z /= 108.883

  x = x**(1/3) if x > 0.008856 else 7.787*x + 16/116
  y = y**(1/3) if y > 0.008856 else 7.787*y + 16/116
  z = z**(1/3) if z > 0.008856 else 7.787*z + 16/116

  L = 116*y - 16
  a = 500*(x - y)
  b = 200*(y - z)

  return (L, a, b)

def rgb2lab(rgb):
  return xyz2lab(rgb2xyz(rgb))

def lab2xyz(lab):
  L, a, b = lab
  y = (L + 16)/116
  x = a/500 + y
  z = y - b/200

  y = y**3 if y**3 > 0.008856 else (y - 16/116)/7.787
  x = x**3 if x**3 > 0.008856 else (x - 16/116)/7.787
  z = z**3 if z**3 > 0.008856 else (z - 16/116)/7.787

  x *= 95.047
  y *= 100
  z *= 108.883

  return (x, y, z)

def xyz2rgb(xyz):
  x, y, z = xyz
  x /= 100
  y /= 100
  z /= 100

  r = x* 3.2406 + y*-1.5372 + z*-0.4986
  g = x*-0.9689 + y* 1.8758 + z* 0.0415
  b = x* 0.0557 + y*-0.2040 + z* 1.0570

  r = 1.055 * (r**(1/2.4)) - 0.055 if r > 0.0031308 else 12.92*r
  g = 1.055 * (g**(1/2.4)) - 0.055 if g > 0.0031308 else 12.92*g
  b = 1.055 * (b**(1/2.4)) - 0.055 if b > 0.0031308 else 12.92*b

  r *= 255
  g *= 255
  b *= 255

  return (r, g, b)

def lab2rgb(lab):
  return xyz2rgb(lab2xyz(lab))

"""
Step 1: Read image and convert to CIELAB
"""

im = Image.open(INFILE)
im = im.convert("RGB")
width, height = prev_size = im.size

pixlab_map = {}

for x in range(width):
    for y in range(height):
        pixlab_map[(x, y)] = rgb2lab(im.getpixel((x, y)))

print("Step 1: Image read and converted")

"""
Step 2: Get feature points
"""

def euclidean(point1, point2):
    return sum((x-y)**2 for x,y in zip(point1, point2))**.5

def neighbours(pixel):
    x, y = pixel
    results = []

    for dx, dy in itertools.product([-1, 0, 1], repeat=2):
        neighbour = (pixel[0] + dx, pixel[1] + dy)

        if (neighbour != pixel and 0 <= neighbour[0] < width
            and 0 <= neighbour[1] < height):
            results.append(neighbour)

    return results

def mse(colors, base):
    return sum(euclidean(x, base)**2 for x in colors)/len(colors)

features = []

for x in range(width):
    for y in range(height):
        pixel = (x, y)
        col = pixlab_map[pixel]
        features.append((mse([pixlab_map[n] for n in neighbours(pixel)], col),
                         random.random(),
                         pixel))

features.sort()
features_copy = [x[2] for x in features]

if DEBUG:
    test_im = Image.new("RGB", im.size)

    for i in range(len(features)):
        pixel = features[i][1]
        test_im.putpixel(pixel, (int(255*i/len(features)),)*3)

    test_im.save(FEATURE_FILE)

print("Step 2a: Edge detection-ish complete")

def random_index(list_):
    r = random.expovariate(2)

    while r > 1:
         r = random.expovariate(2)

    return int((1 - r) * len(list_))

sample_points = set()

while features and len(sample_points) < SAMPLE_POINTS:
    index = random_index(features)
    point = features[index][2]
    sample_points.add(point)
    del features[index]

print("Step 2b: {} feature samples generated".format(len(sample_points)))

if DEBUG:
    test_im = Image.new("RGB", im.size)

    for pixel in sample_points:
        test_im.putpixel(pixel, (255, 255, 255))

    test_im.save(SAMPLE_FILE)

"""
Step 3: Fuzzy k-means
"""

def euclidean(point1, point2):
    return sum((x-y)**2 for x,y in zip(point1, point2))**.5

def get_centroid(points):
    return tuple(sum(coord)/len(points) for coord in zip(*points))

def mean_cell_color(cell):
    return get_centroid([pixlab_map[pixel] for pixel in cell])

def median_cell_color(cell):
    # Pick start point out of mean and up to 10 pixels in cell
    mean_col = get_centroid([pixlab_map[pixel] for pixel in cell])
    start_choices = [pixlab_map[pixel] for pixel in cell]

    if len(start_choices) > 10:
        start_choices = random.sample(start_choices, 10)

    start_choices.append(mean_col)

    best_dist = None
    col = None

    for c in start_choices:
        dist = sum(euclidean(c, pixlab_map[pixel])
                       for pixel in cell)

        if col is None or dist < best_dist:
            col = c
            best_dist = dist

    # Approximate median by hill climbing
    last = None

    while last is None or euclidean(col, last) < 1e-6:
        last = col

        best_dist = None
        best_col = None

        for deviation in itertools.product([-1, 0, 1], repeat=3):
            new_col = tuple(x+y for x,y in zip(col, deviation))
            dist = sum(euclidean(new_col, pixlab_map[pixel])
                       for pixel in cell)

            if best_dist is None or dist < best_dist:
                best_col = new_col

        col = best_col

    return col

def random_point():
    index = random_index(features_copy)
    point = features_copy[index]

    dx = random.random() * 10 - 5
    dy = random.random() * 10 - 5

    return (point[0] + dx, point[1] + dy)

centroids = np.asarray([random_point() for _ in range(N)])
variance = {i:float("inf") for i in range(N)}
cluster_colors = {i:(0, 0, 0) for i in range(N)}

# Initial iteration
tree = KDTree(centroids)
clusters = defaultdict(set)

for point in sample_points:
    nearest = tree.query(point)[1]
    clusters[nearest].add(point)

# Cluster!
for iter_num in range(ITERATIONS):
    if DEBUG:
        test_im = Image.new("RGB", im.size)

        for n, pixels in clusters.items():
            color = 0xFFFFFF * (n/N)
            color = (int(color//256//256%256), int(color//256%256), int(color%256))

            for p in pixels:
                test_im.putpixel(p, color)

        test_im.save(SAMPLE_FILE)

    for cluster_num in clusters:
        if clusters[cluster_num]:
            cols = [pixlab_map[x] for x in clusters[cluster_num]]

            cluster_colors[cluster_num] = mean_cell_color(clusters[cluster_num])
            variance[cluster_num] = mse(cols, cluster_colors[cluster_num])

        else:
            cluster_colors[cluster_num] = (0, 0, 0)
            variance[cluster_num] = float("inf")

    print("Clustering (iteration {})".format(iter_num))

    # Remove useless/high variance
    if iter_num < ITERATIONS - 1:
        delaunay = Delaunay(np.asarray(centroids))
        neighbours = defaultdict(set)

        for simplex in delaunay.simplices:
            n1, n2, n3 = simplex

            neighbours[n1] |= {n2, n3}
            neighbours[n2] |= {n1, n3}
            neighbours[n3] |= {n1, n2}

        for num, centroid in enumerate(centroids):
            col = cluster_colors[num]

            like_neighbours = True

            nns = set() # neighbours + neighbours of neighbours

            for n in neighbours[num]:
                nns |= {n} | neighbours[n] - {num}

            nn_far = sum(euclidean(col, cluster_colors[nn]) > CLOSE_COLOR_THRESHOLD
                         for nn in nns)

            if nns and nn_far / len(nns) < 1/5:
                sample_points -= clusters[num]

                for _ in clusters[num]:
                    if features and len(sample_points) < SAMPLE_POINTS:
                        index = random_index(features)
                        point = features[index][3]
                        sample_points.add(point)
                        del features[index]

                clusters[num] = set()

    new_centroids = []

    for i in range(N):
        if clusters[i]:
            new_centroids.append(get_centroid(clusters[i]))
        else:
            new_centroids.append(random_point())

    centroids = np.asarray(new_centroids)
    tree = KDTree(centroids)

    clusters = defaultdict(set)

    for point in sample_points:
        nearest = tree.query(point, k=6)[1]
        col = pixlab_map[point]

        for n in nearest:
            if n < N and euclidean(col, cluster_colors[n])**2 <= variance[n]:
                clusters[n].add(point)
                break

        else:
            clusters[nearest[0]].add(point)

print("Step 3: Fuzzy k-means complete")

"""
Step 4: Output
"""

for i in range(N):
    if clusters[i]:
        centroids[i] = get_centroid(clusters[i])

centroids = np.asarray(centroids)
tree = KDTree(centroids)
color_clusters = defaultdict(set)

# Throw back on some sample points to get the colors right
all_points = [(x, y) for x in range(width) for y in range(height)]

for pixel in random.sample(all_points, int(min(width*height, 5 * SAMPLE_POINTS))):
    nearest = tree.query(pixel)[1]
    color_clusters[nearest].add(pixel)

with open(OUTFILE, "w") as outfile:
    for i in range(N):
        if clusters[i]:
            centroid = tuple(centroids[i])
            col = tuple(x/255 for x in lab2rgb(median_cell_color(color_clusters[i] or clusters[i])))
            print(" ".join(map(str, centroid + col)), file=outfile)

print("Done! Time taken:", time.time() - start_time)

The algorithm
The core idea is that k-means clustering naturally partitions the image into Voronoi cells, since points are tied to the nearest centroid. However, we need to somehow add in the colours as a constraint.
First we convert each pixel to the Lab color space, for better colour manipulation.
Then we do a sort of "poor man's edge detection". For each pixel, we look at its orthogonal and diagonal neighbours, and calculate the mean-squared difference in colour. We then sort all of the pixels by this difference, with pixels most similar to their neighbours at the front of the list, and pixels most dissimilar to their neighbours at the back (i.e., more likely to be an edge point). Here's an example for the planet, where the brighter the pixel is, the more different it is from its neighbours:

(There's a clear grid-like pattern on the rendered output above. According to @randomra, this is probably due to lossy JPG encoding, or imgur compressing the images.)
Next, we use this pixel ordering to sample a large number of points to be clustered. We use an exponential distribution, giving priority to points which are more edge-like and "interesting".

For the clustering, we first pick N centroids, randomly chosen using the same exponential distribution as above. An initial iteration is performed, and for each of the resulting clusters we assign a mean colour and a colour variance threshold. Then for a number of iterations, we:

Build the Delaunay triangulation of the centroids, so that we can easily query neighbours to centroids.
Use the triangulation to remove centroids which are close in colour to most (> 4/5) of their neighbours and neighbour's neighbours combined. Any associated sample points are also removed, and new replacement centroids and sample points are added. This step tries to force the algorithm to place more clusters where detail is needed.
Construct a kd-tree for the new centroids, so that we can easily query the closest centroids to any sample point.
Use the tree to assign each sample point to one of the 6 closest centroids (6 chosen empirically). A centroid will only accept a sample point if the point's colour is within the centroid's colour variance threshold. We try to assign each sample point to the first accepting centroid, but if that's not possible then we simply assign it to the closest centroid. The "fuzziness" of the algorithm comes from this step, since it's possible for clusters to overlap.
Recompute the centroids.

(Click for full size)
Finally, we sample a large number of points, this time using a uniform distribution. Using another kd-tree, we assign each point to its closest centroid, forming clusters. We then approximate the median colour of each cluster using a hill-climbing algorithm, giving our final cell colours (idea for this step thanks to @PhiNotPi and @MartinBüttner).

Notes
In addition to outputting a text file for the snippet (OUTFILE), if DEBUG is set to True the program will also output and overwrite the images above. The algorithm takes a couple of minutes for each image, so it's a good way of checking on progress which doesn't add an awful lot to the running time.
Sample outputs
N = 32:

N = 100:

N = 1000:

N = 3000:

I really like how well your on-white Yoshis turned out.

I love the intermediate images as much as the final result

Mathematica, Random Cells
This is the baseline solution, so you get an idea of the minimum I'm asking from you. Given the file name (locally or as a URL) in file and N in n, the following code simply picks out N random pixels, and uses the colours found at those pixels. This is really naive and doesn't work incredibly well, but I want you guys to beat this after all. :)
data = ImageData@Import@file;
dims = Dimensions[data][[1 ;; 2]]
{Reverse@#, data[[##]][[1 ;; 3]] & @@ Floor[1 + #]} &[dims #] & /@
 RandomReal[1, {n, 2}]

Here are all the test images for N = 100 (all images link to larger versions):

As you can see, these are essentially useless. While they may have some artistic value, in an expressionist way, the original images are barely recognisable.
For N = 500, the situation is improved a bit, but there are still very odd artefacts, the images look washed out, and a lot of cells are wasted on regions without detail:

Your turn!

I'm not a good coder but my god these images are beautiful looking. Awesome idea!

Any reason for Dimensions@ImageData and not ImageDimensions?  You can avoid the slow ImageData altogether by using PixelValue.

@2012rcampion No reason, I just didn't know either function existed. I might fix that up later and also change the example images to the recommended N-values.

Mathematica
We all know Martin loves Mathematica so let me give it a try with Mathematica.
My algorithm uses random points from the image edges to create an initial voronoi diagram. The diagram is then prettified by an iterative adjustment of the mesh with a simple mean filter. This gives images with high cell density near high contrast regions and visually pleasing cells without crazy angles.
The following images show an example of the process in action. The fun is somewhat spoiled by Mathematicas bad Antialiasing, but we get vector graphics, that must be worth something.
This algorithm, without the random sampling, can be found in the VoronoiMesh documentation here.

Test Images (100,300,1000,3000)

Code
VoronoiImage[img_, nSeeds_, iterations_] := Module[{
    i = img,
    edges = EdgeDetect@img,
    voronoiRegion = Transpose[{{0, 0}, ImageDimensions[img]}],
    seeds, voronoiInitial, voronoiRelaxed
    },
   seeds = RandomChoice[ImageValuePositions[edges, White], nSeeds];
   voronoiInitial = VoronoiMesh[seeds, voronoiRegion];
   voronoiRelaxed =
    Nest[VoronoiMesh[Mean @@@ MeshPrimitives[#, 2], voronoiRegion] &,
     voronoiInitial, iterations];
   Graphics[Table[{RGBColor[ImageValue[img, Mean @@ mp]], mp},
     {mp,MeshPrimitives[voronoiRelaxed, 2]}]]
   ];

Nice job for a first post! :) You might want to try the Voronoi test image with 32 cells (that's the number of cells in the image itself).

Thanks! I guess my algorithm will perform terribly on this example. The seeds will get initialized on the cell edges and the recursion won't make it much better ;)

Despite slower rate to converge to the original image, I find that your algorithm yield a very artistry result! (like an improve version of Georges Seurat works). Great job!

You can also get glassy looking interpolated polygon colours by changing your final lines to Graphics@Table[   Append[mp,     VertexColors -> RGBColor /@ ImageValue[img, First[mp]]], {mp,     MeshPrimitives[voronoiRelaxed, 2]}]

Python + SciPy + emcee
The algorithm I've used is the following:

Resize images to a smallish size (~150 pixels)
Make an unsharp-masked image of the maximum channel values (this helps not pick up white regions too strongly).
Take the absolute value.
Choose random points with a probability proportional to this image. This chooses points either side of discontinuities.
Refine the chosen points to lower a cost function. The function is the maximum of the sum of the squared deviations in the channels (again helping bias to solid colours and not only solid white). I've misused Markov Chain Monte Carlo with the emcee module (highly recommended) as an optimiser. The procedure bails out when no new improvement is found after N chain iterations.

The algorithm appears to work very well. Unfortunately it can only sensibly run on smallish images. I haven't had time to take the Voronoi points and apply them to the larger images. They could be refined at this point. I could have also run the MCMC for longer to get better minima. The algorithm's weak point is that it is rather expensive. I haven't had time to increase beyond 1000 points and a couple of the 1000 point images are actually still being refined.
(right click and view image to get a bigger version)

Links to bigger versions are https://i.stack.imgur.com/waU38.jpg (100 points), https://i.stack.imgur.com/IGTAm.jpg (300 points) and https://i.stack.imgur.com/jSIl8.jpg (1000 points)
#!/usr/bin/env python

import glob
import os

import scipy.misc
import scipy.spatial
import scipy.signal
import numpy as N
import numpy.random as NR
import emcee

def compute_image(pars, rimg, gimg, bimg):
    npts = len(pars) // 2
    x = pars[:npts]
    y = pars[npts:npts*2]
    yw, xw = rimg.shape

    # exit if points are too far away from image, to stop MCMC
    # wandering off
    if(N.any(x > 1.2*xw) or N.any(x < -0.2*xw) or
       N.any(y > 1.2*yw) or N.any(y < -0.2*yw)):
        return None

    # compute tesselation
    xy = N.column_stack( (x, y) )
    tree = scipy.spatial.cKDTree(xy)

    ypts, xpts = N.indices((yw, xw))
    queryxy = N.column_stack((N.ravel(xpts), N.ravel(ypts)))

    dist, idx = tree.query(queryxy)

    idx = idx.reshape(yw, xw)
    ridx = N.ravel(idx)

    # tesselate image
    div = 1./N.clip(N.bincount(ridx), 1, 1e99)
    rav = N.bincount(ridx, weights=N.ravel(rimg)) * div
    gav = N.bincount(ridx, weights=N.ravel(gimg)) * div
    bav = N.bincount(ridx, weights=N.ravel(bimg)) * div

    rout = rav[idx]
    gout = gav[idx]
    bout = bav[idx]
    return rout, gout, bout

def compute_fit(pars, img_r, img_g, img_b):
    """Return fit statistic for parameters."""
    # get model
    retn = compute_image(pars, img_r, img_g, img_b)
    if retn is None:
        return -1e99
    model_r, model_g, model_b = retn

    # maximum squared deviation from one of the chanels
    fit = max( ((img_r-model_r)**2).sum(),
               ((img_g-model_g)**2).sum(),
               ((img_b-model_b)**2).sum() )

    # return fake log probability
    return -fit

def convgauss(img, sigma):
    """Convolve image with a Gaussian."""
    size = 3*sigma
    kern = N.fromfunction(
        lambda y, x: N.exp( -((x-size/2)**2+(y-size/2)**2)/2./sigma ),
        (size, size))
    kern /= kern.sum()
    out = scipy.signal.convolve2d(img.astype(N.float64), kern, mode='same')
    return out

def process_image(infilename, outroot, npts):
    img = scipy.misc.imread(infilename)
    img_r = img[:,:,0]
    img_g = img[:,:,1]
    img_b = img[:,:,2]

    # scale down size
    maxdim = max(img_r.shape)
    scale = int(maxdim / 150)
    img_r = img_r[::scale, ::scale]
    img_g = img_g[::scale, ::scale]
    img_b = img_b[::scale, ::scale]

    # make unsharp-masked image of input
    img_tot = N.max((img_r, img_g, img_b), axis=0)
    img1 = convgauss(img_tot, 2)
    img2 = convgauss(img_tot, 32)
    diff = N.abs(img1 - img2)
    diff = diff/diff.max()
    diffi = (diff*255).astype(N.int)
    scipy.misc.imsave(outroot + '_unsharp.png', diffi)

    # create random points with a probability distribution given by
    # the unsharp-masked image
    yw, xw = img_r.shape
    xpars = []
    ypars = []
    while len(xpars) < npts:
        ypar = NR.randint(int(yw*0.02),int(yw*0.98))
        xpar = NR.randint(int(xw*0.02),int(xw*0.98))
        if diff[ypar, xpar] > NR.rand():
            xpars.append(xpar)
            ypars.append(ypar)

    # initial parameters to model
    allpar = N.concatenate( (xpars, ypars) )

    # set up MCMC sampler with parameters close to each other
    nwalkers = npts*5  # needs to be at least 2*number of parameters+2
    pos0 = []
    for i in xrange(nwalkers):
        pos0.append(NR.normal(0,1,allpar.shape)+allpar)

    sampler = emcee.EnsembleSampler(
        nwalkers, len(allpar), compute_fit,
        args=[img_r, img_g, img_b],
        threads=4)

    # sample until we don't find a better fit
    lastmax = -N.inf
    ct = 0
    ct_nobetter = 0
    for result in sampler.sample(pos0, iterations=10000, storechain=False):
        print ct
        pos, lnprob = result[:2]
        maxidx = N.argmax(lnprob)

        if lnprob[maxidx] > lastmax:
            # write image
            lastmax = lnprob[maxidx]
            mimg = compute_image(pos[maxidx], img_r, img_g, img_b)
            out = N.dstack(mimg).astype(N.int32)
            out = N.clip(out, 0, 255)
            scipy.misc.imsave(outroot + '_binned.png', out)

            # save parameters
            N.savetxt(outroot + '_param.dat', scale*pos[maxidx])

            ct_nobetter = 0
            print(lastmax)

        ct += 1
        ct_nobetter += 1
        if ct_nobetter == 60:
            break

def main():
    for npts in 100, 300, 1000:
        for infile in sorted(glob.glob(os.path.join('images', '*'))):
            print infile
            outroot = '%s/%s_%i' % (
                'outdir',
                os.path.splitext(os.path.basename(infile))[0], npts)

            # race condition!
            lock = outroot + '.lock'
            if os.path.exists(lock):
                continue
            with open(lock, 'w') as f:
                pass

            process_image(infile, outroot, npts)

if __name__ == '__main__':
    main()

Unsharped-masked images look like the following. Random points are selected from the image if a random number is less than the value of the image (normed to 1):

I'll post bigger images and the Voronoi points if I get more time.
Edit: If you increase the number of walkers to 100*npts, change the cost function to be some of the squares of the deviations in all channels, and wait for a long time (increasing the number of iterations to break out of the loop to 200), it's possible to make some good images with just 100 points:

Using image energy as a point-weight map
In my approach to this challenge, I wanted a way to map "relevance" of a particular image area to the probability that a particular point would be chosen as a Voronoi centroid. However, I still wanted to preserve the artistic feel of Voronoi mosaicing by randomly choosing image points. Additionally, I wanted to operate on large images, so I don't lose anything in the downsampling process. My algorithm is roughly like this:

For each image, create a sharpness map. The sharpness map is defined by the normalized image energy (or the square of the high frequency signal of the image). An example looks like this:

Generate a number of points from the image, taking 70 percent from the points in the sharpness map and 30 percent from all the other points. This means that points are sampled more densely from high-detail parts of the image.
Color!

Results
N = 100, 500, 1000, 3000

Would you mind a) including the source code used to generate this, and b) link each thumbnail to the full-size image?


text:
Determine whether rational coordinates are in the right Sierpinski triangle
The Sierpinski triangle is a set of points on the plane which is constructed by starting with a single triangle and repeatedly splitting all triangles into four congruent triangles and removing the centre triangle.  The right Sierpinski triangle has corners at (0,0), (0,1) and (1,0), and looks like this:

Some equivalent definitions of this set are as follows:

Points in the nth iteration of the process described above, for all n.

Points (x,y) with 0 <= x <= 1 and 0 <= y <= 1 such that for all positive integers n, the nth bit in the binary expansion of x and y are not both 1.

Let T = {(0,0),(1,0),(0,1)}

Let f be a function on sets of 2D points defined by the following:
f(X) = {(0,0)} ∪ {(x+t)/2 | x∈X, t∈T}
Then the right Sierpinski triangle is the topological closure of the least fixed point (by set containment) of f.

Let S be the square {(x,y) | 0<=x<=1 and 0<=y<=1}

Let g(X) = S ∩ {(x+t)/2 | x∈(X), t∈T} (where T is as defined above)
Then the right Sierpinski triangle is the greatest fixed point of g.
Challenge
Write a program or function which accepts 4 integers, a,b,c,d and gives a truthy value if (a/b,c/d) belongs to the right Sierpinski triangle, and otherwise gives a falsey value.
Scoring
This is a code golf.  Shortest code in bytes wins.
Test cases
The following are in the right Sierpinski triangle:
0 1 0 1
0 1 12345 123456
27 100 73 100
1 7 2 7
8 9 2 21
8 15 20 63
-1 -7 2 7

The following are not in the right Sierpinski triangle:
1 1 1 1
-1 100 1 3
1 3 1 3
1 23 1 7
4 63 3 66
58 217 4351 7577
-1 -7 3 7

Is -1 -3 1 1 a valid input?

Yes, that is a valid input.  I've added test cases to make this clear.

Python 2, 68
lambda n,d,N,D:1>=n/d>=0<=N/D<=1and(n<<abs(D*d))/d&(N<<abs(D*d))/D<1

A nice way to check for gasket membership made ugly. If we were guaranteed that the inputs are non-negative and in the unit square, we would have 38:
lambda n,d,N,D:(n<<D*d)/d&(N<<D*d)/D<1

The idea is that we check whether a point lies within the gasket by checking whether their binary fraction expansions bitwise-AND to 0. To get the first k character of the expansion, we bit-shift the numerator k bits left before integer-dividing by the denominator. We need to make k sufficiently large to catch a repeat. We note that the binary expansion n/d has period at most d, so the joint expansions have period at most d*D, so k=d*D suffices.
The rest is to check whether the fraction is in the box and insulate against inputs given like -1/-3.

source: Stack Exchange

license: Creative Commons - Attribution Share-Alike - https://creativecommons.org/licenses/by-sa/4.0/

all_licenses: Creative Commons - Attribution Share-Alike - https://creativecommons.org/licenses/by-sa/3.0/

sort: votes


text:
Approximate a given image using only colors from a given list
The Specs

Your program must read the image from one file and write the results to another file, the name of at least the first of which should be taken from the input
You may accept and return any image format or formats, as long as you specify in your answers which ones you are using
An exception should be thrown if the specified source file does not exist
The final image must be as close as possible to the initial image (by per-pixel euclidean metric on RGB values) while using only colors from a list read from the file "colors.txt" (it will not necessarily use all of them)
When determing the nearest color, you should use the euclidean metric on RGB values, Δc=sqrt((Δr)2+(Δg)2+(Δb)2), when two or more colors on the list are equally close, any of those colors is valid
The list of colors given will be formatted with one rgb color code on each line, with spaces seperating the red, green, and blue values
You may assume that at least three but no more than 1,000 colors will be listed

This is code golf, so the shortest program wins
Some test cases:
(if RGB values in the output images do not exactly match those in the list, this is due to a problem with Java's ImageIO.write)
input:

colors.txt:

255 0 0
250 150 0
255 50 0
0 0 200
0 0 80
0 0 0

output:

input:

colors.txt:

100 60 0
0 0 20
0 0 70
0 20 90
15 15 15
100 100 100
120 120 120
150 150 150

output:

input:

colors.txt:

0 0 0
10 20 15
25 25 20
35 30 20
40 35 25
55 55 60
70 60 50
130 100 80
130 140 145
140 100 70
165 130 90
200 160 130
230 210 220
250 250 250

output:

input:

colors.txt:

0 0 0
70 50 40
80 0 0
120 100 90
120 120 120
170 150 90
170 170 170
170 200 195
190 140 100
190 205 165
210 170 120

output:

Again, I put together a java program to generate the sample output. I'll golf it down a little further, and then I may post it as an answer later.

Is there any reason to limit file formats to JPG and requiring I/O to be through files, other than punishing/eliminating languages without image and file system manipulation libraries? Why not allow any true-colour raster format (so that one could use PPM for instance), and I/O as byte streams through STDIN/STDOUT?

@MartinBüttner: Is that better?

Yeah much better, although you still require I/O to be through files.

Per-pixel, it's supposed to use the nearest color in the list by the euclidean metric on RGB values.

All looks clear now. The only other potential ambiguity is if there are more than one permitted colour equally close to the target colour (for a given pixel). Do you want to specify a way of deciding in such cases so all outputs are identical, or do you want to specify that in such cases any closest choice is valid?

Note that this is called posterization or color quantization.

what winning measure is taken here . minimal number of colors ? are we allowed to compress size ? do we use DCT ? DWT ? FFT? can we code output using huffman coding?

Why do you have to choose the closest through RGB distance, and not L*a*b*? It would've made for a more interesting golf and the output would actually be correct (RGB is just incorrect for this).

@SuperJedi224 You should create another challenge where the goal is to take the image and a value n as input, and output the n colors that can most closely approximate the image with the lowest difference, possible with a tie-breaker for colors which equally approximate.

@mbomb007 I did something along those lines (as a popularity contest/ code golf hybrid for scoring) and it was closed as a duplicate of a paint by number generation challenge. For some reason.

@Agawa001: Code size. The list of colors will be read from a file.

Mathematica 109 100 97 111
Update: now outputs the image as required.
The function, f, replaces a single pixel with the closest pixel found in the input file (in this case, "colors.txt").
t is the output file.
The function, g, applies f to each pixel in the input image.
g[i_,h_,t_]:=Export[t,(c=Import[h,"Data"]/255;
f@p_:=Sort[{p~EuclideanDistance~#,#}&/@c][[1,2]];f~ImageApply~i)]

Test Cases
Example 1: Seascape
(*creating "colors.txt"*)
    c = "255 0 0
      250 150 0
      255 50 0
      0 0 200
      0 0 80
      0 0 0";
    Export["colors.txt", c]

Example 2: Earth from moon
I'll skip the part about storing the color info into "colors.txt".
It closely follows example 1.

Example 3: Jupiter

Example 4: American Gothic

Python [3] + Pillow, 276 272 267
Takes 2 parameters: name of source image file and name of file with colors list. Tested on Py3, may work on Py2 too.
My first attempt in code golf, so probably not very impressive (and similar to dieter's code).
import sys
import PIL.Image as P
i,c=sys.argv[1:]
c,z=[tuple(int(i)for i in e.split())for e in open(c).readlines()],P.open(i)
x,y=z.size
p=z.load()
for _ in range(x*y):
    j,k=_//y,_%y
    d={sum((i-j)**2for i,j in zip(p[j,k],a)):a for a in c};p[j,k]=d[min(d)]
z.save(i*2)

Ungolfed:
(displays output image instead of saving, and is a bit safer, because removes empty lines from colors file)
# coding: utf-8

import sys
from PIL import Image

def similarImage(img, colors):
    if isinstance(img, str):
        img = Image.open(img)
    x, y = img.size
    pixels = img.load()
    for i in range(x):
        for j in range(y):
            pixels[i,j] = getClosestColor(pixels[i,j], colors)
    img.show()

def getColorDistanceSq(c1, c2):
    #root not needed
    return sum((i-j)**2 for i,j in zip(c1, c2))

def getClosestColor(col, colBase):
    colDist = {getColorDistanceSq(col, c) : c for c in colBase}
    return  colDist[min(colDist)]

def parseColors(desc):
    desc = desc.split('\n')
    desc = [i for i in desc if i != '']
    return [tuple(int(i) for i in d.split()) for d in desc]

def main():
    imgFile   = sys.argv[1]
    with open(sys.argv[2]) as f:
        colors = parseColors(f.read())
    similarImage(imgFile, colors)

if __name__ == "__main__": main()

You are not allowed to overwrite the source image.

Python 2.7, 279
Reads the name of the input image and the name of the color map for commandline arguments. usage is
python script.py <image file> <colors.txt>

It saves the new image as a file with the same name, preceded by an extra tilde :
from sys import*;from PIL.Image import*;f=argv[1]
c,I,X=[[int(x)for x in l.split()]for l in file(argv[2])],open(f),range
w,h=I.size;P=I.load()
for j in X(h):
 for i in X(w):
  r,g,b=P[i,j];C={(R,G,B):(R-r)**2+(G-g)**2+(B-b)**2for R,G,B in c};P[i,j]=min(C,key=C.get)
I.save('~'+f)

I had to find and install PIL first, but it wound up working.

Java, 726 bytes
Too bad all the imports and main class eat up some many bytes, but it's a standalone program after all. Nevertheless, it was a nice and fun challenge to work on.
This golf can be run by issuing java M.class inputImage.png outputImage.png colours.txt, or equivalent.
import java.awt.*;import java.awt.image.*;import java.io.*;import java.nio.file.*;import javax.imageio.*;class M {public static void main(String[]a)throws Exception{BufferedImage i=ImageIO.read(new File(a[0]));int w,h,k,r,g,b,t,z,q;w=i.getWidth();h=i.getHeight();BufferedImage o=new BufferedImage(w,h,2);for(int x=0;x<w;x++){for(int y=0;y<h;y++){Color v=null;k=i.getRGB(x,y);r=k>>16&0xff;g=k>>8&0xff;b=k&0xff;double e=999;for(String n:Files.readAllLines(Paths.get(a[2]))){String[]m=n.split(" ");t=Integer.parseInt(m[0]);z=Integer.parseInt(m[1]);q=Integer.parseInt(m[2]);double d=Math.sqrt((t-r)*(t-r)+(z-g)*(z-g)+(q-b)*(q-b));if(d<e){e=d;v=new Color(t,z,q);}}o.setRGB(x,y,v.getRGB());}}ImageIO.write(o,"png",new File(a[1]));}}

Actually, I believe the .class extension can be omitted for the java command prompt instruction.

if you want to golf it a little further, use a static block and compile it with java 6

Java, 861 849 842 836 bytes
I think there's probably room for improvement here.
import java.awt.*;import java.awt.image.*;import java.io.*;import java.util.*;import javax.imageio.*;import static java.lang.Math.*;class D{static void g(Color c,ArrayList<Color>h,int x,int y,BufferedImage b){double d=450;Color e=null;for(Color f:h){double g=sqrt(pow(c.getRed()-f.getRed(),2)+pow(c.getGreen()-f.getGreen(),2)+pow(c.getBlue()-f.getBlue(),2));if(g<d){e=f;d=g;}}b.setRGB(x, y, e.getRGB());}public static void main(String[]a)throws Exception{ArrayList<Color>c=new ArrayList<Color>();Scanner b=new Scanner(new File("colors.txt"));while(b.hasNext())c.add(new Color(b.nextInt(),b.nextInt(),b.nextInt()));b=new Scanner(System.in);BufferedImage i=ImageIO.read(new File(b.next()));int x=i.getWidth(),y=i.getHeight();for(;--x>=0;)for(int z=0;z<y;z++)g(new Color(i.getRGB(x,z)),c,x,z,i);ImageIO.write(i,"jpg",new File(b.next()));}}

Takes the names of the source and target image files from user input, reads the color list from colors.txt. Supports .jpg, .png, and probably others for input. Supports only .jpg for output.

license: Creative Commons - Attribution Share-Alike - https://creativecommons.org/licenses/by-sa/3.0/


text:
Me Want Honeycomb
Write the shortest program that prints this ASCII art section of a hexagonal tiling or honeycomb:
       __
    __/  \__
 __/  \__/  \__
/  \__/  \__/  \
\__/  \__/  \__/
/  \__/  \__/  \
\__/  \__/  \__/
/  \__/  \__/  \
\__/  \__/  \__/
   \__/  \__/
      \__/

No input should be taken.
Output to stdout or your language's closest alternative.
Instead of a program, you may write a named function that takes no parameters and prints the result normally or returns it as a string.
The output may have any number of leading and/or trailing newlines and each line in the output may have any number of leading and/or trailing spaces (as long as the pattern lines up properly.)
The shortest code in bytes wins.

Do you play a lot of Catan?

@user3490 In fact I do ^^

Next week on PCG: write a Catan board generator?

It should have an input IMO, in your case it should be 3.

@user3819867 I did consider that but preferred this. It's too late to change but not too late for anyone to make a related challenge.

CJam, 45 43 42 41 40 bytes
741e8 36+Ab"\__/  "38*21/.{G2$-<\S*.e<N}

Try it online in the CJam interpreter.
How it works
"\__/  "38*21/

repeats the pattern \__/ 38 times and splits it into chunks of length 21. If the chunks were separated by linefeeds, this would be the result:
\__/  \__/  \__/  \__
/  \__/  \__/  \__/
\__/  \__/  \__/  \__
/  \__/  \__/  \__/
\__/  \__/  \__/  \__
/  \__/  \__/  \__/
\__/  \__/  \__/  \__
/  \__/  \__/  \__/
\__/  \__/  \__/  \__
/  \__/  \__/  \__/
\__/  \__/  \__/

This clearly contains the desired honeycomb. All that's left to do is to replace some characters with spaces, cut off some others and actually introduce the linefeeds.
741e8 36+Ab

generates the integer 74 100 000 036 and converts it into the array [7 4 1 0 0 0 0 0 0 3 6]. Each array element encodes the number of leading characters of the corresponding line that have to get replaced with spaces. By subtracting this number from 16, we also obtain the correct length for this line.
.{            } e# Execute for each pair of a digit D and a line L:
  G2$-<         e#   Chop off L after 16 - D characters.
       \S*      e#   Push a string of D spaces.
          .e<   e#   Compute the vectorized minimum.
             N  e#   Push a linefeed.

Since a space has a lower code point than the other characters of L and vectorized operators leave the characters of the longer string that do not correspond to one of the shorter one untouched, .e< replaces the first D characters with spaces.

This is brilliant. Nice work.

Python 2, 73
i=0
exec"k=max(7-i,i-24,0);print' '*k+('\__/  '*9)[i:][k:16-k];i+=3;"*11

See it run.
For each of the 11 output rows, computes the number of leading spaces k as a maximum of three linear functions that form the envelope of the left side of the hexagon. Because the diagonal lines have slope 3 and -3, it's better to index the row number as i=0,3,...30.
To create the hexagon mesh, we first tile enough of the unit '\__/  '. Then, shifting [i:] realigns it by 3 for odd rows. Finally, we take the needed portion [k:16-k] of it, leaving a margin of k on the left and right.

CJam, 65 56 55 bytes
"Ý6TNð*¯5"303b4b["/\_")"_ "4*S]f=sB/z{_W%"\/"_W%erN}/

Try it online in the CJam interpreter.
Idea
The right half of each line is a reversed copy of the left half with slashes and backslashes swapped. Therefore, it suffices to encode the left half of the honeycomb:
       _
    __/
 __/  \_
/  \__/
\__/  \_
/  \__/
\__/  \_
/  \__/
\__/  \_
   \__/
      \_

Rather than analyzing this pattern line by line, we can analyze it column by column:
   /\/\/\
  _ _ _ _
  _ _ _ _
  /\/\/\/\
 _ _ _ _ _
 _ _ _ _ _
 /\/\/\/\/\
_ _ _ _ _ _

Obvious patterns emerge:

The string _ _ _ _ occurs five times.
Every / is followed by a \.

By replacing every /\, _, _ _ _ _ and space with a number from 0 to 3, we can convert the resulting array from a base 4 number to a higher base and store the complete pattern in a compact fashion.
Code
"Ý6TNð*¯5"303b4b  e# Convert the string from base 303 to base 4.
["/\_")"_ "4*S]f= e# Replace each resulting digit by the corresponding item of the array
                  e# ["/\" "_" "_ _ _ _ " " "].
sB/               e# Split into strings of length 11.
z                 e# Zip: transpose rows with columns.
{             }/  e# For each string:
 _W%              e#     Push a reversed copy.
    "\/"_W%er     e#     Swap slashes and backslashes.
             N    e#     Push a linefeed.

C, 148 144 140 bytes
k,r,c,d,p;f(){for(;k<187;k++){r=k/17;c=k%17;d=c+r%2*3;p=3*r+c-7<33u&3*r-c+8<33u;putchar(c==16?10:p&(d+5)%6<2?95:p&d%6==3?47:!p|d%6?32:92);}}

With whitespace, without compiler warnings, and before some code tweaks to save a few bytes:
#include <stdio.h>

int k, r, c, d, p;

void f() {
    for ( ; k < 187; k++) {
        r = k / 17;
        c = k % 17;
        d = c + 3 * (r % 2);
        p = 3 * r + c - 7 < 33u && 3 * r - c + 8 < 33u;
        putchar(
            c == 16 ? 10 :
            p && (d + 5) % 6 < 2 ? 95 :
            p && d % 6 == 3 ? 47 :
            p && d % 6 == 0 ? 92 :
            32);
    }
}

This approach does not use any character/string tables. It loops over all 187 (11 rows, 17 columns including newlines), and decides which character to print for each position, based on a combination of conditions.
The conditions include a test for being inside/outside the 4 corners, using 4 line equations, with the result stored in variable p. The rest is then mostly repeating every 6 characters, with the odd rows shifted by 3 characters relative to the even rows.

For the golfed one, you could use implicit-int and drop the int .

A few micro-improvements: k;f(r,c,d,p){for(;k<187;putchar(c>15?10:p&d<2&p?95:p&d<3?47:!p|d-5?32:92))r=k/17,c=k++%17,d=(5+c+r%2*3)%6,p=3*r+c-7<33u&3*r-c+8<33u;}

Ruby - 78 Bytes
(-5..5).map{|i|puts' '*(j=[i.abs*3-i/6-9,0].max)+('/  \__'*4)[j+i%2*3,16-j*2]}

A transcription of xnor's solution (69 Bytes):
11.times{|i|puts' '*(j=[7-i*=3,i-24,0].max)+('\__/  '*9)[i+j,16-j*2]}

11.times saves 2 bytes

JavaScript (ES6), 129 130
That's a pure string replace/replace/replace ... not taking advantage of any geometric property.
Using template strings, all newlines are significant and count.
Run snippet in Firefox to test

f=_=>`70
405
 055
9992 3051
6301`[R='replace'](/9/g,`55123
30551
`)[R](/5/g,1230)[R](/\d/g,n=>['__','/',,'\\'][n]||' '.repeat(n))

// TEST

O.innerHTML = f()

PHP - 139 137 107 101 91 87 bytes
I don't know if this is the best way to  golf it but here is my try at it:
30 36 46 -50 bytes thanks to Ismael Miguel
Test it online here
<?='       __
    __',$a='/  \__',"
 __$a$a
",$c="$a$a/  \
\__$a$a/
","$c$c   \__$a/
      \__/";

old code:
<?php $a="/  \\__";$c=$a.$a."/  \\\n\__".$a.$a."/\n";echo "       __\n    __".$a."\n __".$a,$a."\n".$c,$c,$c."   \\__".$a."/\n      \\__/";

Welcome to Programming Puzzles and Code Golf. Your PHP code can be reduced a lot more. You can use <? instead of  <?php to save 3 bytes. Instead of $c=$a.$a."/  \\\n\__".$a.$a."/\n";, you can write $c="$a$a/  \\\n\__.$a$a/\n"; (since PHP expands variables in strings). You can apply the same logic to the echo to further reduce it's length. Also, you don't need that space between the echo and the string.

Thank you, I learn something new every time I try to golf hehe.

You are welcome. Here is a 99-byte long solution: https://ideone.com/WHWEZS. It looks a mess, but it is a lot shorter. Instead of using \n, you can put a real  newline and save 1 byte per line.

Meanwhile, I've reduced it 2 bytes. Now it has 97 bytes.

Playing around, I've reduced that file to 88 bytes. On it, one thing to look out is that I've replaced each group of 4 spaced by a tab (\t). This means that you slim down 3 bytes. I will paste an hex dump of the file (or similar), just to be sure that everybody has it right and that the browser doesn't replace the tabs with 4 spaces.

Here is a one-liner: <?="\t   __\n\t__",$a='/  \__',"\n __$a$a\n",$c="$a$a/  \\\n\__$a$a/\n","$c$c   \__$a/\n\t  \__/";. Remember to replace the \t with a tab character and \n with a real newline.

For me the \t is 4 spaces too long, I'm running it in the Windows Command Promt.

Tha't alright. As long as you can reproduce it. That means, if it works under a certain environment, and someone else can pick up and reproduce the same results, then it's all good. You can check the ideone link (https://ideone.com/WHWEZS) showing it working fine. It's the same code as before, but with the tabs and newlines replaced. As an alternative, you can click on the icon that looks like <> and paste this: .

Lua 146
T="       __\n    __/  \\__\n __/  \\__/  \\__\n"
S="/  \\__/  \\__/  \\\n\\__/  \\__/  \\__/\n"
B="   \\__/  \\__/\n      \\__/"
io.write(T,S,S,S,B)

(newlines added for clarity)

Isn't this code longer than the desired output because it just contains the output plus escaping and string handling stuff?

@Caleb: I suspect you've not actually counted & just guessed in making that comment. According to my text editor, there are 165 characters in the lattice (newlines included). I have 19 less characters due to repeating S three times.

But you are correct that the escapes and newlines do limit my code's ability to compete in this particular game. But winning also isn't why I do this, it's for the fun & learning.

Dart - 113
main({r:"/  \\__"}){print("       __\n    __$r\n __$r$r\n${"$r$r/  \\\n\\__$r$r/\n"*3}   \\__$r/\n      \\__/");}

Pure string-interpolation solution, nothing fancy.
String operations like "substring" are too verbose to compete in practice.
Run it on DartPad.

Python 3, 100 87 bytes
n=0x85208e08e08dd445
while n:l=n>>2&15;print(' '*(8-l)+('\__/  '*4)[n&3:n&3|l*2]);n>>=6

And a readable version of the code below.
The idea is to hardcode the intervals (begin, length) and then pad with the correct amount of spaces to center the interval.
hexagon_string = '\__/  ' * 4
for i in range(11):
    # Pick the begin and length of the i-th interval using hardcoded numbers
    b = (0x00030303111 >> 4*i) & 3   # n is a combination of these two numbers
    l = (0x25888888741 >> 4*i) & 15  #
    # Print the interval with appropriate whitespace in front of it
    spaces = ' ' * (8-l)
    print(spaces + hexagon_string[b : b+l*2])

Javascript (ES7 Draft), 96 94 93 bytes
Inspiration taken from a few solutions here...
Edit: -1 from edc65

f=z=>[for(x of[741e6]+36)' '.repeat(x)+'/  \\__'.repeat(4,z^=3).slice(+x+z,16-x+z)].join(`
`)

// For snippet demo:
document.write('
'+f()+'
');

Commented:
f=z=>[
        for(x of [741e6] + 36)        // for each character "x" in "74100000036"
            ' '.repeat(x) +           // x spaces
            '/  \\__'.repeat(4,       // build pattern string
            z ^= 3).                  // toggle z offset (3) for even lines
            slice(+x + z, 16 - x + z) // get appropriate substring
    ].join(`                          // join with newlines
    `)

Nice one. .substr(+x+z,16-x-x) -> .slice(+x+z,16-x+z) -1

@edc65 nice catch!

Quick question: Why is it +x instead of just x?

@QPaysTaxes It's to cast x to a number. Without it, it would concatenate z

Retina, 66 bytes
<empty line>
sss __<LF>ssa__<LF> aa__<LF>lll dau<LF>ssdu
l
/daa<LF>\aau<LF>
a
ud
u
__/
d
s\
s

Each line should go to its own file and <LF> means actual newline in the file. 1 byte per extra file added to byte-count.
You can run the code as one file with the -s flag, keeping the <LF> markers and maybe changing them to newlines in the output for readability if you wish.
The algorithm is 5 simple substitute steps (change odd line content to even line content) starting from an empty input string. The results after each steps are (delimited by ='s):
sss __
ssa__
 aa__
lll dau
ssdu
============
sss __
ssa__
 aa__
/daa
\aau
/daa
\aau
/daa
\aau
 dau
ssdu
============
sss __
ssud__
 udud__
/dudud
\ududu
/dudud
\ududu
/dudud
\ududu
 dudu
ssdu
============
sss __
ss__/d__
 __/d__/d__
/d__/d__/d
\__/d__/d__/
/d__/d__/d
\__/d__/d__/
/d__/d__/d
\__/d__/d__/
 d__/d__/
ssd__/
============
sss __
ss__/s\__
 __/s\__/s\__
/s\__/s\__/s\
\__/s\__/s\__/
/s\__/s\__/s\
\__/s\__/s\__/
/s\__/s\__/s\
\__/s\__/s\__/
 s\__/s\__/
sss\__/
============
       __
    __/  \__
 __/  \__/  \__
/  \__/  \__/  \
\__/  \__/  \__/
/  \__/  \__/  \
\__/  \__/  \__/
/  \__/  \__/  \
\__/  \__/  \__/
   \__/  \__/
      \__/

Javascript, 154 151 bytes
b="\n",c="\t",d="/  \\",f="__",g="\\",h="/",a=f+d,j=a+a,i=d+j+b+g+j+f+h+b,k="   ";console.log(c+k+f+b+c+a+f+b+" "+j+f+b+i+i+i+k+g+a+f+h+b+c+"  "+g+f+h)


text:
Thorough though Thoreau threw, Troll throws through tough troughs
To each of these nine confusingly similar words, assign a number 1-9 in any way you like:
though
through
thorough
Thoreau
throw
threw
trough
tough
troll

Write a program that takes in a string. If the input is one of these nine words, output the number you assigned to it. If the input is not one of the words above, the program may do anything (including error or loop forever).
The words are case sensitive, e.g. Thoreau, should produce a number from 1-9 but thoreau will not necessarily do the same.
Example
Suppose you assign the numbers as follows:
though   9
through  2
thorough 7
Thoreau  6
throw    3
threw    5
trough   4
tough    1
troll    8

Then when tough is input, 1 should be output.
When through is input, 2 should be output.
When throw is input, 3 should be output.
. . .
When though is input, 9 should be output.
All other inputs may do anything.
Details

Take the input string via stdin or the command line and output to stdout.
The output may contain a single trailing newline.
Instead of a program, you may write a function that takes in a string and prints the result normally or returns it.
The shortest submission in bytes wins.

Darn it! I had a clever solution to output zero when not found by using the Python string find method. Then the rules changed. Clever idea not so clever now.

@CarpetPython My bad really. Don't hesitate to downvote if you feel unsatisfied with the change. (Though I promise to everyone there will be no more changes.)

That's ok. I think my answer is still valid (though I little verbose).

Can I make it work regardless of capitalization?

@ASCIIThenANSI as long as it works for the 9 cases

CJam, 11 9 7 bytes
q1b2+B%

How it works:
We are making use of the fact that the sum of the ASCII codes + 2 moded with 11 gives very nice order of 1 through 9 and then 10 for the nine concerned words. Here is the ordering:
through -> 1
thorough -> 2
tough -> 3
Thoreau -> 4
throw -> 5
threw -> 6
trough -> 7
though -> 8
troll -> 9

Code explanation:
q               e# Read the input
 1b             e# Sum the ASCII code values of all characters in this word
   2+           e# Increment the sum by 2
     B%         e# Mod by 11 and automatically print the mod result at the end

4 bytes saved thanks to user23013
Try it online here

How are we supposed to try it? (not everyone speaks CJam, some of us speak lolcode)

@Behrooz There is a link. Click on it, put your input in the input section, Run. Not sure how it can be easier :)

Holly shit, I was thinking how am I supposed to give it the list of words. nice one

:i:+ is the same as 1b.

@user23013 Damn! Every single time!

In fact it's probably because I'm currently not posting so many answers. And if my answer is longer (sometimes much longer) than existing answers, I choose not to post in most cases.

Looking at the link: it should work, but it does not. Just a CJAM empty box, no code, no input. In Firefox.

@edc65 some known issues with firefox. Try opening the link in chrome

Or q1b2+B%.

ÇOÌ11% - 05AB1E, your method, 6 bytes.

Pyth, 8 chars
e%Cz8109

Try it online: Demonstration or Test Suite
I'm using the assignment:
though   5
through  9
thorough 4
Thoreau  7
throw    3
threw    2
trough   8
tough    6
troll    1

Explanation:
   z       input()
  C        convert to int (convert each char to their ASCII value
           and interprete the result as number in base 256)
 %  8109   modulo 8109
e          modulo 10

Btw, I found the magic number 8109 by using this script: fqr1 10Sme%CdT.z1.

Will this not have false-positives?

@alexander-brett What exactly do you mean? The output of all other inputs is not specified in the OP. We can output anything we want.

Sorry, I missed that edit to the OP. That's a shame :P

Python 2, 92 54 bytes
print'K{7j)yE<}'.find(chr(hash(raw_input())%95+32))+1

The index string is created with for word in words: print chr(hash(word)%95+32),. As pointed out in Jakube's answer, the hash function will give different results depending on Python version. This index string is computed on 64 bit Python 2.7.6.
Longer (92 bytes) but less cryptic answer:
print'though through thorough Thoreau throw threw trough tough troll'.find(raw_input())/7+1

The programs returns 1-9 for though through thorough Thoreau throw threw trough tough troll in that order. When the input is not found, find will return a -1 which conveniently turns into a zero after the +1.

Note that the 0 stuff is no longer required. Sorry for changing it on you.

Thanks for noticing the effort. It was a good solution for a little while...

@CarpetPython Nice use of floor-division -- it works out surprisingly neatly.

Python 2.7.9 32 bit version, 22 bytes
lambda x:hash(x)%78%10

Notice, the version is really important here. You will get different results if your using a 64 bit version of Python. Since the hash method will compute 64 bit hash values instead of 32 bit.
The assignment is:
though  => 5   through => 6   thorough => 8
Thoreau => 7   throw   => 3   threw    => 1
trough  => 9   tough   => 4   troll    => 2

Try it online: http://ideone.com/Rqp9J8

Wow, so all this time, you were iterating through language versions and operating system bits ? :P

Very nice answer. Did you find the constant 78 through mathematics, automated iteration, or some guesses?

@CarpetPython Just a simple brute-force loop which goes through all possible modules. Once sorted(...)==range(1,10), I stopped.

Pyth, 7 bytes
et%Cz31

I am using the following asignment:
though   8
through  3
thorough 1
Thoreau  5
throw    4
threw    7
trough   6
tough    2
troll    9

Cz interprets the input as a base 256 number. Then, we take this mod 31, subtract 1, and take the result mod 10. Equivalent pseudocode:
((base_256(input()) % 31) - 1) % 10

Demonstration, test harness.

Japt, 6 bytes
nH %BÉ

Try it | Test all words

Explanation
Takes advantage of the fact that, when parsing a base-n string to an integer, JavaScript will stop parsing if it encounters a digit greater than n and return the result up to that point. By using base-32 here (digits 0-v) the ws in "threw" and "throw" are, essentially, ignored.
nH      :Convert from base-32
   %B   :Modulo 11
     É  :Subtract 1

Here's a breakdown of the formula applied to each word:

Input
nH
%B
É

through
31738067473
2
1

tough
31226385
3
2

troll
31318709
4
3

trough
1002207761
5
4

though
991722001
6
5

throw
968568
7
6

threw
968558
8
7

thorough
1015520459281
9
8

Thoreau
31735003486
10
9

JavaScript, 22 bytes
A direct translation - doesn't seem worth posting it separately.

f=
U=>parseInt(U,32)%11-1
o.innerText=["through","tough","troll","trough","though","throw","threw","thorough","Thoreau"].map(s=>f(s)+": "+s).join`\n`

Vyxal, 6 bytes
C∑⇧11%

Try it Online!
Cjam port.

Python 2, 27 bytes
f=lambda w:int(w,34)%444/46

With this assignment:
>>> for w in "though through thorough Thoreau throw threw trough tough troll".split(): print f(w),w
...
9 though
7 through
3 thorough
8 Thoreau
2 throw
5 threw
6 trough
1 tough
4 troll

Several variations are possible, e.g.
f=lambda w:int(w,35)/159%10

C (gcc), 66 bytes
h,k;f(char*s){for(h=33;*s;)h^=*s++;h=strchr(k="(Z5qW]2@H",h)-k+1;}

Try it online!

May require the -O compiler flag. h;f(char*s){for(h=33;*s;)h^=*s++;h=index("(Z5qW]2@H",h)-"H"+9;}

Java 8, 53 25 bytes
s->(s.chars().sum()+2)%11

or
s->-~-~s.chars().sum()%11

Port of @Optimizer's CJam answer, because it (most likely) cannot be done any shorter in Java..
Try it online.

Java has parseInt, doesn't it? Would a port of my solution not be shorter?

@Shaggy Java indeed has parseInt with given base, but unfortunately it's quite byte-excessive due to static class requirement: Long.parseLong(...,32) as shortest variant. In addition, it seems to fail for "throw" (and "threw" as well) in Java for some reason. w is outside the base-32 range it seems (and using 33 gives incorrect results).

Jelly, 7 bytes
OS+2%11

Try it online!
Boring Jelly port of the fantastic CJam answer.

all_licenses: Creative Commons - Attribution Share-Alike - https://creativecommons.org/licenses/by-sa/4.0/


text:
Given an input, move it along the keyboard by N characters
The Challenge:
Given any input that can be typed on a keyboard, move the text along by N chars.
Here is the QWERTY keyboard to be used. You can ignore the modifier keys (Shift, Caps, Enter, Delete and Tab). Once you reach one side (for example |) loop back round, so | goes to Q if N = 1.

Spaces do not have to be moved along (they go back round to space as you skip modifiers). If shift was used to type the character (e.g. ! and @) the changed character should also be typed using shift (i.e. ! goes to @ not 2 if N = 1).
UK keyboards are different to this, but please use this so that we can compare.
Input:
Any sentence that can be typed on the above keyboard followed by a positive integer. There is no maximum to the size of this integer.
Output:
The same sentence, shifted along by N.
Examples:
My name is Tim 3
?o .f/y [g I[/

Hello World 7
Spgge Oe[g;

I Wi5h I h4d b3773r C@d3ing ski{{s 3
{ T[8l { l7h ,6006u N%h6[.k g'[QQg

This is code golf, so shortest code wins.

Can we assume that N<= 13 as you will need to shift at most 13 in order to return to the original character?

@flawr nope, sorry. It can be any positive value.

Shouldn't the "Hello World 7" example be "Spggr Oe[g;"?

Shouldn't the Hello World 7 example be Spgge Oe[g;?  The two o should map to the same character

C, 217 bytes
char*t=" @A$%^*a)_(~.=/z-234567890\"'>`?Z#SNVFRGHJOKL:<MP{WTDYIBECUX]q\\&=1snvfrghjokl;,mp[wtdyibecux}Q|!";l,d,k;f(char*s){for(l=strlen(s);s[--l]-32;);d=atoi(s+l);for(s[l]=0;d--;)for(k=l;k--;s[k]=t[s[k]-32]);puts(s);}

Readable version with whitespace, includes, etc:
#include <string.h>
#include <stdlib.h>
#include <stdio.h>

char* t = " @A$%^*a)_(~.=/z-234567890\"'>`?Z#SNVFRGHJOKL:<MP{WTDYIBECUX]q\\&=1snvfrghjokl;,mp[wtdyibecux}Q|!";
int l, d, k;

void f(char* s) {
    l = strlen(s);
    for( ; s[--l] - 32; );
    d = atoi(s + l);
    s[l] = 0;
    for ( ; d--; ) {
        for (k = l; k--; s[k] = t[s[k] - 32]);
    }
    puts(s);
}

The code pretty much speaks for itself. Just a lookup table that maps from each character to the next character, which is applied the given number of times. Much of the code is actually for parsing the number out of the input.

@Ypnypn You can use undeclared functions in C. So the includes are not needed for it to build. It typically gives compiler warnings, but I've been told that this is allowed as long as it builds and runs.

05AB1E, 61 bytes
-1 byte thanks to Kevin Cruijssen
žVDu«1ú`žhÀ“~!@#$%^&*()_+ `ÿ-=ÿ<>?ÿ:"ÿ{}|ÿ,./ÿ;'ÿ[]\“#vyy¹._‡

Try it online!

Pyth, 126 bytes
XjdPczdsJc"~!@#$%^&*()_+ `1234567890-= qwertyuiop[]\ QWERTYUIOP{}| asdfghjkl;, ASDFGHJKL:\" zxcvbnm,./ ZXCVBNM<>?")sm.<dvecz)J

Try it online: Demonstration or Test Suite
Explanation:
    czd                       split input by spaces
   P                          remove the last element
 jd                           join by spaces (=#1)

          "..."               string with the chars of each row
         c     )              split by spaces
        J                     assign to J
       sJ                     sum of J (=#2)

                       cz)    split input by spaces
                      e       take the last element
                     v        and evaluate it
                 m        J   map each row d of J to:
                  .<d           rotate the row d by value
                s             sum (=#3)

X                             Take #1, and replace the chars in #2 by the chars in #3

Python 3, 311 bytes
*i,s=input().split()
r=["`1234567890-=","qwertyuiop[]\\","asdfghjkl;'","zxcvbnm,./","~!@#$%^&*()_+","QWERTYUIOP{}|",'ASDFGHJKL:"',"ZXCVBNM<>?"]
print("".join([[x[int(s):]+x[:int(s)]for x in r][r.index([x for x in r if c in x][0])][([x for x in r if c in x][0]).index(c)]if c!=" "else " " for c in " ".join(i)]))

Remove the unnecessary spaces in " " for c in " "

Python 3, 271 255 bytes
Baseline, almost ungolfed, used to create the shifted words in the question.
x=input().split()
n=int(x[-1])
x=' '.join(x[:-1])
l=['`1234567890-=','qwertyuiop[]\\',"asdfghjkl;'",'zxcvbnm,./', '~!@#$%^&*()_+','QWERTYUIOP{}|','ASDFGHJKL:"','ZXCVBNM<>?',' ']
y=''
for i in x:
 for q in l:
  if i in q:y+=q[(q.index(i)+n)%len(q)]
print(y)

Explanation:
x=input().split()    # Get input
n=int(x[-1])         # Get N from input
x=' '.join(x[:-1])   # Get the words from input
                     # Create list of letters

l=['`1234567890-=', 'qwertyuiop[]\\',
   "asdfghjkl;'",   'zxcvbnm,./',
   '~!@#$%^&*()_+', 'QWERTYUIOP{}|',
   'ASDFGHJKL:"',   'ZXCVBNM<>?',
   ' ']

y=''                 # Blank string
for i in x:          # Loop through letters in input
    for q in l:      # Loop through items in list
        if i in q:   # Is letter of input in item of list?
            y+=q[                          # Append letter to y
                 (q.index(i)+n)            # locate the letter in item, and add N
                               %len(q)]    # % is modulus, loop to beginning if big
print(y)             # Print out the offset word.

I think you should delete this. Let other people work out their own strategies...

@mbomb007 it's not very golfed, and I used it to create them... I think it's fair enough to post it, personally.

JavaScript (ES6), 200 216
Using template strings, the newlines are significant and counted.
Note about replace: the two snippets
string.split('x').map(w=>...) and
string.replace(/[^x]+/g,w=>...) are equally valid ways to execute a function for each part in a string, using a separator. Using a newline as a separator is handy as the replace regexp becomes /.+/g, because the dot match any non-newline. And using templated strings the newlines have no extra cost.

f=(t,d)=>[for(c of t)`~!@#$%^&*()_+
1234567890-=
QWERTYUIOP{}|
qwertyuiop[]\\
ASDFGHJKL:"
asdfghjkl;'
ZXCVBNM<>?
zxcvbnm,./`.replace(/.+/g,r=>(p=r.indexOf(c))<0?0:q=r[(p+d)%r.length],q=c)&&q].join('')

// less golfed
x=(t,d)=>
  [for(c of t)
    '~!@#$%^&*()_+ 1234567890-= QWERTYUIOP{}| qwertyuiop[]\\ ASDFGHJKL:" asdfghjkl;\' ZXCVBNM<>? zxcvbnm,./'
    .split(' ')
    .map(r=>(p=r.indexOf(c))<0?0:q=r[(p+d)%r.length],q=c)&&q
  ].join('')

// TEST

out=x=>O.innerHTML+=x+'\n'

;[['Hello World',7,],['My name is Tim',3],['I Wi5h I h4d b3773r C@d3ing ski{{s', 3]]
.forEach(p=>out(p+' -> '+f(p[0],p[1])))

CJam, 107 bytes
lS/)~\S*\",./ ;'  <>? :\"  _+~!@#$%^&*() -=`"A,s(++S/"zxcvbnm
asdfghjkl
[]\qwertyuiop"N/_32ff^+.+_@fm>s\ser

Try it online in the CJam interpreter.
How it works
lS/)   e# Read one line from STDIN, split at spaces and pop the last chunk.
~\S*\  e# Evaluate the popped chunk and join the remaining ones back together.
",./ ;'  <>? :\"  _+~!@#$%^&*() -=`"
       e# Push that string.
A,s(++ e# Concatenate it with "1234567890".
S/     e# Split at spaces.
"zxcvbnm asdfghjkl []\qwertyuiop"
       e# Push that string.
S/     e# Split at spaces. (`N/' would split at linefeeds.)
_32ff^ e# XOR each character of a copy with 32.
+      e# Concatenate the copies.
.+     e# Perform vectorized concatenation. This pushes the following array:
          [ ",./zxcvbnm" ";'asdfghjkl" "[]\qwertyuiop" "<>?ZXCVBNM"
           ":\"ASDFGHJKL" "{}|QWERTYUIOP" "_+~!@#$%^&*()" "-=`1234567890" ]
_@fm>  e# Rotate each chunk by the number of character specified in the input.
s\s    e# Flatten this array and the original.
er     e# Perform transliteration.

Jelly, 67 bytes
ØDṙ1ṭØQ;Øq;"“{}|“:"“<>?“-=`“[]\“;'“,./“~!@#$%^&*()_+“ ”
¢œiⱮ+2¦€œị¢

Try it online!
A dyadic link taking the string as its left argument and the number of places to shift as its right argument.

Python 2, 194 bytes
f=lambda s,n:n and f(''.join(B[B.find(c)+1]for c in s),n-1)or s
A='qwertyuiop%sqasdfghjkl%sazxcvbnm%sz'
B='`1234567890-=`~!@#$%^&*()_+~'+A%('[]\\',";'",',./')+(A%('{}|',':"','<>?')).upper()+'  '

Try it online!